From abd0ca8aa78134883fa8fe11ee2d0246253ecd77 Mon Sep 17 00:00:00 2001 From: n3040 Date: Mon, 17 Jan 2022 00:52:02 +0800 Subject: [PATCH] =?UTF-8?q?=E5=A2=9E=E5=8A=A0=E8=BE=93=E5=87=BA=E5=B7=A6?= =?UTF-8?q?=E5=8F=B3=E8=BE=B9=E7=BA=BF=E6=96=AD=E9=9D=A2=E5=8A=9F=E8=83=BD?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- dem.py | 85 ++++++++++++++++++++++++++++++++++++++++------------------ 1 file changed, 59 insertions(+), 26 deletions(-) diff --git a/dem.py b/dem.py index 3d03ffc..68a7035 100644 --- a/dem.py +++ b/dem.py @@ -105,8 +105,8 @@ class Dem: # * delta_y # ) - # 通过空间平面方程乃拟合Z值 - equa_a = np.array( + # 通过空间平面方程乃拟合Z值 参考《数字高程模型教程》 汤国安，李发源，刘学军编著 p89 + equation_a = np.array( [ [lu_xy[0], lu_xy[0] * lu_xy[1], lu_xy[1], 1], [ld_xy[0], ld_xy[0] * ld_xy[1], ld_xy[1], 1], @@ -114,15 +114,15 @@ class Dem: [rd_xy[0], rd_xy[0] * rd_xy[1], rd_xy[1], 1], ] ) - equa_b = np.array( + equation_b = np.array( [lu_elevation[0], ld_elevation[0], ru_elevation[0], rd_elevation[0]] ) - equa_c = np.linalg.solve(equa_a, equa_b) + equation = np.linalg.solve(equation_a, equation_b) point_z = ( - Xpixel * equa_c[0] - + equa_c[1] * Xpixel * Yline - + equa_c[2] * Yline - + equa_c[3] + Xpixel * equation[0] + + equation[1] * Xpixel * Yline + + equation[2] * Yline + + equation[3] ) site_ele[i] = point_z print(f"row:{Yline} col:{Xpixel} elev:{site_ele[i]}") @@ -132,9 +132,43 @@ class Dem: def to_line_coordination(point_x_s, point_y_s, point_x_e, point_y_e): path_length = ((point_x_s - point_x_e) ** 2 + (point_y_s - point_y_e) ** 2) ** 0.5 n = round(path_length / 1) - point_x = np.linspace(point_x_s, point_x_e, n, endpoint=True) - point_y = np.linspace(point_y_s, point_y_e, n, endpoint=True) - r = np.array([point_x, point_y]).transpose() + center_point_x = np.linspace(point_x_s, point_x_e, n, endpoint=True) + center_point_y = np.linspace(point_y_s, point_y_e, n, endpoint=True) + # 计算左右边线 + # 计算角度 + side_width = 20 # 边线宽度 + _line_angel = math.atan((point_y_e - point_y_s) / (point_x_e - point_x_s)) + if point_x_e < point_x_s: + line_angel = _line_angel + math.pi / 2 + else: + line_angel = _line_angel + left_offset_x = side_width * math.cos(line_angel + math.pi / 2) + left_offset_y = side_width * math.sin(line_angel + math.pi / 2) + left_offset_point_x = center_point_x + left_offset_x + left_offset_point_y = center_point_y + left_offset_y + right_offset_point_x = center_point_x - left_offset_x # 向左偏移和向右偏移正好是相反的 + right_offset_point_y = center_point_y - left_offset_y + # r = np.concatenate( + # ( + # np.array(left_offset_point_x), + # np.array(left_offset_point_y), + # center_point_x, + # center_point_y, + # np.array(right_offset_point_x), + # np.array(right_offset_point_y), + # ), + # axis=0, + # ) + r = np.array( + [ + np.array(left_offset_point_x), + np.array(left_offset_point_y), + center_point_x, + center_point_y, + np.array(right_offset_point_x), + np.array(right_offset_point_y), + ] + ).T return r @@ -144,23 +178,14 @@ def distance(a, b): if __name__ == "__main__": - dem = Dem(r"e:\西北院\天都山750kV线路\天都山-妙岭卫片\ASC\tianmiaoN03-dem.tif") + dem = Dem(r"d:\工程\金上线\DEM\四川-金上105-CGCS2000.tif") dem.get_dem_info(if_print=True) line_coordination = to_line_coordination( - 573284.886, 4104898.933, 574543.845, 4105825.047 + 450077.359310936, 3304781.49970352, 451068.154964846, 3304604.32630216 ) - elevation = dem.get_elevation(line_coordination) - - # fdsfd = np.array([ - # [573284.886, 4104898.933], - # [573286.584, 4104900.182], - # [573288.281, 4104901.431], - # [573289.979 , 4104902.68], - # [573291.677 , 4104903.929], - # [573293.374 , 4104905.178], - # - # ]) - # elevation = dem.get_elevation(elevation) + left_elevation = dem.get_elevation(line_coordination[:, 0:2]) + center_elevation = dem.get_elevation(line_coordination[:, 2:4]) + right_elevation = dem.get_elevation(line_coordination[:, 4:6]) doc = ezdxf.new(dxfversion="R2010") msp = doc.modelspace() x_axis = [0] @@ -168,7 +193,15 @@ if __name__ == "__main__": for cord in line_coordination[1:]: x_axis.append(distance(cord, cord_0)) msp.add_polyline2d( - [(x_axis[i] / 5, elevation[i] * 2) for i in range(len(elevation))] + [(x_axis[i] / 5, left_elevation[i] * 2) for i in range(len(left_elevation))], + dxfattribs={"color": 1}, + ) # 红色 + msp.add_polyline2d( + [(x_axis[i] / 5, center_elevation[i] * 2) for i in range(len(center_elevation))] ) + msp.add_polyline2d( + [(x_axis[i] / 5, right_elevation[i] * 2) for i in range(len(right_elevation))], + dxfattribs={"color": 5}, + ) # 蓝色 doc.saveas("d.dxf") print("Finished.")