facat
/
gdal
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

优化了dem插值部分代码。

This commit is contained in:
n3040 2022-11-30 20:53:49 +08:00
parent ed76be5f7c
commit 3864d128d0
2 changed files with 77 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

105
dem.py
View File

@ -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):

16
main.py
View File

@ -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.")