重构了DEM模块，dem只有核心功能，无业务功能。
This commit is contained in:
n3040
parent
abd0ca8aa7
commit
bece894a0c
205
dem.py
205
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 <osgeo.gdal.Dataset> -- 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
|
||||
|
|
|
|||
|
|
@ -0,0 +1,6 @@
|
|||
import numpy as np
|
||||
|
||||
|
||||
def distance(a, b):
|
||||
r = np.sum((a - b) ** 2) ** 0.5
|
||||
return r
|
||||
|
|
@ -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.")
|
||||
Loading…
Reference in New Issue