建站
咨詢
售后
建站咨詢
新聞中心
坡度坡向分析方法
為科爾沁左翼等地區(qū)用戶提供了全套網(wǎng)頁設(shè)計(jì)制作服務(wù)，及科爾沁左翼網(wǎng)站建設(shè)行業(yè)解決方案。主營業(yè)務(wù)為做網(wǎng)站、網(wǎng)站制作、科爾沁左翼網(wǎng)站設(shè)計(jì)，以傳統(tǒng)方式定制建設(shè)網(wǎng)站，并提供域名空間備案等一條龍服務(wù)，秉承以專業(yè)、用心的態(tài)度為用戶提供真誠的服務(wù)。我們深信只要達(dá)到每一位用戶的要求，就會得到認(rèn)可，從而選擇與我們長期合作。這樣，我們也可以走得更遠(yuǎn)！
坡度(slope)是地面特定區(qū)域高度變化比率的量度。坡度的表示方法有百分比法、度數(shù)法、密位法和分?jǐn)?shù)法四種，其中以百分比法和度數(shù)法較為常用。本文計(jì)算的為坡度百分比數(shù)據(jù)。如當(dāng)角度為45度（弧度為π/4）時，高程增量等于水平增量，高程增量百分比為100%。
坡向(aspect)是指地形坡面的朝向。坡向用于識別出從每個像元到其相鄰像元方向上值的變化率最大的下坡方向。坡向可以被視為坡度方向。坡向是一個角度，將按照順時針方向進(jìn)行測量，角度范圍介于 0（正東）到 360（仍是正東）之間，即完整的圓。不具有下坡方向的平坦區(qū)域?qū)①x值為-1（arcgis處理時為-1，其他可能為0）。
坡度、坡向計(jì)算一般采用擬合曲面法。擬合曲面一般采用二次曲面，即3×3的窗口，如下圖所示。每個窗口的中心為一個高程點(diǎn)。圖中的中心點(diǎn)e坡度和坡向計(jì)算過程如下。
參考鏈接：
[1]https://blog.csdn.net/zhouxuguang236/article/details/
[2]https://blog.csdn.net/weixin_/article/details/
[3]https://www.cnblogs.com/gispathfinder/p/.html
注意：DEM的空間坐標(biāo)系一定要為投影坐標(biāo)系。
ArcGIS坡度坡向分析
打開DEM數(shù)據(jù)
坡度分析
坡度結(jié)果
坡向分析
坡向結(jié)果
python-gdal坡度坡向分析
from osgeo import gdal
demfile = r"D:\微信公眾號\坡度坡向\N40E117_Albers.tif"
# 獲取DEM信息
infoDEM = gdal.Info(demfile)
# 計(jì)算坡度
slopfile = r"D:\微信公眾號\坡度坡向\N40E117_Albers_gdal_Slope.tif"
slope = gdal.DEMProcessing(slopfile, demfile, "slope", format='GTiff', slopeFormat="percent", zeroForFlat=1, computeEdges=True)
# 計(jì)算坡向
aspectfile = r"D:\微信公眾號\坡度坡向\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:\微信公眾號\坡度坡向\N40E117_Albers.tif"
slopefile = r"D:\微信公眾號\坡度坡向\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)
# 將無效值或背景值臨近像元填充
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:\微信公眾號\坡度坡向\N40E117_Albers.tif"
aspectfile = r"D:\微信公眾號\坡度坡向\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)
# 將無效值或背景值臨近像元填充
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.
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é)果
測試數(shù)據(jù)：
鏈接：https://pan.baidu.com/s/1PODbTJn1JOqOA4qeaJq4Gg
提取碼：l3fw
?
歡迎關(guān)注個人wx_gzh: 小Rser
網(wǎng)站名稱：基于DEM的坡度坡向分析
網(wǎng)頁地址：http://www.ef60e0e.cn/article/dsogpoe.html
