坡度坡向分析方法

坡度(slope)是地面特定區(qū)域高度變化比率的量度。坡度的表示方法有百分比法、度數(shù)法、密位法和分?jǐn)?shù)法四種，其中以百分比法和度數(shù)法較為常用。本文計(jì)算的為坡度百分比數(shù)據(jù)。如當(dāng)角度為45度（弧度為π/4）時(shí)，高程增量等于水平增量，高程增量百分比為100%。

坡向(aspect)是指地形坡面的朝向。坡向用于識(shí)別出從每個(gè)像元到其相鄰像元方向上值的變化率最大的下坡方向。坡向可以被視為坡度方向。坡向是一個(gè)角度，將按照順時(shí)針?lè)较蜻M(jìn)行測(cè)量，角度范圍介于 0（正東）到 360（仍是正東）之間，即完整的圓。不具有下坡方向的平坦區(qū)域?qū)①x值為-1（arcgis處理時(shí)為-1，其他可能為0）。

坡度、坡向計(jì)算一般采用擬合曲面法。擬合曲面一般采用二次曲面，即3×3的窗口，如下圖所示。每個(gè)窗口的中心為一個(gè)高程點(diǎn)。圖中的中心點(diǎn)e坡度和坡向計(jì)算過(guò)程如下。

參考鏈接：

[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的空間坐標(biāo)系一定要為投影坐標(biāo)系。

ArcGIS坡度坡向分析

打開DEM數(shù)據(jù)

坡度分析

坡度結(jié)果

坡向分析

坡向結(jié)果

python-gdal坡度坡向分析

from?osgeo?import?gdal

demfile?=?r"D:\微信公眾號(hào)\坡度坡向\N40E117_Albers.tif"

#?獲取DEM信息
infoDEM?=?gdal.Info(demfile)

#?計(jì)算坡度
slopfile?=?r"D:\微信公眾號(hào)\坡度坡向\N40E117_Albers_gdal_Slope.tif"
slope?=?gdal.DEMProcessing(slopfile,?demfile,?"slope",?format='GTiff',?slopeFormat="percent",?zeroForFlat=1,?computeEdges=True)

#?計(jì)算坡向
aspectfile?=?r"D:\微信公眾號(hào)\坡度坡向\N40E117_Albers_gdal_Aspect.tif"
b?=?gdal.DEMProcessing(aspectfile,?demfile,?"aspect",?format='GTiff',?trigonometric=0,?zeroForFlat=1,?computeEdges=True)

坡度結(jié)果

坡向結(jié)果

python坡度坡向分析

import?gdal
import?numpy?as?np
from?scipy?import?ndimage?as?nd
from?copy?import?deepcopy

demfile?=?r"D:\微信公眾號(hào)\坡度坡向\N40E117_Albers.tif"
slopefile?=?r"D:\微信公眾號(hào)\坡度坡向\N40E117_Albers_python_Slope.tif"

#讀取DEM數(shù)據(jù)
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)
#?將無(wú)效值或背景值臨近像元填充
if?np.sum(mask)?>?0:
????ind?=?nd.distance_transform_edt(mask,?return_distances=False,?return_indices=True)
????data?=?data[tuple(ind)]

#?計(jì)算坡度
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
#?寫出結(jié)果
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)

坡度結(jié)果

import?gdal
import?numpy?as?np
from?scipy?import?ndimage?as?nd
from?copy?import?deepcopy

demfile?=?r"D:\微信公眾號(hào)\坡度坡向\N40E117_Albers.tif"
aspectfile?=?r"D:\微信公眾號(hào)\坡度坡向\N40E117_Albers_python_Aspect.tif"

#讀取DEM數(shù)據(jù)
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)
#?將無(wú)效值或背景值臨近像元填充
if?np.sum(mask)?>?0:
????ind?=?nd.distance_transform_edt(mask,?return_distances=False,?return_indices=True)
????data?=?data[tuple(ind)]

#?計(jì)算坡向
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
#?寫出結(jié)果
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)

坡向結(jié)果

測(cè)試數(shù)據(jù)：

鏈接：https://pan.baidu.com/s/1PODbTJn1JOqOA4qeaJq4Gg?

提取碼：l3fw?

歡迎關(guān)注個(gè)人wx_gzh: 小Rser