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計算的結果稍稍偏大，可能計算方式存在差別吧！

歡迎大家對代碼進行評測，有問題可以在評論中反饋，轉載請注明來源！