基于python实现利用DEM数据计算坡度、坡向

您所在的位置:网站首页 地形起伏和坡度 基于python实现利用DEM数据计算坡度、坡向

基于python实现利用DEM数据计算坡度、坡向

2024-07-12 00:55| 来源: 网络整理| 查看: 265

基本概念

DEM数据

DataMark:CNSDTF-DEM Version:1.0 Unit:M Alpha:0.000000 Compress:0.000000 X0:258000.000 Y0:324000.000 DX:22.500 DY:22.500 Row:321 Col:481 ValueType:Integer Hzoom:1000 192743 191068 187814 181388 173357 165286 157185 152266 147335 142391 139052 137327 135608 130633 124014 117372 118908 123706 128495 133265 138026 139540 139440 139337 135994 131025 127669 124312 120947 119208 125633 135280 141643 149598 155913 160612 162106 163600 160300 160203 158500 156795 155094 151780 148466 143534 138590 133632 130287 126935 125208 126733 133122 136259 137773 139280 142407 143916 142199 142101

坡度 在这里插入图片描述 坡度是法线与铅垂线之间的夹角。如图2.2所示,利用坡度计算公式 α = a r c t a n f x 2 + f y 2 \alpha=arctan \sqrt{f_x^2+f_y^2} α=arctanfx2​+fy2​ ​ ,进行求解。 坡向 在这里插入图片描述 坡向是法线在水平面上的投影与正北方向之间的夹角(顺时针度量)。如图2.3所示,利用坡向计算公式 A = arctan ⁡ f y f x A = \arctan \frac{{f_y }}{{f_x }} A=arctanfx​fy​​ , 求解。

程序流程

在这里插入图片描述 在这里插入图片描述 2)计算各点梯度 在这里插入图片描述 如图3.5所示,可以利用像下面的公式写一个函数进行计算, s l o p e x = e 1 − e 3 2 × d s i z e x          s l o p e y = e 4 − e 2 2 × d s i z e y slope_x = \frac{{e_1 - e_3 }}{{2 \times dsizex}}\;\;\;\;slope_y = \frac{{e_4 - e_2 }}{{2 \times dsizey}} slopex​=2×dsizexe1​−e3​​slopey​=2×dsizeye4​−e2​​ 。首先计算3*3DEM格网数据的X,Y方向的变化。对于整幅图而言,由于本图像元点不是特别多,故对于边缘的点处理时,可以在栅格点外再加一圈栅格便于求梯度。而添加点的方式只需要把原本边缘的栅格点的高程值复制一遍即可,如图3.6所示。 在这里插入图片描述 3)解算坡度、坡向 如图3.5所示,可以利用下面的公式计算求得坡度和坡向 s l o p = arctan ⁡ s l o p e x 2 + s l o p e y 2        A = arctan ⁡ s l o p e y s l o p e x slop = \arctan \sqrt {slope_x^2 + slope_y^2 } \;\;\;A = \arctan \frac{{slope_y }}{{slope_x }} slop=arctanslopex2​+slopey2​ ​A=arctanslopex​slopey​​。在这一过程中可以整体进行计算,从而较为方便地获得坡度和坡向的二维矩阵。 4)可视化 利用python中的绘图库matplotlib,编写了一个可视化函数,利用plot里面的imshow,将之前计算得到的二维矩阵存进去,最后以栅格的方式显示。对于三维DEM显示,利用Axes3D的方法,利用DEM文件中的起始坐标和分辨率推算出各个栅格的坐标,对应出高程值,绘制出三维图形。

程序代码 DEMclass模块 from mpl_toolkits.mplot3d import Axes3D from matplotlib import cbook from matplotlib import cm from matplotlib.colors import LightSource import matplotlib.pyplot as plt import numpy as np import re import math #####读取文件并处理成numpy并返回 def readfile(): with open('datas.DEM', 'r', encoding='utf-8')as file: datas = file.readlines()[13:] list1 = [] strs = "" row = 321 col = 481 npdata = np.zeros((row, col), dtype=np.int16) for data in datas: data = data.strip() if len(data) > 20: strs = strs + " " + data if len(data) 0.: a[i,j]=0. else: a[i,j]=180. elif y==0.: if x>0: a[i,j]=90. else: a[i,j]=270. else: a[i,j]=float(math.atan(y/x))*57.29578 if a[i,j]90.: a[i,j]=450.-a[i,j] else: a[i,j]=90.-a[i,j] return slope,a ####绘制平面栅格图 def Drawgrid(judge,pre=[],A=[],strs=""): if judge==0: if strs == "": plt.imshow(A, interpolation='nearest', cmap=plt.cm.hot, origin='lower') # cmap='bone' cmap=plt.cm.hot plt.colorbar(shrink=0.8) plt.xticks(()) plt.yticks(()) plt.show() else: plt.imshow(A, interpolation='nearest', cmap=strs, origin='lower') # cmap='bone' cmap=plt.cm.hot plt.colorbar(shrink=0.8) xt=range(258000, 268822,22) yt=range(324000, 331222,22) plt.xticks(()) plt.yticks(()) plt.show() elif judge==1: fig = plt.figure() ax = Axes3D(fig) # X = np.arange(1,482,1) # Y = np.arange(1,322,1) X = np.arange(258000.000, 268822.500, 22.5) Y = np.arange(324000.000, 331222.500, 22.5) X, Y = np.meshgrid(X, Y) Z = pre ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=plt.get_cmap('rainbow')) # cmap=plt.get_cmap('rainbow') ax.contourf(X, Y, Z, zdir='z', offset=-2, cmap=plt.cm.hot) ax.set_zlim(0, 200000) plt.show() 主函数 import DEMclass as dem from DEMclass import Drawgrid import numpy as np ####程序入口 if __name__=='__main__': npgrid=dem.readfile() pre=npgrid npgrid=dem.AddRound(npgrid) dx,dy=dem.Cacdxdy(npgrid,22.5,22.5) slope,arf=dem.CacSlopAsp(dx,dy) dem.np.savetxt("slope.csv",slope,delimiter=",") #绘制三维DEM Drawgrid(judge=1,pre=pre) #绘制二维DEM Drawgrid(judge=0,A=pre,strs="bone") #绘制坡度图 Drawgrid(judge=0,A=slope,strs="rainbow") #绘制坡向图 Drawgrid(judge=0,A=arf) 效果图

三维图: 在这里插入图片描述 二维DEM: 在这里插入图片描述 坡度图: 在这里插入图片描述 坡向图: 在这里插入图片描述

参考资料

1.Python的地形三维可视化——简介Matplotlib和gdal https://blog.csdn.net/allenlu2008/article/details/51880333 2.Pycharm中配置GDAL库的一种方法 https://blog.csdn.net/qq_35960361/article/details/96438650 3.python 调用gdal 处理dem数据 https://blog.csdn.net/qq_15642411/article/details/79123677 4.适用于Python扩展程序包的非官方Windows二进制文件——GDAL https://www.lfd.uci.edu/~gohlke/pythonlibs/#gdal 5.python gdal安装与简单使用 https://blog.csdn.net/nima1994/article/details/79207805/ pip install GDAL-2.4.1-cp37-cp37m-win_amd64.whl 6.(目测有用)【GDAL】python读取DEM计算坡度与坡向 https://blog.csdn.net/qq_37935516/article/details/85028979 7.arcgis计算坡度、坡向 http://help.arcgis.com/zh-cn/arcgisdesktop/10.0/help/index.html#//009z000000vz000000 8.C语言实现GDAL使用DEM数据计算坡度坡向 https://blog.csdn.net/liminlu0314/article/details/8498985 9.gdal在python功能(官方) https://gdal.org/python/index.html 10.玩转python的正则表达式|提取字符串中的所有数字 https://blog.csdn.net/guanyonglai/article/details/89512659 https://blog.csdn.net/u011436427/article/details/82628597 可视化 1.三维图(目测有用) https://blog.51cto.com/9291927/2435621 2.地图+热力图(可以用于坡度坡向) https://www.zhihu.com/question/33783546 3.统计 https://www.cnblogs.com/dudududu/p/9149762.html 4.turtle库 https://blog.csdn.net/weixin_44078196/article/details/10065651



【本文地址】


今日新闻


推荐新闻


CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3