.gitignore

@@ -1,16 +0,0 @@
-build
-dist
-out_dxf
-__pycache__
-.idea
-*.spec
-canva
-*.xls*
-*.toml
-control_file_template.txt
-.idea
-.vscode
-__pycache__
-Lib
-Scripts
-example

Makefile

@@ -1,2 +0,0 @@
-all:
-	pyinstaller -F main_dem.py -n Dem_NWEPDI

cad.py

@@ -1,50 +0,0 @@
-from typing import List
-import os
-import subprocess
-
-
-def convert_dxf_to_dwg(dxf_files: List[str]):
-    cad_file_path = r"D:/Program Files/Cad2022/AutoCAD 2022/accoreconsole.exe"
-    first_dwg_file = dxf_files[0]
-    dir_path = os.path.dirname(first_dwg_file)
-    script_path = f"{dir_path}/batch.scr"
-    for _, dxf in enumerate(dxf_files):
-        with open(script_path, "w", encoding="ansi") as f:
-            file_name_path, _ = os.path.splitext(dxf)
-            new_dwg_name = f"{file_name_path}.dwg"
-            if os.path.exists(new_dwg_name):
-                os.remove(new_dwg_name)
-            new_dwg_name = new_dwg_name.replace("//", "/")
-            f.write(f'saveas 2004 "{new_dwg_name}"\n')
-            cmd = rf'"{cad_file_path}" /s "{script_path}" /i "{dxf}" /iso'
-        cmd = cmd.replace("//", "/")
-        subprocess.call(
-            cmd,
-            stderr=subprocess.DEVNULL,
-            stdin=subprocess.DEVNULL,
-            stdout=subprocess.DEVNULL,
-        )
-
-
-def convert_dxf_to_dwg1(dxf_files: List[str]):
-    cad_file_path = r"D:/Program Files/Cad2022/AutoCAD 2022/acad.exe"
-    first_dwg_file = dxf_files[0]
-    dir_path = os.path.dirname(first_dwg_file)
-    script_path = f"{dir_path}/batch.scr"
-    with open(script_path, "w", encoding="ansi") as f:
-        f.write("FILEDIA 0\n")
-        for ind, dxf in enumerate(dxf_files):
-            f.write(f'_open "{dxf}"\n')
-            file_name_path, _ = os.path.splitext(dxf)
-            new_dwg_name = f"{file_name_path}.dwg"
-            if os.path.exists(new_dwg_name):
-                os.remove(new_dwg_name)
-            f.write(f'saveas 2004 "{new_dwg_name}"\n')
-            if ind == len(dxf_files) - 1:
-                f.write("FILEDIA 1\n")  # 最后一个回复打开dialog
-                f.write("qquit\n")
-            else:
-                f.write(f"_close\n")
-        # f.write("qquit\n")
-    cmd = rf'"{cad_file_path}" /b "{script_path}"'
-    subprocess.Popen(cmd)

dem.py

@@ -1,26 +1,13 @@
 import math
-import os
+
-import shutil
-import tomli
 import ezdxf
 from osgeo import gdal
 import numpy as np
-import pandas as pd
-from pw import DFile, ControlFile
-import dem_utils
-from nwed import Nwed
-import cad
 
 
 class Dem:
-    def __init__(self, toml_path):
+    def __init__(self, filepath):
-        with open(toml_path, "rb") as tf:
+        self._dataset = gdal.Open(filepath)
-            toml_dict = tomli.load(tf)
-        self._toml_dict = toml_dict
-        dem_file_path = toml_dict["parameter"]["dem_file_path"]
-        self._tree_height = toml_dict["parameter"]["tree_height"]
-        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.
@@ -47,249 +34,7 @@ class Dem:
 
         return dem_row, dem_col, dem_band, dem_gt, dem_proj
 
-    def write(self):
+    def get_elevation(self, site_latlng):
-        # TODO:不应该设置缩放因数
-        dxfs = []
-        zoom_factor = 1
-        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"]
-        # 写整个断面
-        dm_whole_doc = ezdxf.new(dxfversion="R2010")
-        dm_whole_accumulative_distance = 0  # 累加里程
-        dm_whole_msp = dm_whole_doc.modelspace()
-        for foo in range(len(excel_pfs) - 1):
-            start_point_name: str = excel_pfs.iloc[foo, 3]
-            end_point_name: str = excel_pfs.iloc[foo + 1, 3]
-            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])
-            dm_doc = ezdxf.new(dxfversion="R2010")
-            # 设置线形
-            # for name, desc, pattern in linetypes():
-            #     if name not in dm_doc.linetypes:
-            #         dm_doc.linetypes.add(
-            #             name=name,
-            #             pattern=pattern,
-            #             description=desc,
-            #         )
-            dm_msp = dm_doc.modelspace()
-            x_axis = [0]
-            cord_0 = line_coordination[0, 2:4]  # 取中线的x轴作为横坐标
-            for cord in line_coordination[1:, 2:4]:
-                x_axis.append(dem_utils.distance(cord, cord_0))
-            # 左边线
-            left_line = [
-                (x_axis[i] / 5 / zoom_factor, left_elevation[i] * 2 / zoom_factor)
-                for i in range(len(left_elevation))
-            ]
-            # dm_whole_msp.add_polyline2d([[0,0],[100,100]])
-            # dm_doc.saveas('f.dxf')
-            dm_whole_msp.add_polyline2d(
-                np.array(left_line)
-                + np.hstack(
-                    (
-                        dm_whole_accumulative_distance * np.ones((len(x_axis), 1)) / 5,
-                        np.zeros((len(x_axis), 1)),
-                    )
-                ),
-                dxfattribs={"color": 1},
-            )  # 红色
-            dm_msp.add_polyline2d(left_line, dxfattribs={"color": 1})  # 红色
-            mid_line = [
-                (x_axis[i] / 5 / zoom_factor, center_elevation[i] * 2 / zoom_factor)
-                for i in range(len(center_elevation))
-            ]  # 中线
-            dm_whole_msp.add_polyline2d(
-                np.array(mid_line)
-                + np.hstack(
-                    (
-                        dm_whole_accumulative_distance * np.ones((len(x_axis), 1)) / 5,
-                        np.zeros((len(x_axis), 1)),
-                    )
-                )
-            )
-            dm_msp.add_polyline2d(mid_line)
-            # 右边线
-            right_line = [
-                (x_axis[i] / 5 / zoom_factor, right_elevation[i] * 2 / zoom_factor)
-                for i in range(len(right_elevation))
-            ]
-            dm_whole_msp.add_polyline2d(
-                np.array(right_line)
-                + np.hstack(
-                    (
-                        dm_whole_accumulative_distance * np.ones((len(x_axis), 1)) / 5,
-                        np.zeros((len(x_axis), 1)),
-                    )
-                ),
-                dxfattribs={"color": 5},  # 蓝色
-            )
-            dm_msp.add_polyline2d(
-                right_line,
-                dxfattribs={"color": 5},
-            )  # 蓝色
-            # 树的线
-            # 考虑用最高边线的情况
-            tree_line = [
-                (
-                    x_axis[i] / 5 / zoom_factor,
-                    (
-                        np.max(
-                            (
-                                center_elevation[i],
-                                left_elevation[i],
-                                right_elevation[i],
-                            )
-                        )
-                        + self._tree_height
-                    )
-                    * 2
-                    / zoom_factor,
-                )
-                for i in range(len(center_elevation))
-            ]
-            if self._tree_height > 0:
-                dm_whole_msp.add_polyline2d(
-                    np.array(tree_line)
-                    + np.hstack(
-                        (
-                            dm_whole_accumulative_distance
-                            * np.ones((len(x_axis), 1))
-                            / 5,
-                            np.zeros((len(x_axis), 1)),
-                        )
-                    ),
-                    dxfattribs={"color": 5},
-                )
-                dm_msp.add_polyline2d(tree_line, dxfattribs={"color": 5})
-            dm_whole_accumulative_distance += x_axis[-1]
-            os.makedirs(out_dxf_file_dir, exist_ok=True)
-            # ezdxf.options.set(
-            #     "odafc-addon",
-            #     "win_exec_path",
-            #     "d:/ProgramFiles/ODAFileConverter/ODAFileConverter.exe",
-            # )
-            # TODO 去掉odafc依赖
-            # from ezdxf.addons.odafc import export_dwg
-
-            # export_dwg(
-            #     dm_doc,
-            #     f"{out_dxf_file_dir}/D{100+int(start_point_name[1:])}.dwg",
-            #     replace=True,
-            # )  # 写断面文件
-            dm_doc_dxf = f"{out_dxf_file_dir}/D{100+int(start_point_name[1:])}.dxf"
-            dm_doc.saveas(dm_doc_dxf)
-            dxfs.append(dm_doc_dxf)
-            # 写Z文件
-            z_file_path = f"{out_dxf_file_dir}/ZD{100+int(start_point_name[1:])}"
-            with open(z_file_path, "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}")
-            if foo == 0:
-                # copy file a to dist b
-                shutil.copy(z_file_path, f"{out_dxf_file_dir}/ZDA")
-            # 写平面文件
-            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]:
-                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,
-                self._tree_height,
-            )
-            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()
-            ##################
-            ##################
-            # 写nwed
-            section_name = f"{100+int(start_point_name[1:])}"
-            licheng = f"{mileage[-1]}"
-            remarks_start = f"{int(start_point_name[1:])}"
-            pole_name_start = f"{start_num}"
-            remarks_end = f"{int(end_point_name[1:])}"
-            pole_name_end = f"{end_num}"
-            toml_dict = self._toml_dict
-            side_width = toml_dict["parameter"]["side_width"]  # 边线宽度
-            verticalExtent2 = f"{side_width}"
-            mid_vec_list = list(map(lambda x: [f"{x[0]*5}", f"{x[1]/2}"], mid_line))
-            right_vec_list = list(map(lambda x: [f"{x[0]*5}", f"{x[1]/2}"], right_line))
-            left_vec_list = list(map(lambda x: [f"{x[0]*5}", f"{x[1]/2}"], left_line))
-            nwed = Nwed(
-                section_name,
-                licheng,
-                remarks_start,
-                pole_name_start,
-                remarks_end,
-                pole_name_end,
-                verticalExtent2,
-                mid_vec_list,
-                right_vec_list,
-                left_vec_list,
-            )
-            nwed.write(
-                f"{out_dxf_file_dir}/{100+int(start_point_name[1:])}.nwed",
-            )
-        # 写整个断面文件
-        # export_dwg(
-        #     dm_whole_doc,
-        #     f"{out_dxf_file_dir}/DA.dwg",
-        #     replace=True,
-        # )
-        dm_whole_doc_dxf = f"{out_dxf_file_dir}/DA.dxf"
-        dm_whole_doc.saveas(dm_whole_doc_dxf)
-        # cad.convert_dxf_to_dwg([f"{out_dxf_file_dir}/DA.dxf"])
-        dxfs.append(dm_whole_doc_dxf)
-        # 写平面文件
-        plate_doc_dxf = f"{out_dxf_file_dir}/plate.dxf"
-        plate_doc.saveas(plate_doc_dxf)
-        dxfs.append(plate_doc_dxf)
-        cad.convert_dxf_to_dwg(dxf_files=dxfs)
-
-    def get_elevation(self, site_x_y):
         """Get the elevation of given locations from DEM in GCS.
         Parameters:
             dem_gcs <osgeo.gdal.Dataset> -- The input DEM (in GCS).
@@ -307,9 +52,9 @@ class Dem:
         gt = gdal_data.GetGeoTransform()
         # print("\nThe 6 GeoTransform parameters of DEM are:\n", gt)
 
-        N_site = site_x_y.shape[0]
+        N_site = site_latlng.shape[0]
-        Xgeo = site_x_y[:, 0]
+        Xgeo = site_latlng[:, 0]
-        Ygeo = site_x_y[:, 1]
+        Ygeo = site_latlng[:, 1]
         site_ele = np.zeros(N_site)
         for i in range(N_site):
             # Note:
@@ -333,10 +78,10 @@ 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.ceil(Yline)])  # 左上
+            lu_xy = np.array([math.floor(Xpixel), math.floor(Yline)])  # 左上
-            ld_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.ceil(Yline)])  # 右上
+            ru_xy = np.array([math.ceil(Xpixel), math.floor(Yline)])  # 右上
-            rd_xy = np.array([math.ceil(Xpixel), math.floor(Yline)])  # 右下
+            rd_xy = np.array([math.ceil(Xpixel), math.ceil(Yline)])  # 右下
             lu_elevation = gdal_data.ReadAsArray(
                 int(lu_xy[0]), int(lu_xy[1]), 1, 1
             ).astype(float)
@@ -349,85 +94,114 @@ class Dem:
             rd_elevation = gdal_data.ReadAsArray(
                 int(rd_xy[0]), int(rd_xy[1]), 1, 1
             ).astype(float)
-            # 依据https://blog.csdn.net/jameschen9051/article/details/109469228中的公司
+            # delta_x = (Xpixel - ld_xy[0]) / 1  # 距离是1
+            # delta_y = (ld_xy[1]-Yline) / 1
+            # site_ele[i] = (
+            #     ld_elevation
+            #     + (lu_elevation - ld_elevation) * delta_y
+            #     + (rd_elevation - ld_elevation) * delta_x
+            #     + (ld_elevation - lu_elevation + ru_elevation - rd_elevation)
+            #     * delta_x
+            #     * delta_y
+            # )
+
+            # 通过空间平面方程乃拟合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 = np.linalg.solve(equation_a, equation_b)
             point_z = (
-                lu_elevation[0]
+                Xpixel * equation[0]
-                * (
+                + equation[1] * Xpixel * Yline
-                    (1 - math.fabs(Xpixel - lu_xy[0]))
+                + equation[2] * Yline
-                    * (1 - math.fabs(Yline - lu_xy[1]))
+                + equation[3]
-                )
-                + 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
-            pass
+            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
-        # dem_resolution = self._dem_resolution  # dem的精度
-        dem_resolution = 5
-        # TODO:设定为5m 1个点。
-        n = round(path_length / dem_resolution)
-        # n = round(path_length / 5)
-        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
-        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 _read_path_file(self):
+def to_line_coordination(point_x_s, point_y_s, point_x_e, point_y_e):
-        toml_dict = self._toml_dict
+    path_length = ((point_x_s - point_x_e) ** 2 + (point_y_s - point_y_e) ** 2) ** 0.5
-        path_file = toml_dict["parameter"]["path_file"]
+    n = round(path_length / 1)
-        excel_pds = pd.read_excel(path_file, header=None)
+    center_point_x = np.linspace(point_x_s, point_x_e, n, endpoint=True)
-        return excel_pds
+    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 _copy_db_file(self):
+
-        toml_dict = self._toml_dict
+def distance(a, b):
-        db_file_dir = toml_dict["parameter"]["db_file_dir"]
+    r = np.sum((a - b) ** 2) ** 0.5
-        out_dxf_file_dir = toml_dict["parameter"]["out_dxf_file_dir"]
+    return r
-        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")
+if __name__ == "__main__":
-        shutil.copy(f"{db_file_dir}/RULE.db", f"{out_dxf_file_dir}/RULE.db")
+    dem = Dem(r"d:\工程\金上线\DEM\四川-金上105-CGCS2000.tif")
-        shutil.copy(f"{db_file_dir}/WIRE.db", f"{out_dxf_file_dir}/WIRE.db")
+    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.")

dem_utils.py

@@ -1,6 +0,0 @@
-import numpy as np
-
-
-def distance(a, b):
-    r = np.sum((a - b) ** 2) ** 0.5
-    return r

main.py

@@ -1,13 +1,7 @@
-import sys
+
-from dem import Dem
+
-from loguru import logger
+
 if __name__ == "__main__":
-    if len(sys.argv) < 2:
+    import transformer.exporter as exporter
-        toml_file_path = r"db_山地750.toml"
+    exporter.main()
-    else:
+    print("PyCharm")
-        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()
-    print("Finished.")

nwed.py

@@ -1,182 +0,0 @@
-from xml.etree.ElementTree import Element
-from xml.etree.ElementTree import ElementTree
-from xml.etree.ElementTree import SubElement
-import attrs
-from typing import List, Tuple
-
-
-@attrs.define
-class Nwed:
-    # section_name = "002"
-    # licheng = "100"
-    # remarks_start = "002"
-    # pole_name_start = "4001"
-    # remarks_end = "003"
-    # pole_name_end = "4002"
-    # verticalExtent2 = "20"  # 边线宽度
-    section_name: str
-    licheng: str
-    remarks_start: str
-    pole_name_start: str
-    remarks_end: str
-    pole_name_end: str
-    verticalExtent2: str  # 边线宽度
-    mid_vec_list: List[List[str]]
-    right_vec_list: List[List[str]]
-    left_vec_list: List[List[str]]
-
-    def indent(self, elem, level=0):
-        i = "\n" + level * "\t"
-        if len(elem):
-            if not elem.text or not elem.text.strip():
-                elem.text = i + "\t"
-            if not elem.tail or not elem.tail.strip():
-                elem.tail = i
-            for elem in elem:
-                self.indent(elem, level + 1)
-            if not elem.tail or not elem.tail.strip():
-                elem.tail = i
-        else:
-            if level and (not elem.tail or not elem.tail.strip()):
-                elem.tail = i
-
-    def AddEleElevationCurve(
-        self, elevationCurves, id, vec_list, curve_id, curve_side=None
-    ):
-        if curve_side:
-            elevationCurve_attribs = {"ID": id, "CurveSide": curve_side}
-        else:
-            elevationCurve_attribs = {"ID": id}
-        elevationCurve = SubElement(
-            elevationCurves, "EleElevationCurve", attrib=elevationCurve_attribs
-        )  # 中线
-        SubElement(elevationCurve, "GrowthHeight").text = "0"
-        SubElement(elevationCurve, "LinkedSectionObjectID_Group")
-        SubElement(elevationCurve, "wallheight").text = "0"
-        segmentList = SubElement(elevationCurve, "SegmentList")
-        eleElevationCurveSegment = SubElement(
-            segmentList,
-            "EleElevationCurveSegment",
-            # TODO 这个ID其实依据中、边导线会变化。
-            attrib={
-                "ID": curve_id,
-                "StartPointClass": "NotSet",
-                "EndPointClass": "NotSet",
-            },
-        )
-        SubElement(eleElevationCurveSegment, "GrowthHeight").text = "0"
-        SubElement(eleElevationCurveSegment, "LinkedSectionObjectID_Group")
-        vec_list_element = SubElement(eleElevationCurveSegment, "VecList")
-        for vec in vec_list:
-            SubElement(vec_list_element, "Vector3", attrib={"xyz": ",".join(vec)})
-        vec_list_rendering_element = SubElement(
-            eleElevationCurveSegment, "VecListRendering"
-        )
-        for vec in vec_list:
-            SubElement(
-                vec_list_rendering_element, "Vector3", attrib={"xyz": ",".join(vec)}
-            )
-        return elevationCurve
-
-    def AddAnnotationElements(self, sectionData, pole_attrib, xyz):
-        pdmGraphicsObject = SubElement(
-            sectionData, "PdmGraphicsObject", attrib=pole_attrib
-        )
-        SubElement(pdmGraphicsObject, "GrowthHeight").text = "0"
-        SubElement(pdmGraphicsObject, "LinkedSectionObjectID_Group")
-        SubElement(pdmGraphicsObject, "Position", attrib={"xyz": xyz})
-
-    def write(self, xml_path):
-        root = Element(
-            "OverheadTransmissionLineDrawing",
-            attrib={
-                "Name": self.section_name,
-                "xmlns:xsi": "http://www.w3.org/2001/XMLSchema-instance",
-                "xmlns:xsd": "http://www.w3.org/2001/XMLSchema",
-                "Description": "",
-                "Version": "1.0",
-                "StartAngle": "0",
-                "EndAngle": "0",
-                "NextID": "9",
-            },
-        )
-        SubElement(root, "lstCrossingCable")
-        drawingFrames = SubElement(root, "DrawingFrames", attrib={"PaperSizeId": "0"})
-        SubElement(drawingFrames, "GrowthHeight").text = "0"
-        SubElement(drawingFrames, "LinkedSectionObjectID_Group")
-        frames = SubElement(drawingFrames, "Frames")
-        frame = SubElement(
-            frames,
-            "Frame",
-            attrib={
-                "Name": "完整图纸",
-                "SplitterPosition": "0",
-                "PaperOrientation": "Landscape",
-                "UserVerticalOffset": "0",
-            },
-        )
-        SubElement(frame, "GrowthHeight").text = "0"
-        SubElement(frame, "LinkedSectionObjectID_Group")
-        elevAdjSections = SubElement(frame, "ElevAdjSections")
-        SubElement(elevAdjSections, "GrowthHeight").text = "0"
-        SubElement(elevAdjSections, "LinkedSectionObjectID_Group")
-        sections = SubElement(elevAdjSections, "Sections")
-        SubElement(
-            sections,
-            "ElevationAdjSection",
-            attrib={"Position": "0", "ElevationAjusting": "0"},
-        )
-        SubElement(root, "GPSPointGroups")
-        SubElement(root, "Orthoimages")
-        SubElement(root, "GPSXianGao")
-        SubElement(root, "HouseSwing")
-        sectionData = SubElement(root, "SectionData")
-        SubElement(sectionData, "licheng").text = self.licheng
-        SubElement(sectionData, "lstlicheng")
-        SubElement(sectionData, "lstSctnRcd")
-        SubElement(sectionData, "ForestElevationCurves")
-        elevationCurves = SubElement(sectionData, "ElevationCurves")
-        self.AddEleElevationCurve(elevationCurves, "1", self.mid_vec_list, "4")
-        self.AddEleElevationCurve(elevationCurves, "2", self.left_vec_list, "5", "Left")
-        self.AddEleElevationCurve(
-            elevationCurves, "3", self.right_vec_list, "6", "Right"
-        )
-        annotationElements = SubElement(sectionData, "AnnotationElements")
-        self.AddAnnotationElements(
-            annotationElements,
-            {
-                "xsi:type": "ElePowerLinePole",
-                "ID": "7",
-                "CurveSide": "All",
-                "Remarks": self.remarks_start,
-                "PoleType": "AnglePole",
-                "LeftElevation": "0",
-                "RightElevation": "0",
-                "CableAngle": "0",
-                "PoleName": self.pole_name_start,
-            },
-            xyz=",".join(self.mid_vec_list[0]),
-        )
-        self.AddAnnotationElements(
-            annotationElements,
-            {
-                "xsi:type": "ElePowerLinePole",
-                "ID": "8",
-                "CurveSide": "All",
-                "Remarks": self.remarks_end,
-                "PoleType": "AnglePole",
-                "LeftElevation": "0",
-                "RightElevation": "0",
-                "CableAngle": "0",
-                "PoleName": self.pole_name_end,
-            },
-            xyz=",".join(self.mid_vec_list[-1]),
-        )
-        mapData = SubElement(root, "MapData")
-        SubElement(mapData, "GrowthHeight").text = "0"
-        SubElement(mapData, "LinkedSectionObjectID_Group")
-        SubElement(mapData, "verticalExtent2").text = self.verticalExtent2
-        SubElement(mapData, "MapElements")
-        tree = ElementTree(root)
-        self.indent(root)
-        tree.write(xml_path, encoding="utf-8", xml_declaration=True)

pw.py

@@ -1,102 +0,0 @@
-import numpy as np
-
-
-# D文件
-class DFile:
-    def __init__(
-        self,
-        start_num,
-        end_num,
-        mileage: np.ndarray,
-        center_z,
-        left_z,
-        right_z,
-        tree_height,
-    ):
-        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
-        self._tree_height = tree_height
-        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
-        if self._tree_height > 0:
-            mat = np.concatenate(
-                (
-                    np.array([mat[0, :]]),
-                    np.array([[197, 0, 0, 0, self._tree_height]]),
-                    mat[1:, :],
-                )
-            )  # 197
-            mat = np.concatenate(
-                (
-                    mat[0:-2, :],
-                    np.array([[198, mat[-1,1], 0, 0, self._tree_height]]),
-                    np.array([mat[-1, :]]),
-                ),
-            )  # 198
-        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}/PW{segments[0]-4000+1}.dat", "w") as c_file:
-        with open(f"{out_dir}/PW01.dat", "w") as c_file:
-            c_file.write(c_template)

transformer/__init__.py

@@ -12,9 +12,8 @@ class Projector:
 
 
 EPSG_WGS_84 = "epsg:4326"
-EPSG_CGCS_2000_105E = "epsg:4507"
+EPSG_Xian_1980_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")
 

transformer/exporter.py

@@ -16,15 +16,17 @@ def prettify(elem):
 
 
 def main():
-    df = pd.read_excel("杆塔成果表(左侧）-复测.xls")
+    df = pd.read_excel("成果表-0120.xlsx")
-    needed_df = df.iloc[:85, :]
+    needed_df = df.iloc[:120, :3]
-    trans = transformer.Projector(transformer.EPSG_Xian_1980_114E, transformer.EPSG_WGS_84)
+    trans = transformer.Projector(
-    y_array = needed_df["北坐标"].to_numpy()
+        transformer.EPSG_CGCS_2000_108E, transformer.EPSG_WGS_84
-    x_array = needed_df["东坐标"].to_numpy()
+    )
+    y_array = needed_df["北坐标(m)"].to_numpy()
+    x_array = needed_df["东坐标(m)"].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 = {
@@ -59,5 +61,5 @@ def main():
         "{lon},{lat},{elevation}".format(lon=foo[1], lat=foo[0], elevation=foo[2])
         for foo in line_coordination_list
     )
-    with open("周口.kml", "w", encoding="utf-8") as f:
+    with open("out.kml", "w", encoding="utf-8") as f:
         f.write(prettify(xml_root).decode("utf-8"))

transmission_line/plane.py

@@ -4,101 +4,26 @@ import ezdxf
 from data_types import WGS84Point
 import transformer
 import numpy as np
-from functools import singledispatch
-from typing import Tuple
 
 
-class Plane:
+def generate_wgs_84_point(point_start: WGS84Point, point_end: WGS84Point, n_point):
-    # dem文件坐标系统
+    _transformer_84_yx = transformer.Projector("epsg:4326", "epsg:4507")
-    def __init__(self, dem_system, plane_system, dem_path):
+    start_y, start_x = _transformer_84_yx.transform(
-        self._dem_path = dem_path
+        point_start.latitude, point_start.longitude
-        self._dem_system = dem_system
+    )
-        self._start_84_point = None
+    end_y, end_x = _transformer_84_yx.transform(point_end.latitude, point_end.longitude)
-        self._end_84_point = None
+    # 插入点
-        self._plane_system = plane_system
+    x_linspace = np.linspace(start_x, end_x, n_point)
-        self._transformer_84_to_dem_ = transformer.Projector(
+    y_linspace = np.linspace(start_y, end_y, n_point)
-            transformer.EPSG_WGS_84, dem_system
+    elevation_linspace = np.linspace(
-        )
+        point_start.elevation, point_end.elevation, n_point
-        self._transformer_dem_to_84 = transformer.Projector(
+    )
-            dem_system, transformer.EPSG_WGS_84
+    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(
-    @singledispatch
+        [np.array(_transformer_yx_84.transform(y, x)) for y, x, _elevation in yxz_point]
-    def set_line(self, start: WGS84Point, end: WGS84Point):
+    )
-        self._start_84_point = start
+    return wgs_84_point
-        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):
@@ -115,26 +40,28 @@ 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 main():
+# 计算档距
-    _plane = Plane(
+def cal_span(start: WGS84Point, end: WGS84Point):
-        transformer.EPSG_Xian_1980_114E, transformer.EPSG_Xian_1980_114E, r"d:\307.tif"
+    _transformer_84_yx = transformer.Projector("epsg:4326", "epsg:4507")
+    start_point = np.array(
+        _transformer_84_yx.transform(start.latitude, start.longitude)
     )
-    # start_point = WGS84Point(33.24643289, 114.58574497, 28)
+    end_point = np.array(_transformer_84_yx.transform(end.latitude, end.longitude))
-    # end_point = WGS84Point(33.22873308, 114.57951634, 28)
+    return np.round(np.sum((start_point - end_point) ** 2) ** 0.5)
-    _plane.set_line_yx((3680231.27, 554470.91), (3678264.93, 553901.32))
+
-    span = _plane.cal_span()
+
+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))
     print("档距 {span}".format(span=span))
-    elevation_in_line = _plane.get_elevation(int(span / 5))
+    wgs_84_point_in_line = generate_wgs_84_point(start_point, end_point, span)
-    # dem_gcs = dem.Dem()
+    elevation_in_line = dem_gcs.get_elevation(wgs_84_point_in_line)
-    # dem_gcs.get_dem_info(if_print=False)
+    draw_dxf(elevation_in_line, 1, "a.dxf")
-    #
-    # 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")

wgs84toxy.py

@@ -1,11 +0,0 @@
-import transformer
-import pandas as pd
-if __name__=='__main__':
-    proj=transformer.Projector(transformer.EPSG_WGS_84,transformer.EPSG_CGCS_2000_108E)
-    bzx=pd.read_excel(r'd:\3040\Desktop\宝中线.xlsx')
-    lat=bzx['纬度']
-    lon=bzx['经度']
-    y,x=proj.transform(lat.to_numpy(),lon.to_numpy())
-    xy_df=pd.DataFrame({'y':y,'x':x})
-    print(xy_df)
-    xy_df.to_excel(r'd:\3040\Desktop\宝中线xy.xlsx')