这个是虚坏叔叔的主站。牢记域名: [xuhss.com](拼音首字母)。首页可以搜索或: 虚幻账号

基于DEM的坡度坡向分析

🚀 优质资源分享 🚀

学习路线指引(点击解锁) 知识定位 人群定位
🧡 Python实战微信订餐小程序 🧡 进阶级 本课程是python flask+微信小程序的完美结合,从项目搭建到腾讯云部署上线,打造一个全栈订餐系统。
💛Python量化交易实战💛 入门级 手把手带你打造一个易扩展、更安全、效率更高的量化交易系统

坡度坡向分析方法

坡度(slope)是地面特定区域高度变化比率的量度。坡度的表示方法有百分比法、度数法、密位法和分数法四种,其中以百分比法和度数法较为常用。本文计算的为坡度百分比数据。如当角度为45度(弧度为π/4)时,高程增量等于水平增量,高程增量百分比为100%。

4f3ac20cd7c337c87f18186a6bb06ef7 - 基于DEM的坡度坡向分析

坡向(aspect)是指地形坡面的朝向。坡向用于识别出从每个像元到其相邻像元方向上值的变化率最大的下坡方向。坡向可以被视为坡度方向。坡向是一个角度,将按照顺时针方向进行测量,角度范围介于 0(正东)到 360(仍是正东)之间,即完整的圆。不具有下坡方向的平坦区域将赋值为-1(arcgis处理时为-1,其他可能为0)。

坡度、坡向计算一般采用拟合曲面法。拟合曲面一般采用二次曲面,即3×3的窗口,如下图所示。每个窗口的中心为一个高程点。图中的中心点e坡度和坡向计算过程如下。

7562d13016aeeb9ca70040ae7a9c3d2d - 基于DEM的坡度坡向分析

afbb58b9fde23fae8d292e8ad2da64b4 - 基于DEM的坡度坡向分析

5715a76d3a3263e9c12e4391c53246e0 - 基于DEM的坡度坡向分析

参考链接:

[1]https://blog.csdn.net/zhouxuguang236/article/details/40017219

[2]https://blog.csdn.net/weixin\_45561357/article/details/106677574

[3]https://blog.csdn.net/gispathfinder/p/5790469.html

注意:DEM的空间坐标系一定要为投影坐标系

ArcGIS坡度坡向分析

打开DEM数据

f1e72285dd5bd93750b3e1f8137778b0 - 基于DEM的坡度坡向分析

坡度分析

bfae9cf6225524ca6df8a2d4787c0476 - 基于DEM的坡度坡向分析

坡度结果

e7fee373095c9df8f63cb60abd7bdf86 - 基于DEM的坡度坡向分析

坡向分析

830ace88f7b70123d1a0a3c0985bb332 - 基于DEM的坡度坡向分析

坡向结果

ea09a7338c95a22760a95f5412600b47 - 基于DEM的坡度坡向分析

python-gdal坡度坡向分析

from osgeo import gdaldemfile = r"D:\微信公众号\坡度坡向\N40E117\_Albers.tif"# 获取DEM信息infoDEM = gdal.Info(demfile)# 计算坡度slopfile = r"D:\微信公众号\坡度坡向\N40E117\_Albers\_gdal\_Slope.tif"slope = gdal.DEMProcessing(slopfile, demfile, "slope", format='GTiff', slopeFormat="percent", zeroForFlat=1, computeEdges=True)# 计算坡向aspectfile = r"D:\微信公众号\坡度坡向\N40E117\_Albers\_gdal\_Aspect.tif"b = gdal.DEMProcessing(aspectfile, demfile, "aspect", format='GTiff', trigonometric=0, zeroForFlat=1, computeEdges=True)

坡度结果

fc6ee04ded76dd57072b519fac8cc655 - 基于DEM的坡度坡向分析

坡向结果

db6afe06e1671eac47dd00a50eee8726 - 基于DEM的坡度坡向分析

python坡度坡向分析

import gdalimport numpy as npfrom scipy import ndimage as ndfrom copy import deepcopydemfile = r"D:\微信公众号\坡度坡向\N40E117\_Albers.tif"slopefile = r"D:\微信公众号\坡度坡向\N40E117\_Albers\_python\_Slope.tif"#读取DEM数据ds = gdal.Open(demfile)cols = ds.RasterXSizerows = ds.RasterYSizegeo = ds.GetGeoTransform()proj = ds.GetProjection()dem\_data = ds.ReadAsArray()data = deepcopy(dem\_data).astype(np.float32)band = ds.GetRasterBand(1)nodata = band.GetNoDataValue()data[data == nodata] = np.nan# data[data<-999]=np.nanmask = np.isnan(data)# 将无效值或背景值临近像元填充if np.sum(mask) > 0:    ind = nd.distance\_transform\_edt(mask, return\_distances=False, return\_indices=True)    data = data[tuple(ind)]# 计算坡度xsize = np.abs(geo[1])ysize = np.abs(geo[5])x = ((data[:-2, 2:] - data[:-2, :-2]) + 2 * (data[1:-1, 2:] - data[1:-1, :-2]) + (data[2:, 2:] - data[2:, :-2])) / (8 * xsize)y = ((data[2:, :-2] - data[:-2, :-2]) + 2 * (data[2:, 1:-1] - data[:-2, 1:-1]) + (data[2:, 2:] - data[:-2, 2:])) / (8 * ysize)s\_data = np.full((rows, cols), -999, dtype=np.float32)s\_data[1:-1, 1:-1] = (np.arctan(np.sqrt((np.power(x, 2) + np.power(y, 2)))))s\_data[1:-1, 1:-1] = np.abs(np.tan(s\_data[1:-1, 1:-1])) * 100s\_mask = s\_data==-999# 边缘填充if np.sum(s\_mask) > 0:    ind = nd.distance\_transform\_edt(s\_mask, return\_distances=False, return\_indices=True)    s\_data = s\_data[tuple(ind)]# 掩膜s\_data[dem\_data==nodata] = -999# 写出结果driver = gdal.GetDriverByName("gtiff")outds = driver.Create(slopefile, cols, rows, 1, gdal.GDT\_Float32)outds.SetGeoTransform(geo)outds.SetProjection(proj)outband = outds.GetRasterBand(1)outband.WriteArray(s\_data)outband.SetNoDataValue(-999)

坡度结果

2b7439240364ccef39d58cace4b240d5 - 基于DEM的坡度坡向分析

import gdalimport numpy as npfrom scipy import ndimage as ndfrom copy import deepcopydemfile = r"D:\微信公众号\坡度坡向\N40E117\_Albers.tif"aspectfile = r"D:\微信公众号\坡度坡向\N40E117\_Albers\_python\_Aspect.tif"#读取DEM数据ds = gdal.Open(demfile)cols = ds.RasterXSizerows = ds.RasterYSizegeo = ds.GetGeoTransform()proj = ds.GetProjection()dem\_data = ds.ReadAsArray()data = deepcopy(dem\_data).astype(np.float32)band = ds.GetRasterBand(1)nodata = band.GetNoDataValue()data[data == nodata] = np.nan# data[data<-999]=np.nanmask = np.isnan(data)# 将无效值或背景值临近像元填充if np.sum(mask) > 0:    ind = nd.distance\_transform\_edt(mask, return\_distances=False, return\_indices=True)    data = data[tuple(ind)]# 计算坡向xsize = np.abs(geo[1])ysize = np.abs(geo[5])x = ((data[:-2, 2:] - data[:-2, :-2]) + 2 * (data[1:-1, 2:] - data[1:-1, :-2]) + (data[2:, 2:] - data[2:, :-2])) / (8 * xsize)y = ((data[2:, :-2] - data[:-2, :-2]) + 2 * (data[2:, 1:-1] - data[:-2, 1:-1]) + (data[2:, 2:] - data[:-2, 2:])) / (8 * ysize)a\_data = np.full((rows, cols), -999, dtype=np.float32)a\_data[1:-1, 1:-1] = np.arctan2(y, -1 * x) * 57.29578a\_data\_ = deepcopy(a\_data[1:-1, 1:-1])a\_data[1:-1, 1:-1][a\_data\_ < 0] = 90 - a\_data[1:-1, 1:-1][a\_data\_ < 0]a\_data[1:-1, 1:-1][a\_data\_ >90] = 450 - a\_data[1:-1, 1:-1][a\_data\_ > 90]a\_data[1:-1, 1:-1][(a\_data\_ >= 0) & (a\_data\_ <= 90)] = 90 - a\_data[1:-1, 1:-1][(a\_data\_ >= 0) & (a\_data\_ <= 90)]a\_data[1:-1, 1:-1][(x==0.)& (y==0.)] = -1a\_data[1:-1, 1:-1][(x==0.)& (y>0.)] = 0a\_data[1:-1, 1:-1][(x==0.)& (y<0.)] = 180a\_data[1:-1, 1:-1][(x>0.)& (y==0.)] = 90a\_data[1:-1, 1:-1][(x<0.)& (y==0.)] = 270.# 边缘填充a\_mask = a\_data==-999if np.sum(a\_mask) > 0:    ind = nd.distance\_transform\_edt(a\_mask, return\_distances=False, return\_indices=True)    a\_data = a\_data[tuple(ind)]# 掩膜a\_data[dem\_data==nodata] = -999# 写出结果driver = gdal.GetDriverByName("gtiff")outds = driver.Create(aspectfile, cols, rows, 1, gdal.GDT\_Float32)outds.SetGeoTransform(geo)outds.SetProjection(proj)outband = outds.GetRasterBand(1)outband.WriteArray(a\_data)outband.SetNoDataValue(-999)

坡向结果

dd9a2296e7cdea39e1e98c1cfc620f91 - 基于DEM的坡度坡向分析

测试数据:

链接:https://pan.baidu.com/s/1PODbTJn1JOqOA4qeaJq4Gg

提取码:l3fw

欢迎关注个人wx_gzh: 小Rser

转载请注明:虚坏叔叔 » 基于DEM的坡度坡向分析

喜欢 (0)

您必须 登录 才能发表评论!