141 lines
5.2 KiB
Python
141 lines
5.2 KiB
Python
from scipy import interpolate
|
|
import dem
|
|
import ezdxf
|
|
from data_types import WGS84Point
|
|
import transformer
|
|
import numpy as np
|
|
from functools import singledispatch
|
|
from typing import Tuple
|
|
|
|
|
|
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):
|
|
doc = ezdxf.new(dxfversion="R2004")
|
|
msp = doc.modelspace()
|
|
x_axis = (
|
|
np.linspace(0, len(elevation_in_line) * interval, len(elevation_in_line)) / 5
|
|
)
|
|
points = np.array([x_axis, (elevation_in_line - elevation_in_line[0]) * 2])
|
|
new_x_axis = (
|
|
np.linspace(
|
|
0, len(elevation_in_line) * interval, int(len(elevation_in_line) / 12.5)
|
|
)
|
|
/ 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, y) for x, y in points.transpose()])
|
|
doc.saveas(dxf_path)
|
|
|
|
|
|
def main():
|
|
_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))
|
|
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")
|