diff --git a/dem.py b/dem.py index fe110f1..f4515d5 100644 --- a/dem.py +++ b/dem.py @@ -1,15 +1,21 @@ import math import os - +import shutil import tomli import ezdxf -from ezdxf.tools.standards import linetypes from osgeo import gdal import numpy as np import pandas as pd - +from pw import DFile, ControlFile import dem_utils +ezdxf.options.set( + "odafc-addon", + "win_exec_path", + "d:/ProgramFiles/ODAFileConverter/ODAFileConverter.exe", +) +from ezdxf.addons.odafc import export_dwg + class Dem: def __init__(self, toml_path): @@ -18,6 +24,7 @@ class Dem: self._toml_dict = toml_dict dem_file_path = toml_dict["parameter"]["dem_file_path"] self._dataset = gdal.Open(dem_file_path) + self._dem_resolution = self._dataset.GetGeoTransform()[1] def get_dem_info(self, if_print=False): """Get the information of DEM data. @@ -46,6 +53,7 @@ class Dem: def write_dxf(self): excel_pfs = self._read_path_file() + segments = [] 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] @@ -59,34 +67,34 @@ class Dem: 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") + dm_doc = ezdxf.new(dxfversion="R2004") # 设置线形 # for name, desc, pattern in linetypes(): - # if name not in doc.linetypes: - # doc.linetypes.add( + # if name not in dm_doc.linetypes: + # dm_doc.linetypes.add( # name=name, # pattern=pattern, # description=desc, # ) - msp = doc.modelspace() + dm_msp = dm_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( + dm_msp.add_polyline2d( [ (x_axis[i] / 5, left_elevation[i] * 2) for i in range(len(left_elevation)) ], dxfattribs={"color": 1}, ) # 红色 - msp.add_polyline2d( + dm_msp.add_polyline2d( [ (x_axis[i] / 5, center_elevation[i] * 2) for i in range(len(center_elevation)) ] ) - msp.add_polyline2d( + dm_msp.add_polyline2d( [ (x_axis[i] / 5, right_elevation[i] * 2) for i in range(len(right_elevation)) @@ -96,10 +104,66 @@ class Dem: 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" + # dm_doc.saveas( + # f"{out_dxf_file_dir}/D{400+int(start_point_name[1:])}.dxf" + # ) # 写断面文件 + export_dwg( + dm_doc, + f"{out_dxf_file_dir}/D{100+int(start_point_name[1:])}.dwg", + replace=True, + ) # 写断面文件 + # 写Z文件 + with open( + f"{out_dxf_file_dir}/ZD{100+int(start_point_name[1:])}", "w" + ) as z_file: + z_file.write("0 ") + z_file.write(f"{center_elevation[0]*2} ") + z_file.write("0 ") + 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") + mileage = [0] + start_point = line_coordination[0, 2:4] + for bar in line_coordination[1:, 2:4]: + mileage.append(round(dem_utils.distance(start_point, bar))) + mileage = np.array(mileage) + start_num = 4000 + int(start_point_name[1:]) + segments.append(start_num) + end_num = 4000 + int(end_point_name[1:]) + if foo == len(excel_pfs) - 2: + segments.append(end_num) + d_file = DFile( + start_num, + end_num, + mileage, + center_elevation, + left_elevation, + right_elevation, ) + d_file.save(out_dxf_file_dir) + self._copy_db_file() + tower_start_num = toml_dict["parameter"]["tower_number_start"] + control_file_template_path = toml_dict["parameter"][ + "control_file_template_path" + ] + c_file = ControlFile( + segments, tower_start_num, control_file_template_path, out_dxf_file_dir + ) + c_file.save() def get_elevation(self, site_x_y): """Get the elevation of given locations from DEM in GCS. @@ -192,14 +256,15 @@ class Dem: + equation_c[3] ) site_ele[i] = point_z - print(f"row:{Yline} col:{Xpixel} elev:{site_ele[i]}") + # 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) + dem_resolution = self._dem_resolution # dem的精度 + n = round(path_length / dem_resolution) 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) # 计算左右边线 @@ -234,3 +299,13 @@ class Dem: path_file = toml_dict["parameter"]["path_file"] excel_pds = pd.read_excel(path_file) return excel_pds + + def _copy_db_file(self): + toml_dict = self._toml_dict + db_file_dir = toml_dict["parameter"]["db_file_dir"] + out_dxf_file_dir = toml_dict["parameter"]["out_dxf_file_dir"] + shutil.copy(f"{db_file_dir}/Tower.db", f"{out_dxf_file_dir}/Tower.db") + shutil.copy(f"{db_file_dir}/ATMOS.db", f"{out_dxf_file_dir}/ATMOS.db") + shutil.copy(f"{db_file_dir}/Fit.db", f"{out_dxf_file_dir}/Fit.db") + shutil.copy(f"{db_file_dir}/RULE.db", f"{out_dxf_file_dir}/RULE.db") + shutil.copy(f"{db_file_dir}/WIRE.db", f"{out_dxf_file_dir}/WIRE.db") diff --git a/main.py b/main.py index f8daa0a..7de8980 100644 --- a/main.py +++ b/main.py @@ -1,7 +1,5 @@ - - - if __name__ == "__main__": import transformer.exporter as exporter + exporter.main() print("PyCharm") diff --git a/pw.py b/pw.py new file mode 100644 index 0000000..079c188 --- /dev/null +++ b/pw.py @@ -0,0 +1,72 @@ +import numpy as np + + +# D文件 +class DFile: + def __init__( + self, start_num, end_num, mileage: np.ndarray, center_z, left_z, right_z + ): + self.start_num = start_num + self.end_num = end_num + self._mileage = mileage + self._center_z = center_z + self._left_z = left_z + self._right_z = right_z + pass + + def save(self, out_dir): + start_num = self.start_num + end_num = self.end_num + mileage = self._mileage + center_z = self._center_z + left_z = self._left_z + right_z = self._right_z + mat = np.zeros((mileage.size, 5)) + mat[0, 0] = start_num + mat[-1, 0] = end_num + mat[:, 1] = mileage + mat[:, 2] = center_z + mat[:, 3] = left_z + mat[:, 4] = right_z + with open(f"{out_dir}/D{100+start_num-4000}.txt", "w") as d_file: + for foo in mat: + d_file.write(f"{foo[0]:.0f}\t") + d_file.write(f"{foo[1]:.0f}\t") + d_file.write(f"{foo[2]}\t") + d_file.write(f"{foo[3]}\t") + d_file.write(f"{foo[4]}\t") + d_file.write("\n") + with open(f"{out_dir}/X{100+start_num-4000}.txt", "w") as x_file: + x_file.write("dm\n") + x_file.write(f"{mileage.size}\n") + x_file.write("0\n") + x_file.write(f"{mileage.size*3}\n") + x_file.write("1\n") + x_file.write("0\n") + x_file.write("0\n") + pass + + +class ControlFile: + def __init__(self, segments, tower_start_number, c_file_template_path, out_dir): + self._c_file_template_path = c_file_template_path + self._out_dir = out_dir + self._segments = segments + self._tower_start_number = tower_start_number + pass + + def save(self): + c_file_template_path = self._c_file_template_path + out_dir = self._out_dir + segments = self._segments + tower_start_number = self._tower_start_number + with open(c_file_template_path) as c_file_template: + c_template = c_file_template.read() + c_template=c_template.replace("{SEGMENT_N}", str(len(segments))) + c_template=c_template.replace("{TOWER_NUM_START}", f'{tower_start_number}') + c_template=c_template.replace("{SEGMENT}", "\n".join([str(seg-4000+100) for seg in segments])) + c_template=c_template.replace("{SEGMENT_S}", f"{segments[0]-4000+100}") + c_template=c_template.replace("{SEGMENT_E}", f"{segments[-1]-4000+100}") + c_template=c_template.replace("{CONTROL_FILE}", f"S{segments[0]-4000+100}") + with open(f"{out_dir}/C{segments[0]-4000+100}.dat", "w") as c_file: + c_file.write(c_template) diff --git a/transformer/__init__.py b/transformer/__init__.py index ab50b67..6887633 100644 --- a/transformer/__init__.py +++ b/transformer/__init__.py @@ -12,8 +12,9 @@ class Projector: EPSG_WGS_84 = "epsg:4326" -EPSG_Xian_1980_105E = "epsg:4507" +EPSG_CGCS_2000_105E = "epsg:4507" EPSG_CGCS_2000_108E = "epsg:4545" +EPSG_Xian_1980_114E = "epsg:2383" # _transformer_84_yx = Transformer.from_crs("epsg:4326", "epsg:4507") # _transformer_yx_84 = Transformer.from_crs("epsg:4507", "epsg:4326") diff --git a/transformer/exporter.py b/transformer/exporter.py index 542261d..59be894 100644 --- a/transformer/exporter.py +++ b/transformer/exporter.py @@ -16,17 +16,15 @@ def prettify(elem): def main(): - df = pd.read_excel("成果表-0120.xlsx") - needed_df = df.iloc[:120, :3] - trans = transformer.Projector( - transformer.EPSG_CGCS_2000_108E, transformer.EPSG_WGS_84 - ) - y_array = needed_df["北坐标(m)"].to_numpy() - x_array = needed_df["东坐标(m)"].to_numpy() + df = pd.read_excel("杆塔成果表(左侧）-复测.xls") + needed_df = df.iloc[:85, :] + trans = transformer.Projector(transformer.EPSG_Xian_1980_114E, transformer.EPSG_WGS_84) + y_array = needed_df["北坐标"].to_numpy() + x_array = needed_df["东坐标"].to_numpy() yx_array = np.array([y_array, x_array]).transpose() wgs_84_point = np.array([trans.transform(y, x) for y, x in yx_array]) name_lat_lon_df = pd.DataFrame( - {"桩号": needed_df["点号"], "纬度": wgs_84_point[:, 0], "经度": wgs_84_point[:, 1]} + {"桩号": needed_df["杆塔"], "纬度": wgs_84_point[:, 0], "经度": wgs_84_point[:, 1]} ) ET.register_namespace("", "http://www.opengis.net/kml/2.2") ns = { @@ -61,5 +59,5 @@ def main(): "{lon},{lat},{elevation}".format(lon=foo[1], lat=foo[0], elevation=foo[2]) for foo in line_coordination_list ) - with open("out.kml", "w", encoding="utf-8") as f: + with open("周口.kml", "w", encoding="utf-8") as f: f.write(prettify(xml_root).decode("utf-8")) diff --git a/transmission_line/plane.py b/transmission_line/plane.py index b42b219..994e8d1 100644 --- a/transmission_line/plane.py +++ b/transmission_line/plane.py @@ -4,26 +4,101 @@ import ezdxf from data_types import WGS84Point import transformer import numpy as np +from functools import singledispatch +from typing import Tuple -def generate_wgs_84_point(point_start: WGS84Point, point_end: WGS84Point, n_point): - _transformer_84_yx = transformer.Projector("epsg:4326", "epsg:4507") - start_y, start_x = _transformer_84_yx.transform( - point_start.latitude, point_start.longitude - ) - end_y, end_x = _transformer_84_yx.transform(point_end.latitude, point_end.longitude) - # 插入点 - x_linspace = np.linspace(start_x, end_x, n_point) - y_linspace = np.linspace(start_y, end_y, n_point) - elevation_linspace = np.linspace( - point_start.elevation, point_end.elevation, n_point - ) - yxz_point = np.array([y_linspace, x_linspace, elevation_linspace]).transpose() - _transformer_yx_84 = transformer.Projector("epsg:4507", "epsg:4326") - wgs_84_point = np.array( - [np.array(_transformer_yx_84.transform(y, x)) for y, x, _elevation in yxz_point] - ) - return wgs_84_point +class Plane: + # dem文件坐标系统 + def __init__(self, dem_system, plane_system, dem_path): + self._dem_path = dem_path + self._dem_system = dem_system + self._start_84_point = None + self._end_84_point = None + self._plane_system = plane_system + self._transformer_84_to_dem_ = transformer.Projector( + transformer.EPSG_WGS_84, dem_system + ) + self._transformer_dem_to_84 = transformer.Projector( + dem_system, transformer.EPSG_WGS_84 + ) + + @singledispatch + def set_line(self, start: WGS84Point, end: WGS84Point): + self._start_84_point = start + self._end_84_point = end + pass + + # def set_line_xy(self, start: Tuple[float], end: Tuple[float]): + @set_line.register(tuple) + def set_line_yx(self, start: Tuple[float, float], end: Tuple[float, float]): + _transformer = self._transformer_dem_to_84 + start_wgs_84 = self._transformer_dem_to_84.transform(*start) + end_wgs_84 = self._transformer_dem_to_84.transform(*end) + self._start_84_point = WGS84Point(start_wgs_84[0], start_wgs_84[1], 0) + self._end_84_point = WGS84Point(end_wgs_84[0], end_wgs_84[1], 0) + pass + + def get_elevation(self, n): + point_start = self._start_84_point + point_end = self._end_84_point + all_wgs_84_point = self._generate_wgs_84_point(point_start, point_end, n) + # 转换为DEM的坐标 + transformer_dem_to_84 = self._transformer_84_to_dem_ + all_dem_point = np.array( + [ + transformer_dem_to_84.transform(y_lat, x_lon) + for y_lat, x_lon in all_wgs_84_point[:, [0, 1]] + ] + ) + dem_gcs = dem.Dem(self._dem_path) + elevation_in_line = dem_gcs.get_elevation(all_dem_point) + return np.array(elevation_in_line) + + # 计算档距 + + def cal_span(self): + start = self._start_84_point + end = self._end_84_point + _transformer_84_yx = transformer.Projector( + transformer.EPSG_WGS_84, transformer.EPSG_Xian_1980_114E + ) + start_point = np.array( + _transformer_84_yx.transform(start.latitude, start.longitude) + ) + end_point = np.array(_transformer_84_yx.transform(end.latitude, end.longitude)) + return int(np.round(np.sum((start_point - end_point) ** 2) ** 0.5)) + + def _generate_wgs_84_point( + self, point_start: WGS84Point, point_end: WGS84Point, n_point + ): + plane_system = self._plane_system + _transformer_84_yx = transformer.Projector( + transformer.EPSG_WGS_84, plane_system + ) + start_y, start_x = _transformer_84_yx.transform( + point_start.latitude, point_start.longitude + ) + end_y, end_x = _transformer_84_yx.transform( + point_end.latitude, point_end.longitude + ) + # 插入点 + x_linspace = np.linspace(start_x, end_x, n_point) + y_linspace = np.linspace(start_y, end_y, n_point) + elevation_linspace = np.linspace( + point_start.elevation, point_end.elevation, n_point + ) + yxz_point = np.array([y_linspace, x_linspace, elevation_linspace]).transpose() + _transformer_yx_84 = transformer.Projector( + plane_system, transformer.EPSG_WGS_84 + ) + wgs_84_point = np.array( + [ + np.array(_transformer_yx_84.transform(y, x)) + for y, x, _elevation in yxz_point + ] + ) + return wgs_84_point def draw_dxf(elevation_in_line, interval, dxf_path): @@ -40,28 +115,26 @@ def draw_dxf(elevation_in_line, interval, dxf_path): / 5 ) # 平滑用 func = interpolate.interp1d(points[0, :], points[1, :], kind="cubic") - msp.add_polyline2d([(x + 10, func(x)) for x in new_x_axis]) + # msp.add_polyline2d([(x + 10, func(x)) for x in new_x_axis]) msp.add_polyline2d([(x + 10, y) for x, y in points.transpose()]) doc.saveas(dxf_path) -# 计算档距 -def cal_span(start: WGS84Point, end: WGS84Point): - _transformer_84_yx = transformer.Projector("epsg:4326", "epsg:4507") - start_point = np.array( - _transformer_84_yx.transform(start.latitude, start.longitude) - ) - end_point = np.array(_transformer_84_yx.transform(end.latitude, end.longitude)) - return np.round(np.sum((start_point - end_point) ** 2) ** 0.5) - - def main(): - dem_gcs = dem.Dem(r"d:\宁夏自治区.tif") - dem_gcs.get_dem_info(if_print=False) - start_point = WGS84Point(36.74138039, 105.58954375, 1670) - end_point = WGS84Point(36.72869989, 105.61023855, 1718) - span = int(cal_span(start_point, end_point)) + _plane = Plane( + transformer.EPSG_Xian_1980_114E, transformer.EPSG_Xian_1980_114E, r"d:\307.tif" + ) + # start_point = WGS84Point(33.24643289, 114.58574497, 28) + # end_point = WGS84Point(33.22873308, 114.57951634, 28) + _plane.set_line_yx((3680231.27, 554470.91), (3678264.93, 553901.32)) + span = _plane.cal_span() print("档距 {span}".format(span=span)) - wgs_84_point_in_line = generate_wgs_84_point(start_point, end_point, span) - elevation_in_line = dem_gcs.get_elevation(wgs_84_point_in_line) - draw_dxf(elevation_in_line, 1, "a.dxf") + elevation_in_line = _plane.get_elevation(int(span / 5)) + # dem_gcs = dem.Dem() + # dem_gcs.get_dem_info(if_print=False) + # + # span = int(cal_span(start_point, end_point)) + # print("档距 {span}".format(span=span)) + # wgs_84_point_in_line = generate_wgs_84_point(start_point, end_point, span) + # elevation_in_line = dem_gcs.get_elevation(wgs_84_point_in_line) + draw_dxf(elevation_in_line, 5, "a.dxf")