python 计算坡度

本文介绍了如何使用Python的Numpy库在地理信息系统中进行批量坡度计算,详细解释了通过九宫格法计算坡度的过程,以及如何处理边缘数据和生成坡度数组。
摘要由CSDN通过智能技术生成

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

  • 7
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
坡度计算是根据地形高程数据计算某一区域的坡度变化情况。在使用Python进行坡度计算时,可以利用数字高程模型(Digital Elevation Model,DEM)数据。DEM数据是描述地形高程的数字化表示,可以使用DEM数据来进行坡度计算。 一种常用的计算坡度的方法是使用坡度公式:坡度 = arctan((dz/dx)^2 + (dz/dy)^2)的平方根。其中dz是垂直距离的变化,dx和dy是水平距离的变化。 在Python中可以使用一些库来进行坡度计算,例如使用gdal库来读取DEM数据,然后根据坡度公式计算每个像素点的坡度。以下是一个简单的示例代码: ```python import gdal import numpy as np import math # 读取DEM数据 dataset = gdal.Open('dem.tif') band = dataset.GetRasterBand(1) dem = band.ReadAsArray() # 计算DEM数据的形状 rows, cols = dem.shape # 定义坡度数组 slope = np.zeros((rows, cols)) # 计算坡度 for row in range(1, rows-1): for col in range(1, cols-1): dz_dx = (dem[row+1, col] - dem[row-1, col]) / 2 dz_dy = (dem[row, col+1] - dem[row, col-1]) / 2 slope[row, col] = math.sqrt((dz_dx ** 2) + (dz_dy ** 2)) # 保存坡度数据 driver = gdal.GetDriverByName('GTiff') out_dataset = driver.Create('slope.tif', cols, rows, 1, gdal.GDT_Float32) out_dataset.GetRasterBand(1).WriteArray(slope) # 清除缓存并关闭数据集 band.FlushCache() dataset = None out_dataset = None ``` 上述代码中,首先使用gdal库读取DEM数据,然后根据坡度公式计算每个像素点的坡度值,并将坡度数据保存为一个新的.tif文件。 以上是一个简单的Python代码示例,用于计算坡度。根据实际需求,可能需要进行更复杂的处理,例如计算平均坡度或者根据阈值划分斜坡等。因此,根据具体任务需要,可以对代码进行进一步的修改和扩展。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值