完成基本功能。
This commit is contained in:
commit
f3a1dda9a0
|
|
@ -0,0 +1,5 @@
|
|||
class WGS84Point:
|
||||
def __init__(self, latitude, longitude, elevation):
|
||||
self.latitude = latitude
|
||||
self.longitude = longitude
|
||||
self.elevation = elevation
|
||||
|
|
@ -0,0 +1,78 @@
|
|||
from osgeo import gdal
|
||||
import numpy as np
|
||||
|
||||
|
||||
class Dem:
|
||||
def __init__(self, filepath):
|
||||
self._dataset = gdal.Open(filepath)
|
||||
|
||||
def get_dem_info(self, if_print=False):
|
||||
"""Get the information of DEM data.
|
||||
Parameters:
|
||||
dem_data <osgeo.gdal.Dataset> -- The data of DEM.
|
||||
if_print <bool> -- If print the information of DEM. Default is False (not print).
|
||||
Return:
|
||||
... <...> -- The information and parameters of DEM.
|
||||
"""
|
||||
dem_data = self._dataset
|
||||
dem_row = dem_data.RasterYSize # height
|
||||
dem_col = dem_data.RasterXSize # width
|
||||
dem_band = dem_data.RasterCount
|
||||
dem_gt = dem_data.GetGeoTransform()
|
||||
dem_proj = dem_data.GetProjection()
|
||||
|
||||
if if_print:
|
||||
print("\nThe information of DEM:")
|
||||
print("The number of row (height) is: %d" % dem_row)
|
||||
print("The number of column (width) is: %d" % dem_col)
|
||||
print("The number of band is: %d" % dem_band)
|
||||
print("The 6 GeoTransform parameters are:\n", dem_gt)
|
||||
print("The GCS/PCS information is:\n", dem_proj)
|
||||
|
||||
return dem_row, dem_col, dem_band, dem_gt, dem_proj
|
||||
|
||||
def get_elevation(self, site_latlng):
|
||||
"""Get the elevation of given locations from DEM in GCS.
|
||||
Parameters:
|
||||
dem_gcs <osgeo.gdal.Dataset> -- The input DEM (in GCS).
|
||||
site_latlng <numpy.ndarray> -- The latitude and longitude of given locations.
|
||||
Return:
|
||||
site_ele <numpy.ndarray> -- The elevation and other information of given locations.
|
||||
"""
|
||||
gdal_data = self._dataset
|
||||
|
||||
# gdal_band = gdal_data.GetRasterBand(1)
|
||||
# nodataval = gdal_band.GetNoDataValue()
|
||||
# if np.any(gdal_array == nodataval):
|
||||
# gdal_array[gdal_array == nodataval] = np.nan
|
||||
|
||||
gt = gdal_data.GetGeoTransform()
|
||||
# print("\nThe 6 GeoTransform parameters of DEM are:\n", gt)
|
||||
|
||||
N_site = site_latlng.shape[0]
|
||||
Xgeo = site_latlng[:, 1] # longitude
|
||||
Ygeo = site_latlng[:, 0] # latitude
|
||||
site_ele = np.zeros(N_site)
|
||||
for i in range(N_site):
|
||||
# Note:
|
||||
# Xgeo = gt[0] + Xpixel * gt[1] + Yline * gt[2]
|
||||
# Ygeo = gt[3] + Xpixel * gt[4] + Yline * gt[5]
|
||||
#
|
||||
# Xpixel - Pixel/column of DEM
|
||||
# Yline - Line/row of DEM
|
||||
#
|
||||
# Xgeo - Longitude
|
||||
# Ygeo - Latitude
|
||||
#
|
||||
# [0] = Longitude of left-top pixel
|
||||
# [3] = Latitude of left-top pixel
|
||||
#
|
||||
# [1] = + Pixel width
|
||||
# [5] = - Pixel height
|
||||
#
|
||||
# [2] = 0 for north up DEM
|
||||
# [4] = 0 for north up DEM
|
||||
Xpixel = int(round((Xgeo[i] - gt[0]) / gt[1]))
|
||||
Yline = int(round((Ygeo[i] - gt[3]) / gt[5]))
|
||||
site_ele[i] = gdal_data.ReadAsArray(Xpixel, Yline, 1, 1).astype(np.float)
|
||||
return site_ele
|
||||
|
|
@ -0,0 +1,34 @@
|
|||
from pyproj import Transformer
|
||||
|
||||
|
||||
class Projector:
|
||||
def __init__(self, from_epsg, to_epsg):
|
||||
self._transformer = Transformer.from_crs(from_epsg, to_epsg)
|
||||
|
||||
def transform(self, y_lat, x_lon):
|
||||
transformer = self._transformer
|
||||
lat_y, lon_x = transformer.transform(y_lat, x_lon)
|
||||
return lat_y, lon_x
|
||||
|
||||
|
||||
EPSG_WGS_84 = "epsg:4326"
|
||||
EPSG_Xian_1980_105E = "epsg:4507"
|
||||
EPSG_CGCS_2000_108E = "epsg:4545"
|
||||
# _transformer_84_yx = Transformer.from_crs("epsg:4326", "epsg:4507")
|
||||
# _transformer_yx_84 = Transformer.from_crs("epsg:4507", "epsg:4326")
|
||||
|
||||
# WGS 84 坐标转平面
|
||||
# latitude 纬度
|
||||
# longitude 经度
|
||||
|
||||
|
||||
# def wgs_84_yx(latitude, longitude, _elevation):
|
||||
# y, x = _transformer_84_yx.transform(latitude, longitude) # lon 和lat 可以是元组
|
||||
# return y, x
|
||||
#
|
||||
#
|
||||
# def yx_wgs_84(y, x, _elevation):
|
||||
# latitude, longitude = _transformer_yx_84.transform(y, x) # lon 和lat 可以是元组
|
||||
# return latitude, longitude
|
||||
#
|
||||
# def yx_2000_wgs_84(y, x, _elevation):
|
||||
|
|
@ -0,0 +1,66 @@
|
|||
from xml.dom import minidom
|
||||
import numpy as np
|
||||
import transformer
|
||||
import pandas as pd
|
||||
import xml.etree.ElementTree as ET
|
||||
|
||||
|
||||
def prettify(elem):
|
||||
"""Return a pretty-printed XML string for the Element.
|
||||
"""
|
||||
rough_string = ET.tostring(
|
||||
elem, encoding="utf-8", method="xml", xml_declaration=True
|
||||
)
|
||||
reparsed = minidom.parseString(rough_string)
|
||||
return reparsed.toprettyxml(indent="\t")
|
||||
|
||||
|
||||
def main():
|
||||
df = pd.read_excel("成果表-0120.xlsx")
|
||||
needed_df = df.iloc[:120, :3]
|
||||
# print(needed_df)
|
||||
trans = transformer.Projector(
|
||||
transformer.EPSG_CGCS_2000_108E, transformer.EPSG_WGS_84
|
||||
)
|
||||
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]}
|
||||
)
|
||||
xml_f = ET.parse("sgt.kml")
|
||||
ET.register_namespace("", "http://www.opengis.net/kml/2.2")
|
||||
root = xml_f.getroot()
|
||||
ns = {
|
||||
"NS": r"http://www.opengis.net/kml/2.2",
|
||||
}
|
||||
folders_has_placemark = root.findall(r"NS:Document//NS:Folder[NS:Placemark]", ns)
|
||||
point_folder = None
|
||||
for fld in folders_has_placemark:
|
||||
if fld.find(r"NS:name", ns).text == "点":
|
||||
point_folder = fld
|
||||
break
|
||||
# remove all placemark
|
||||
# for child in point_folder:
|
||||
# if child.tag == "{" + ns["NS"] + "}Placemark":
|
||||
# point_folder.remove(child)
|
||||
point_folder.clear()
|
||||
for _, row in name_lat_lon_df.iterrows():
|
||||
name = row["桩号"]
|
||||
lon = row["经度"]
|
||||
lat = row["纬度"]
|
||||
ele_placemark = ET.SubElement(point_folder, "Placemark")
|
||||
ele_name = ET.SubElement(ele_placemark, "name")
|
||||
ele_name.text = name
|
||||
ele_point = ET.SubElement(ele_placemark, "Point")
|
||||
coordinates = ET.SubElement(ele_point, "coordinates")
|
||||
coordinates.text = "{lon},{lat},0".format(lon=lon, lat=lat)
|
||||
# xml_f.write("out.kml",encoding='utf-8',method='xml',pretty_print=True)
|
||||
# print(xml_f.toString())
|
||||
with open("out.kml", "w") as f:
|
||||
f.write(prettify(root))
|
||||
|
||||
# with open('84.txt','w') as f:
|
||||
# for point in wgs_84_point:
|
||||
# f.write('{point[1]},{point[0]},0 \n'.format(point=point))
|
||||
|
|
@ -0,0 +1,67 @@
|
|||
from scipy import interpolate
|
||||
import dem
|
||||
import ezdxf
|
||||
from data_types import WGS84Point
|
||||
import transformer
|
||||
import numpy as np
|
||||
|
||||
|
||||
def generate_wgs_84_point(point_start: WGS84Point, point_end: WGS84Point, n_point):
|
||||
_transformer_84_yx = transformer.Projector("epsg:4326", "epsg:4507")
|
||||
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("epsg:4507", "epsg:4326")
|
||||
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 cal_span(start: WGS84Point, end: WGS84Point):
|
||||
_transformer_84_yx = transformer.Projector("epsg:4326", "epsg:4507")
|
||||
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 np.round(np.sum((start_point - end_point) ** 2) ** 0.5)
|
||||
|
||||
|
||||
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))
|
||||
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, 1, "a.dxf")
|
||||
Loading…
Reference in New Issue