如何利用Python对DEM数据进行显示（对之前代码的更正），附带利用Python计算坡度

1、更正---图像显示代码

在之前的某一篇文章中，《如何利用Python对遥感影像进行显示》，阐述了如何利用Python对遥感影像进行显示。最近利用代码对DEM数据进行显示时，输入DEM数据，显示的结果却是全黑的，仔细看了一下之前的代码，发现了一个小的BUG，如下：

这句代码有点小问题，当 (im_max-im_min) 值比较大，大于255时，255/(im_max-im_min) 的结果就默认为0，因为默认为int型。如此，后面和任何数字相乘都为0，导致整个显示图像都为0。因此，该部分代码应更正为：

def img_processing(im_band,img_data):
    if im_band == 1:
        data_jpg = np.zeros((img_data.shape[0],img_data.shape[1]),dtype='uint8')
        im_max = np.amax(img_data)
        im_min = np.amin(img_data)
        for m in range(0, img_data.shape[0]):
            for n in range(0, img_data.shape[1]):
                data_jpg[m,n] = float(255./(im_max-im_min))*(img_data[m,n]-im_min)  ## 需要注意

    else:
        data_jpg = np.zeros((img_data.shape[1],img_data.shape[2],3),dtype='uint8')
        for i in range(3):
            im_max = np.amax(img_data[i,:,:])
            im_min = np.amin(img_data[i,:,:])
            for m in range(0, img_data.shape[1]):
                for n in range(0, img_data.shape[2]):
                    data_jpg[m,n,i] = float(255./(im_max-im_min))*(img_data[i,m,n]-im_min)  ## 需要注意
    return data_jpg

具体就是，将 255/(im_max-im_min) 计算结果设为浮点型，这样后面就没问题啦。

利用更正后的代码，对DEM影像显示结果如下：

2、利用Python计算坡度

关于坡度如何计算，网上应该有很多文章和公式可以查阅，这里不详细介绍。

直接给出Python计算坡度的代码，输入数据是DEM数据，输出数据是计算后的坡度数据，包括弧度和角度两种形式。

注：DEM数据来源于ASTER DEM数据，空间分辨率为30米。

# -*- coding: utf-8 -*-
import os
import glob
import math
import cv2
import numpy as np
from osgeo import ogr,osr,gdal


### GDAL影像读取函数 ###
class IMAGE:

    #读图像文件
    def read_img(self,filename):
        dataset = gdal.Open(filename)    #打开文件
        if dataset == "None":
            QtWidgets.QMessageBox.critical(None, "错误", "影像读取错误！")
        im_width = dataset.RasterXSize  #栅格矩阵的列数
        im_height = dataset.RasterYSize  #栅格矩阵的行数

        im_geotrans = dataset.GetGeoTransform() #仿射矩阵
        im_proj = dataset.GetProjection() #地图投影信息
        im_data = dataset.ReadAsArray(0,0,im_width,im_height) #将数据写成数组，对应栅格矩阵
        del dataset 
        return im_proj,im_geotrans,im_data

    #写文件，以写成tiff为例
    def write_img(self,filename,im_proj,im_geotrans,im_data):
        
        #判断栅格数据的数据类型
        if 'int8' in im_data.dtype.name:
            datatype = gdal.GDT_Byte
        elif 'int16' in im_data.dtype.name:
            datatype = gdal.GDT_UInt16
        else:
            datatype = gdal.GDT_Float32

        #判读数组维数
        if len(im_data.shape) == 3:
            im_bands, im_height, im_width = im_data.shape
        else:
            im_bands, (im_height, im_width) = 1,im_data.shape 

        #创建文件
        driver = gdal.GetDriverByName("GTiff")      #数据类型必须有，因为要计算需要多大内存空间
        dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)

        dataset.SetGeoTransform(im_geotrans)       #写入仿射变换参数
        dataset.SetProjection(im_proj)          #写入投影

        if im_bands == 1:
            dataset.GetRasterBand(1).WriteArray(im_data) #写入数组数据
        else:
            for i in range(im_bands):
                dataset.GetRasterBand(i+1).WriteArray(im_data[i])

        del dataset

### 计算坡度 ###
def calculate_slope(filename, save_path):
    image_grid = IMAGE()
    im_proj, im_geotrans, im_data = image_grid.read_img(filename)  ## 读DEM
    image_expand = np.pad(im_data, 1, 'edge')  ## DEM扩充, 避免边缘值损失
    size_x = abs(im_geotrans[1])  ## 绝对值
    size_y = abs(im_geotrans[-1])  ## 绝对值

    dx = ((image_expand[1:-1,:-2])-(image_expand[1:-1,2:]))/(size_x*2)
    dy = ((image_expand[2:,1:-1])-(image_expand[:-2,1:-1]))/(size_y*2)


    slope_radian = np.arctan(np.sqrt(dx**2 + dy**2))  ## 求反正切值
    slope_degree = slope_radian*180/(math.pi)  ## 转换成度

    image_grid.write_img(save_path+"/slope_radian.tif",im_proj,im_geotrans,slope_radian)
    image_grid.write_img(save_path+"/slope_degree.tif",im_proj,im_geotrans,slope_degree)

if __name__ == '__main__':
    calculate_slope('DEM.tif', 'results')

计算结果如下

可以看到，代码计算的坡度和利用ArcGIS Toolbox计算的坡度，目视上差不多，但数值上还稍微有些差距，代码得到的最大坡度要比ArcGIS计算的结果稍稍偏大，可能计算方式存在差别吧！

 

欢迎大家对代码进行评测，有问题可以在评论中反馈，转载请注明来源！