坡度坡向分析方法 坡度(slope)是地面特定区域高度变化比率的量度。坡度的表示方法有百分比法、度数法、密位法和分数法四种,其中以百分比法和度数法较为常用。本文计算的为坡度百分比数据。如当角度为45度(弧度为π/4)时,高程增量等于水平增量,高程增量百分比为100%。 坡向(aspect)是指地形坡面的朝向。坡向用于识别出从每个像元到其相邻像元方向上值的变化率最大的下坡方向。坡向可以被视为坡度方向。坡向是一个角度,将按照顺时针方向进行测量,角度范围介于 0(正东)到 360(仍是正东)之间,即完整的圆。不具有下坡方向的平坦区域将赋值为-1(arcgis处理时为-1,其他可能为0)。 坡度、坡向计算一般采用拟合曲面法。拟合曲面一般采用二次曲面,即3×3的窗口,如下图所示。每个窗口的中心为一个高程点。图中的中心点e坡度和坡向计算过程如下。 参考链接: [1]https://blog.csdn.net/zhouxuguang236/article/details/40017219 [2]https://blog.csdn.net/weixin_45561357/article/details/106677574 [3]https://www.cnblogs.com/gispathfinder/p/5790469.html 注意:DEM的空间坐标系一定要为投影坐标系。 ArcGIS坡度坡向分析 打开DEM数据 坡度分析 坡度结果 坡向分析 坡向结果 python-gdal坡度坡向分析 from osgeo import gdal demfile = 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) 坡度结果 坡向结果 python坡度坡向分析 import gdal import numpy as np from scipy import ndimage as nd from copy import deepcopy demfile = r"D:\微信公众号\坡度坡向\N40E117_Albers.tif" slopefile = r"D:\微信公众号\坡度坡向\N40E117_Albers_python_Slope.tif" #读取DEM数据 ds = gdal.Open(demfile) cols = ds.RasterXSize rows = ds.RasterYSize geo = 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.nan mask = 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])) * 100 s_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) 坡度结果 import gdal import numpy as np from scipy import ndimage as nd from copy import deepcopy demfile = r"D:\微信公众号\坡度坡向\N40E117_Albers.tif" aspectfile = r"D:\微信公众号\坡度坡向\N40E117_Albers_python_Aspect.tif" #读取DEM数据 ds = gdal.Open(demfile) cols = ds.RasterXSize rows = ds.RasterYSize geo = 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.nan mask = 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.29578 a_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.)] = -1 a_data[1:-1, 1:-1][(x==0.)& (y>0.)] = 0 a_data[1:-1, 1:-1][(x==0.)& (y<0.)] = 180 a_data[1:-1, 1:-1][(x>0.)& (y==0.)] = 90 a_data[1:-1, 1:-1][(x<0.)& (y==0.)] = 270. # 边缘填充 a_mask = a_data==-999 if 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) 坡向结果 测试数据: 链接:https://pan.baidu.com/s/1PODbTJn1JOqOA4qeaJq4Gg 提取码:l3fw 欢迎关注个人wx_gzh: 小Rser |
|