python-CAD绘制等高线

基于python的CAD等高线绘制 这是我本科在做《路基工程课程设计》挡土墙设计的时候,老师只给了原始地面数据点,因为疫情在家时间充裕,我就摸索了python画CAD,最终实现了绘制不同等高距的等高线绘制。 python面向对象编程。 基础知识 我这个代码只适合地面趋势简单的,通过差值得到固定整数高程点,然后将相同高程点依次连接。 note: 当时是在windows先写的代码,因为要用到win32com.client这个库,所以我估计 Mac下应该用不了,但可以修改代码绘制dxf文件,不调用CAD。 后续我再补充 Excel文件、输出CAD的文件及效果图 源代码 # _*_ coding:utf-8- _*_ from pyautocad import APoint, Autocad, aDouble import win32com.client import numpy as np dgj = 0.5 # 等高距 ctfs = 10 # 放大倍数,本身采取m为单位绘制,可能放大10倍后输出效果好些 """ 等高线计算 """ #连接excel读取数据,这个地方可以做一个前处理,直接把excel数据作为数组或者numpy数组进行处理 xlapp = win32com.client.Dispatch("Excel.Application") wb = xlapp.ActiveWorkbook ws = xlapp.ActiveWorkbook.ActiveSheet """ K58+70向右走为前进方向,绘制在 """ """ First 定义一个类,输入 list,其元素为APoint 内部函数得到插值的数据点,打印点数据 """ class nei_cha: def __init__(self, pnts = [APoint(0,1)], dgj = 0.5, gu_ding_cha_zhi = "X" ): # gu_ding_cz = "X" 意思是不进行固定插值操作 self.pnts = pnts self.dgj = dgj self.gu_ding_cha_zhi = gu_ding_cha_zhi # 输入两个三维坐标,x/y是相等的,内插得到所需高程 # 返回的是,该区间内得到的点组成的list def do_two_points(self, p1 = APoint(0), p2 = APoint(0)): if self.gu_ding_cha_zhi == "X": nmax = int(max(p1.z, p2.z) // self.dgj) nmin = int(min(p1.z, p2.z) // self.dgj) if nmax == nmin: # 这两个点之间内插不出来 return 0 else: nmax = self.gu_ding_cha_zhi nmin = self.gu_ding_cha_zhi - 1 if (p1.z - nmax) * (p2.z - nmin) > 0: # 固定插值不在区间内 return 0 if p1.x == p2.x: # x坐标一样,需要内插Z到y """ 求直线方程,y """ k = (p1.y - p2.y)/(p1.z - p2.z) b = p1.y - k * p1.z if self.gu_ding_cha_zhi == "X": ni = np.array(range(nmin, nmax)).reshape(int(nmax - nmin), 1) zi = (ni + 1) * self.dgj else: ni = np.array([self.gu_ding_cha_zhi]).reshape(int(nmax - nmin), 1) zi = ni yi = np.around(k * zi + b, 3) #yi = np.interp(zi, [p1.z, p2.z], [p1.y, p2.y]) xi = np.full((len(yi), 1), p1.x) pts_np = np.hstack((xi, yi, zi)) return pts_np else: # y坐标一样,需要内插Z得到x return 0 def do_it(self): pnts_np = np.array([]) for i in range(len(self.pnts) - 1): p_np_i = self.do_two_points(self.pnts[i], self.pnts[i+1]) if type(p_np_i) != type(0): try: pnts_np = np.vstack((pnts_np, p_np_i)) except: pnts_np = p_np_i if len(pnts_np) != len(np.array([])): return pnts_np else: return 0 # points是原始地面点高程 ground_points = [] row = 2 while row <= 20: i = 6 # 左侧数据 j = 8 # 右侧数据 # 桩号点绘图坐标 x0 = ((row - 2) / 3) * 10 + 70 y0 = 0 z0 = ws.cells(row, 7).value p0 = APoint(x0, y0, z0) pts = [] pts.append(p0) # 拾取左侧点数据 while i >= 1: try: z_ = ws.cells(row, i).value y_ = ws.cells(row + 1, i).value pti = pts[-1] + (0, y_, z_) pts.append(pti) except: pass i -= 1 # 拾取右侧点数据 while ws.cells(row, j).value: try: z_ = ws.cells(row, j).value y_ = - ws.cells(row + 1, j).value ptj = pts[0] + (0, y_, z_) pts.insert(0, ptj) except: print("右侧点拾取失败了:cells(%d,%d)"% (row, j)) pass j += 1 ground_points.append(pts) for i in pts: print(i) print("") row = row + 3 # 利用上方点开始使用类进行,插值 all_np = [] for i in ground_points: cha_zhi_la = nei_cha(i, dgj) a = cha_zhi_la.do_it() try: all_np = np.vstack((all_np, a)) except: all_np = a pts_all = [APoint(list(i)) for i in all_np] pts_all.sort(key = lambda x:x.z) del all_np h_min = pts_all[0].z h_max = pts_all[-1].z zzzi = h_min sorted_pts = [] while zzzi <= h_max: _ = [i for i in pts_all if i.z == zzzi] _.sort(key = lambda x:x[0]) sorted_pts.append(_) zzzi += dgj del pts_all for i in sorted_pts: elevation = i[0].z print("高程为 %s 的数据点"% elevation) for l_ in i: print(l_) # 连接cad,并输出结果 del sorted_pts[0:11] del sorted_pts[-10:-1] acad = Autocad(create_if_not_exists=True) acad.prompt("connect with the Autocad successfully.") print("connect with the Autocad successfully.") dwg = acad.ActiveDocument # 定义等高线图层 dgx = acad.doc.layers.add("等高线") dgx.LineWeight = 25 dwg.ActiveLayer = dgx dgxbz = acad.doc.layers.add("等高线标注") for m in range(len(sorted_pts)): i = sorted_pts[m] elevation = i[0].z print("高程为 %s 的数据点"% elevation) for l_ in i: print(l_) try: for j in i: j.z = 0 pt = aDouble([j for _ in i for j in _ * ctfs]) hl = acad.model.AddPolyLine(pt) hl.LineWeight = 25 hl.Elevation = elevation hl.Thickness = 900 if m % 5 == 0: hl.LineWeight = 100 hl.color = 1 text_string = elevation insertpt = i[0] * ctfs textobj = acad.model.addMtext(insertpt, 26, text_string) textobj.Height = 8 textobj.BackgroundFill = True textobj.color = 1 textobj.layer = "等高线标注" except: print('*' * 10, "高程为 %s 的数据点不足,跳过绘制此等高线"% elevation) pass

2022年5月24日 · 3 分钟 · Vinkey

结构可靠度理论代码共享

本文将记录《结构可靠度理论》的一些代码及课堂笔记,采用Matlab编程 1 概率分布密度(pdf、cdf) 自己编写函数 clc clear close all % unif a=1; b=3; x=1:0.1:3; for i=1:length(x) fx(i)=1/(b-a); Fx(i)=(x(i)-a)/(b-a); end h1=figure(1); ax1=axes('Parent',h1,'Position',[0.1,0.4,0.35,0.35],'FontSize',10); plot(x,fx,'r-') xlabel('x','fontsize',10) ylabel('pdf','FontSize',10) ax2=axes('Parent',h1,'Position',[0.55,0.4,0.35,0.35],'FontSize',10); plot(x,Fx,'r-') xlabel('x','fontsize',10) ylabel('cdf','FontSize',10) %------- mux=1;sigx=1; x=-2:0.1:4; syms t for i=1:length(x) tx=x(i); a1=1/(sigx*(2*pi)^0.5); b1=exp(-0.5*((tx-mux)/sigx)^2); fx(i)=a1*b1; c1=int(exp(-0.5*((t-mux)/sigx)^2),t,-inf,tx); Fx(i)=a1*c1; end h1=figure(2); ax1=axes('Parent',h1,'Position',[0.1,0.4,0.35,0.35],'FontSize',10); plot(x,fx,'r-') xlabel('x','fontsize',10) ylabel('pdf','FontSize',10) ax2=axes("Parent",h1,"Position",[0.55,0.4,0.35,0.35],"FontSize",10); plot(x,Fx,'r-') xlabel('x','fontsize',10) ylabel('cdf','FontSize',10) 调用Matlab自带函数 clc clear close all %----- a=1;b=3; x=-2:0.1:5; fx=unifpdf(x,a,b); Fx=unifcdf(x,a,b); figure(1) plot(x,fx,'-r') hold on plot(x,Fx,'--b') %--------- mux=0;sigx=1; x=-3:0.1:3; fx=normpdf(x,mux,sigx); Fx=normcdf(x,mux,sigx); figure(2) plot(x,fx,'r-') hold on plot(x,Fx,'r-') ...

2022年5月14日 · 1 分钟 · Vinkey