根据dem计算坡度是地理信息系统中的基础操作,获得准确的坡度数据是进行相关表面分析的前提。为了实现批量化的坡度计算,我利用Python的Numpy库来实现。
1.导入相关库
from osgeo import gdal
import numpy as np
gdal库不能直接import,会报错。
如果没有安装以上库,请记得pip install 或者conda install
2. 设计坡度计算函数
坡度计算是通过九宫格来计算的,原理简单来说就是X方向和Y方向的坡度的平方和开根号。
代码如下:
def cal_slope(filename, outname):
'''
1. 设置一个比data数组大一圈的数组,将data赋值到该数组
2. 掩膜,设置无效数据为nan
nan(NAN,Nan):not a number 表示不是一个数字,np.nan是一个float类型的数据.因为nan不是一个数,所以相关计算都无法得到数字。
所有涉及nan的操作,返回的都是nan
np.where(condition, x, y)返回值是一个和condition的shape相同的numpy数组;当满足条件condition时,返回值中的元素从x中取,否则从y中取
3. 取出符合条件的ei(1<=i<=8)
'''
ds = gdal.Open(filename)
data = ds.ReadAsArray()
nrows, ncols = data.shape
trasnformation = ds.GetGeoTransform()
cell_size = trasnformation[1]
#扩充生成一个比data数组大一圈的数组
extend_arr = np.full((nrows + 2, ncols + 2), -9999.)
extend_arr[1:nrows + 1, 1:ncols + 1] = data
extend_arr[np.where(extend_arr == -9999.)] = np.nan
# 不用双重循环,取出ei(1<=i<=8)的值,省略步长1
e5 = extend_arr[0:nrows, 0:ncols]
e2 = extend_arr[0:nrows, 1:ncols + 1]
e6 = extend_arr[0:nrows, 2:ncols + 2]
e1 = extend_arr[1:nrows + 1, 0:ncols]
e3 = extend_arr[1:nrows + 1, 2:ncols + 2]
e8 = extend_arr[2:nrows + 2, 0:ncols]
e4 = extend_arr[2:nrows + 2, 1:ncols + 1]
e7 = extend_arr[2:nrows + 2, 2:ncols + 2]
'''
处理数据:尤其是边缘的数据
'''
e1[np.where((np.isnan(e1)) & (data != -9999))] = data[np.where((np.isnan(e1)) & (data != -9999))]
e2[np.where((np.isnan(e2)) & (data != -9999))] = data[np.where((np.isnan(e2)) & (data != -9999))]
e3[np.where((np.isnan(e3)) & (data != -9999))] = data[np.where((np.isnan(e3)) & (data != -9999))]
e4[np.where((np.isnan(e4)) & (data != -9999))] = data[np.where((np.isnan(e4)) & (data != -9999))]
e5[np.where((np.isnan(e5)) & (data != -9999))] = data[np.where((np.isnan(e5)) & (data != -9999))]
e6[np.where((np.isnan(e6)) & (data != -9999))] = data[np.where((np.isnan(e6)) & (data != -9999))]
e7[np.where((np.isnan(e7)) & (data != -9999))] = data[np.where((np.isnan(e7)) & (data != -9999))]
e8[np.where((np.isnan(e8)) & (data != -9999))] = data[np.where((np.isnan(e8)) & (data != -9999))]
def algorithm():
we = ((e8 + 2 * e1 + e5) - (e7 + 2 * e3 + e6)) / (8 * cell_size)
sn = ((e8 + 2 * e4 + e7) - (e5 + 2 * e2 + e6)) / (8 * cell_size)
s = np.around(slope(we, sn), 4) # 保留4位小数
return s
def slope(we, sn): # 坡度
num = np.square(we) + np.square(sn)
return transform(np.sqrt(num)) # 转成角度
def transform(n): # 计算反正切后转化为角度
return np.degrees(np.arctan(n))
slope_array = algorithm()
driver = gdal.GetDriverByName('GTiff')
slope_dataset = driver.Create(outname, ncols, nrows, 1, gdal.GDT_Float32)
slope_dataset.SetGeoTransform(trasnformation)
slope_dataset.SetProjection(ds.GetProjection())
slope_band = slope_dataset.GetRasterBand(1)
slope_band.WriteArray(slope_array)
slope_band.FlushCache()
# 关闭数据集
ds = None
slope_dataset = None
3. 调用函数
if __name__ == '__main__':
# 输入文件名
filename = 'your/path/name'
# 输出文件名
outname = 'your/path/name'
# 调用函数计算坡度
slope = cal_slope(filename,outname)