From e56399cccb84e7345a64b70838adbeb74fb52afa Mon Sep 17 00:00:00 2001 From: n3040 Date: Sat, 15 Jan 2022 23:42:49 +0800 Subject: [PATCH] =?UTF-8?q?=E5=A2=9E=E5=8A=A0=E4=BA=86=E5=9F=BA=E4=BA=8E?= =?UTF-8?q?=E5=91=A8=E5=9B=B4=E5=87=A0=E4=B8=AA=E7=82=B9=E5=B7=AE=E5=80=BC?= =?UTF-8?q?=E6=B1=82=E9=AB=98=E7=A8=8B=E7=9A=84=E5=8A=9F=E8=83=BD=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- dem.py | 106 ++++++++++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 101 insertions(+), 5 deletions(-) diff --git a/dem.py b/dem.py index b8558d3..3d03ffc 100644 --- a/dem.py +++ b/dem.py @@ -1,3 +1,6 @@ +import math + +import ezdxf from osgeo import gdal import numpy as np @@ -50,8 +53,8 @@ class Dem: # 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 + Xgeo = site_latlng[:, 0] + Ygeo = site_latlng[:, 1] site_ele = np.zeros(N_site) for i in range(N_site): # Note: @@ -72,7 +75,100 @@ class Dem: # # [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) + 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_elevation = gdal_data.ReadAsArray( + int(lu_xy[0]), int(lu_xy[1]), 1, 1 + ).astype(float) + ld_elevation = gdal_data.ReadAsArray( + int(ld_xy[0]), int(ld_xy[1]), 1, 1 + ).astype(float) + ru_elevation = gdal_data.ReadAsArray( + int(ru_xy[0]), int(ru_xy[1]), 1, 1 + ).astype(float) + rd_elevation = gdal_data.ReadAsArray( + int(rd_xy[0]), int(rd_xy[1]), 1, 1 + ).astype(float) + # delta_x = (Xpixel - ld_xy[0]) / 1 # 距离是1 + # delta_y = (ld_xy[1]-Yline) / 1 + # site_ele[i] = ( + # ld_elevation + # + (lu_elevation - ld_elevation) * delta_y + # + (rd_elevation - ld_elevation) * delta_x + # + (ld_elevation - lu_elevation + ru_elevation - rd_elevation) + # * delta_x + # * delta_y + # ) + + # 通过空间平面方程乃拟合Z值 + equa_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], + ] + ) + equa_b = np.array( + [lu_elevation[0], ld_elevation[0], ru_elevation[0], rd_elevation[0]] + ) + equa_c = np.linalg.solve(equa_a, equa_b) + point_z = ( + Xpixel * equa_c[0] + + equa_c[1] * Xpixel * Yline + + equa_c[2] * Yline + + equa_c[3] + ) + site_ele[i] = point_z + print(f"row:{Yline} col:{Xpixel} elev:{site_ele[i]}") return site_ele + + +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) + point_x = np.linspace(point_x_s, point_x_e, n, endpoint=True) + point_y = np.linspace(point_y_s, point_y_e, n, endpoint=True) + r = np.array([point_x, point_y]).transpose() + return r + + +def distance(a, b): + r = np.sum((a - b) ** 2) ** 0.5 + return r + + +if __name__ == "__main__": + dem = Dem(r"e:\西北院\天都山750kV线路\天都山-妙岭卫片\ASC\tianmiaoN03-dem.tif") + dem.get_dem_info(if_print=True) + line_coordination = to_line_coordination( + 573284.886, 4104898.933, 574543.845, 4105825.047 + ) + elevation = dem.get_elevation(line_coordination) + + # fdsfd = np.array([ + # [573284.886, 4104898.933], + # [573286.584, 4104900.182], + # [573288.281, 4104901.431], + # [573289.979 , 4104902.68], + # [573291.677 , 4104903.929], + # [573293.374 , 4104905.178], + # + # ]) + # elevation = dem.get_elevation(elevation) + 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, elevation[i] * 2) for i in range(len(elevation))] + ) + doc.saveas("d.dxf") + print("Finished.")