From f3a1dda9a0bfc1838ea8b49e4b0fabe454647692 Mon Sep 17 00:00:00 2001 From: facat Date: Sun, 17 Jan 2021 20:05:01 +0800 Subject: [PATCH] =?UTF-8?q?=E5=AE=8C=E6=88=90=E5=9F=BA=E6=9C=AC=E5=8A=9F?= =?UTF-8?q?=E8=83=BD=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- data_types.py | 5 +++ dem.py | 78 ++++++++++++++++++++++++++++++++++++++ main.py | 6 +++ transformer/__init__.py | 34 +++++++++++++++++ transformer/exporter.py | 66 ++++++++++++++++++++++++++++++++ transmission_line/plane.py | 67 ++++++++++++++++++++++++++++++++ 6 files changed, 256 insertions(+) create mode 100644 data_types.py create mode 100644 dem.py create mode 100644 main.py create mode 100644 transformer/__init__.py create mode 100644 transformer/exporter.py create mode 100644 transmission_line/plane.py diff --git a/data_types.py b/data_types.py new file mode 100644 index 0000000..136a9f3 --- /dev/null +++ b/data_types.py @@ -0,0 +1,5 @@ +class WGS84Point: + def __init__(self, latitude, longitude, elevation): + self.latitude = latitude + self.longitude = longitude + self.elevation = elevation diff --git a/dem.py b/dem.py new file mode 100644 index 0000000..b8558d3 --- /dev/null +++ b/dem.py @@ -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 -- The data of DEM. + if_print -- 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 -- The input DEM (in GCS). + site_latlng -- The latitude and longitude of given locations. + Return: + site_ele -- 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 diff --git a/main.py b/main.py new file mode 100644 index 0000000..41ff204 --- /dev/null +++ b/main.py @@ -0,0 +1,6 @@ + + + +if __name__ == "__main__": + # main() + print("PyCharm") diff --git a/transformer/__init__.py b/transformer/__init__.py new file mode 100644 index 0000000..ab50b67 --- /dev/null +++ b/transformer/__init__.py @@ -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): diff --git a/transformer/exporter.py b/transformer/exporter.py new file mode 100644 index 0000000..490c394 --- /dev/null +++ b/transformer/exporter.py @@ -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)) diff --git a/transmission_line/plane.py b/transmission_line/plane.py new file mode 100644 index 0000000..b42b219 --- /dev/null +++ b/transmission_line/plane.py @@ -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")