# 基于DEM的坡度坡向分析

[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

ArcGIS坡度坡向分析

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)`

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)`

`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)`