Python计算DEM(tif格式)坡度和坡向

最新推荐文章于 2024-05-21 08:00:00 发布

见贤思齐547

最新推荐文章于 2024-05-21 08:00:00 发布

阅读量7.2k

点赞数 15

分类专栏： Python地理数据处理文章标签： python

本文链接：https://blog.csdn.net/weixin_40501429/article/details/114894497

版权

slopeaspect将背景值设为-300slopeArcgis

摘要由CSDN通过智能技术生成

本文根据山东省DEM图获取坡度、坡向图，使用了三种方式：Python GDAL工具自带的函数处理、Python中自己编写函数实现和arcgis中实现。

一.Python中实现（针对TIF格式的DEM数据）

1.利用gdal工具处理

（1）代码

from osgeo import gdal, osr

# 获取影像信息
infoDEM = gdal.Info(r"D:\ProfessionalProfile\DEMdata\2_OutputArcMap\demsd0330UTM.tif")


# 计算坡度、坡向
slope = gdal.DEMProcessing(r'D:\ProfessionalProfile\DEMdata\6_DEMXuanRan\slopePy.tif',
                           r"D:\ProfessionalProfile\DEMdata\2_OutputArcMap\demsd0330UTM.tif", "slope")
aspect = gdal.DEMProcessing(r'D:\ProfessionalProfile\DEMdata\6_DEMXuanRan\aspectPy.tif',
                            r"D:\ProfessionalProfile\DEMdata\2_OutputArcMap\demsd0330UTM.tif", "aspect", format='GTiff',
                            trigonometric=0, zeroForFlat=1)

（2）结果

原图

slope图

aspect图

2.参考了网上的部分资料，自己写了一下利用DEM.tif格式的DEM数据获取坡度坡向。

（1）代码

大概分为几个步骤：读取DEM影像——计算梯度——计算坡度坡向——输出TIF影像，每一步都对应着相应的函数。

from osgeo import gdal,ogr,osr
import numpy as np
import math
import datetime


# 读取TIFF遥感影像
def read_img(filename):

    dataset = gdal.Open(filename)  # 打开文件
    # dataset = gdal.Open(r'D:\ProfessionalProfile\DEMdata\slopeAspectPython0322\test.tif')
    im_width = dataset.RasterXSize  # 栅格矩阵的列数
    im_height = dataset.RasterYSize  # 栅格矩阵的行数
    im_bands = dataset.RasterCount  # 波段数
    im_geotrans = dataset.GetGeoTransform()  # 仿射矩阵，左上角像素的大地坐标和像素分辨率
    im_proj = dataset.GetProjection()  # 地图投影信息，字符串表示
    im_data = dataset.ReadAsArray(0, 0, im_width, im_height)
    datatype = im_data.dtype
    del dataset  # 关闭对象dataset，释放内存

    return   im_data, im_proj, im_geotrans, im_width, im_height, im_bands, datatype


# 为便于后续坡度计算，需要在原图像的周围添加一圈数值
def AddRound(npgrid):

    nx, ny = npgrid.shape[0], npgrid.shape[1]   # ny:行数&#