python 计算坡度

根据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)