dem.py

@@ -100,11 +100,11 @@ class Dem:
                 dxfattribs={"color": 5},
             )  # 蓝色
             # 树的线
-            # 考虑用最高边线的情况
+            # TODO 没有考虑用最高边线的情况
             if self._tree_height > 0:
                 dm_msp.add_polyline2d(
                     [
-                        (x_axis[i] / 5, ( np.max( (center_elevation[i],left_elevation[i],right_elevation[i])) + self._tree_height) * 2)
+                        (x_axis[i] / 5, (center_elevation[i] + self._tree_height) * 2)
                         for i in range(len(center_elevation))
                     ]
                 )
@@ -130,7 +130,7 @@ class Dem:
                 z_file.write("0     ")
                 z_file.write("0     ")
                 z_file.write(f"{center_elevation[0]*2-50}")
-            # 写平面文件
+            # TODO:写平面文件
             plate_msp.add_polyline2d(
                 line_coordination[:, 0:2],
                 dxfattribs={"color": 1},
@@ -217,6 +217,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.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)])  # 右上
@@ -233,7 +237,36 @@ 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_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 = (
                 lu_elevation[0]
                 * (

main_dem.py

@@ -0,0 +1,13 @@
+import sys
+from dem import Dem
+from loguru import logger
+if __name__ == "__main__":
+    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.")