diff --git a/dem.py b/dem.py index 68a7035..fe110f1 100644 --- a/dem.py +++ b/dem.py @@ -1,13 +1,23 @@ import math +import os +import tomli import ezdxf +from ezdxf.tools.standards import linetypes from osgeo import gdal import numpy as np +import pandas as pd + +import dem_utils class Dem: - def __init__(self, filepath): - self._dataset = gdal.Open(filepath) + def __init__(self, toml_path): + with open(toml_path, "rb") as tf: + toml_dict = tomli.load(tf) + self._toml_dict = toml_dict + dem_file_path = toml_dict["parameter"]["dem_file_path"] + self._dataset = gdal.Open(dem_file_path) def get_dem_info(self, if_print=False): """Get the information of DEM data. @@ -34,7 +44,64 @@ class Dem: return dem_row, dem_col, dem_band, dem_gt, dem_proj - def get_elevation(self, site_latlng): + def write_dxf(self): + excel_pfs = self._read_path_file() + 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] + point_x_s = float(excel_pfs.iloc[foo, 1]) + point_y_s = float(excel_pfs.iloc[foo, 2]) + point_x_e = float(excel_pfs.iloc[foo + 1, 1]) + point_y_e = float(excel_pfs.iloc[foo + 1, 2]) + line_coordination = self.to_line_coordination( + point_x_s, point_y_s, point_x_e, point_y_e + ) + left_elevation = self.get_elevation(line_coordination[:, 0:2]) + center_elevation = self.get_elevation(line_coordination[:, 2:4]) + right_elevation = self.get_elevation(line_coordination[:, 4:6]) + doc = ezdxf.new(dxfversion="R2010") + # 设置线形 + # for name, desc, pattern in linetypes(): + # if name not in doc.linetypes: + # doc.linetypes.add( + # name=name, + # pattern=pattern, + # description=desc, + # ) + msp = doc.modelspace() + x_axis = [0] + cord_0 = line_coordination[0, 2:4] # 取中线的 + for cord in line_coordination[1:, 2:4]: + x_axis.append(dem_utils.distance(cord, cord_0)) + msp.add_polyline2d( + [ + (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}, + ) # 蓝色 + 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) + dem_prefix = toml_dict["parameter"]["dem_prefix"] + doc.saveas( + f"{out_dxf_file_dir}/{dem_prefix}{start_point_name}_{end_point_name}.dxf" + ) + + def get_elevation(self, site_x_y): """Get the elevation of given locations from DEM in GCS. Parameters: dem_gcs -- The input DEM (in GCS). @@ -52,9 +119,9 @@ class Dem: gt = gdal_data.GetGeoTransform() # print("\nThe 6 GeoTransform parameters of DEM are:\n", gt) - N_site = site_latlng.shape[0] - Xgeo = site_latlng[:, 0] - Ygeo = site_latlng[:, 1] + N_site = site_x_y.shape[0] + Xgeo = site_x_y[:, 0] + Ygeo = site_x_y[:, 1] site_ele = np.zeros(N_site) for i in range(N_site): # Note: @@ -105,7 +172,7 @@ class Dem: # * delta_y # ) - # 通过空间平面方程乃拟合Z值 参考《数字高程模型教程》 汤国安，李发源，刘学军编著 p89 + # 通过空间平面方程拟合Z值 参考《数字高程模型教程》 汤国安，李发源，刘学军编著 p89 equation_a = np.array( [ [lu_xy[0], lu_xy[0] * lu_xy[1], lu_xy[1], 1], @@ -117,91 +184,53 @@ class Dem: equation_b = np.array( [lu_elevation[0], ld_elevation[0], ru_elevation[0], rd_elevation[0]] ) - equation = np.linalg.solve(equation_a, equation_b) + equation_c = np.linalg.solve(equation_a, equation_b) point_z = ( - Xpixel * equation[0] - + equation[1] * Xpixel * Yline - + equation[2] * Yline - + equation[3] + Xpixel * equation_c[0] + + equation_c[1] * Xpixel * Yline + + equation_c[2] * Yline + + equation_c[3] ) site_ele[i] = point_z print(f"row:{Yline} col:{Xpixel} elev:{site_ele[i]}") return site_ele + def to_line_coordination(self, 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) + 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) + # 计算左右边线 + # 计算角度 + toml_dict = self._toml_dict + side_width = toml_dict["parameter"]["side_width"] # 边线宽度 + _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.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 -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) - 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 - - -def distance(a, b): - r = np.sum((a - b) ** 2) ** 0.5 - return r - - -if __name__ == "__main__": - dem = Dem(r"d:\工程\金上线\DEM\四川-金上105-CGCS2000.tif") - dem.get_dem_info(if_print=True) - line_coordination = to_line_coordination( - 450077.359310936, 3304781.49970352, 451068.154964846, 3304604.32630216 - ) - 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] - cord_0 = line_coordination[0, :] - for cord in line_coordination[1:]: - x_axis.append(distance(cord, cord_0)) - msp.add_polyline2d( - [(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.") + def _read_path_file(self): + toml_dict = self._toml_dict + path_file = toml_dict["parameter"]["path_file"] + excel_pds = pd.read_excel(path_file) + return excel_pds diff --git a/dem_utils.py b/dem_utils.py new file mode 100644 index 0000000..0e41502 --- /dev/null +++ b/dem_utils.py @@ -0,0 +1,6 @@ +import numpy as np + + +def distance(a, b): + r = np.sum((a - b) ** 2) ** 0.5 + return r diff --git a/main_dem.py b/main_dem.py new file mode 100644 index 0000000..13b552b --- /dev/null +++ b/main_dem.py @@ -0,0 +1,8 @@ +from dem_utils import distance +from dem import Dem +import ezdxf +if __name__ == "__main__": + dem = Dem('db.toml') + dem.get_dem_info(if_print=True) + dem.write_dxf() + print("Finished.") \ No newline at end of file