diff --git a/dem.py b/dem.py index 2620856..76be2c2 100644 --- a/dem.py +++ b/dem.py @@ -48,6 +48,10 @@ class Dem: def write_dxf(self): excel_pfs = self._read_path_file() segments = [] + plate_doc = ezdxf.new(dxfversion="R2010") + plate_msp = plate_doc.modelspace() + toml_dict = self._toml_dict + out_dxf_file_dir = toml_dict["parameter"]["out_dxf_file_dir"] for foo in range(len(excel_pfs) - 1): start_point_name: str = excel_pfs.iloc[foo, 0] end_point_name: str = excel_pfs.iloc[foo + 1, 0] @@ -96,6 +100,7 @@ class Dem: dxfattribs={"color": 5}, ) # 蓝色 # 树的线 + # TODO 没有考虑用最高边线的情况 if self._tree_height > 0: dm_msp.add_polyline2d( [ @@ -103,12 +108,7 @@ class Dem: for i in range(len(center_elevation)) ] ) - toml_dict = self._toml_dict - out_dxf_file_dir = toml_dict["parameter"]["out_dxf_file_dir"] os.makedirs(out_dxf_file_dir, exist_ok=True) - # dm_doc.saveas( - # f"{out_dxf_file_dir}/D{400+int(start_point_name[1:])}.dxf" - # ) # 写断面文件 ezdxf.options.set( "odafc-addon", "win_exec_path", @@ -131,20 +131,17 @@ class Dem: z_file.write("0 ") z_file.write(f"{center_elevation[0]*2-50}") # TODO:写平面文件 - # plate_doc = ezdxf.new(dxfversion="R2010") - # plate_msp = plate_doc.modelspace() - # plate_msp.add_polyline2d( - # line_coordination[:, 0:2], - # dxfattribs={"color": 1}, - # ) # 红色 - # plate_msp.add_polyline2d( - # line_coordination[:, 2:4], - # ) - # plate_msp.add_polyline2d( - # line_coordination[:, 4:6], - # dxfattribs={"color": 5}, - # ) # 蓝色 - # plate_doc.saveas(f"{out_dxf_file_dir}/plate.dxf") + plate_msp.add_polyline2d( + line_coordination[:, 0:2], + dxfattribs={"color": 1}, + ) # 红色 + plate_msp.add_polyline2d( + line_coordination[:, 2:4], + ) + plate_msp.add_polyline2d( + line_coordination[:, 4:6], + dxfattribs={"color": 5}, + ) # 蓝色 mileage = [0] start_point = line_coordination[0, 2:4] for bar in line_coordination[1:, 2:4]: @@ -162,7 +159,7 @@ class Dem: center_elevation, left_elevation, right_elevation, - self._tree_height + self._tree_height, ) d_file.save(out_dxf_file_dir) self._copy_db_file() @@ -174,6 +171,7 @@ class Dem: segments, tower_start_num, control_file_template_path, out_dxf_file_dir ) c_file.save() + plate_doc.saveas(f"{out_dxf_file_dir}/plate.dxf") def get_elevation(self, site_x_y): """Get the elevation of given locations from DEM in GCS. @@ -219,10 +217,14 @@ class Dem: Xpixel = (Xgeo[i] - gt[0]) / gt[1] Yline = (Ygeo[i] - gt[3]) / gt[5] # 寻找左上，左下，右上，右下4个点 - lu_xy = np.array([math.floor(Xpixel), math.floor(Yline)]) # 左上 - ld_xy = np.array([math.floor(Xpixel), math.ceil(Yline)]) # 左下 - ru_xy = np.array([math.ceil(Xpixel), math.floor(Yline)]) # 右上 - rd_xy = np.array([math.ceil(Xpixel), math.ceil(Yline)]) # 右下 + # lu_xy = np.array([math.floor(Xpixel), math.floor(Yline)]) # 左上 + # ld_xy = np.array([math.floor(Xpixel), math.ceil(Yline)]) # 左下 + # ru_xy = np.array([math.ceil(Xpixel), math.floor(Yline)]) # 右上 + # rd_xy = np.array([math.ceil(Xpixel), math.ceil(Yline)]) # 右下 + lu_xy = np.array([math.floor(Xpixel), math.ceil(Yline)]) # 左上 + ld_xy = np.array([math.floor(Xpixel), math.floor(Yline)]) # 左下 + ru_xy = np.array([math.ceil(Xpixel), math.ceil(Yline)]) # 右上 + rd_xy = np.array([math.ceil(Xpixel), math.floor(Yline)]) # 右下 lu_elevation = gdal_data.ReadAsArray( int(lu_xy[0]), int(lu_xy[1]), 1, 1 ).astype(float) @@ -247,26 +249,49 @@ class Dem: # ) # 通过空间平面方程拟合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], - [ru_xy[0], ru_xy[0] * ru_xy[1], ru_xy[1], 1], - [rd_xy[0], rd_xy[0] * rd_xy[1], rd_xy[1], 1], - ] - ) - equation_b = np.array( - [lu_elevation[0], ld_elevation[0], ru_elevation[0], rd_elevation[0]] - ) - equation_c = np.linalg.solve(equation_a, equation_b) + # 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], + # [ru_xy[0], ru_xy[0] * ru_xy[1], ru_xy[1], 1], + # [rd_xy[0], rd_xy[0] * rd_xy[1], rd_xy[1], 1], + # ] + # ) + # equation_b = np.array( + # [lu_elevation[0], ld_elevation[0], ru_elevation[0], rd_elevation[0]] + # ) + # equation_c = np.linalg.solve(equation_a, equation_b) + # point_z = ( + # Xpixel * equation_c[0] + # + equation_c[1] * Xpixel * Yline + # + equation_c[2] * Yline + # + equation_c[3] + # ) point_z = ( - Xpixel * equation_c[0] - + equation_c[1] * Xpixel * Yline - + equation_c[2] * Yline - + equation_c[3] + lu_elevation[0] + * ( + (1 - math.fabs(Xpixel - lu_xy[0])) + * (1 - math.fabs(Yline - lu_xy[1])) + ) + + ld_elevation[0] + * ( + (1 - math.fabs(Xpixel - ld_xy[0])) + * (1 - math.fabs(Yline - ld_xy[1])) + ) + + ru_elevation[0] + * ( + (1 - math.fabs(Xpixel - ru_xy[0])) + * (1 - math.fabs(Yline - ru_xy[1])) + ) + + rd_elevation[0] + * ( + (1 - math.fabs(Xpixel - rd_xy[0])) + * (1 - math.fabs(Yline - rd_xy[1])) + ) ) site_ele[i] = point_z # print(f"row:{Yline} col:{Xpixel} elev:{site_ele[i]}") + pass return site_ele def to_line_coordination(self, point_x_s, point_y_s, point_x_e, point_y_e): diff --git a/main.py b/main.py index 7de8980..b3c035e 100644 --- a/main.py +++ b/main.py @@ -1,5 +1,13 @@ +import sys +from dem import Dem +from loguru import logger if __name__ == "__main__": - import transformer.exporter as exporter - - exporter.main() - print("PyCharm") + if len(sys.argv) < 2: + toml_file_path = r"db_JS.toml" + else: + toml_file_path = sys.argv[1] + logger.info(f'读取配置文件{toml_file_path}') + dem = Dem(toml_file_path) + dem.get_dem_info(if_print=True) + dem.write_dxf() + print("Finished.")