【GDAL】python讀取DEM計算坡度與坡向

利用GDAL讀入DEM與Landsat影像，由於DEM是WG84坐標系，Landsat是WGS84坐標系UTM投影，因此處理在實際應用中需要將DEM進行投影轉換。

大概分為以下幾個步驟：

1. 讀取DEM，讀取Landsat影像
2. 獲取Landsat影像的投影信息，將其賦給DEM，並對DEM進行重采樣
3. 計算dx和dy
4. 計算坡度和坡向
5. 輸出坡度和坡向的影像

from osgeo import gdal
import os
import sys
import osr
import math
import numpy as np
from matplotlib import pyplot as plt
from gdalconst import *
 
 
# 給柵格最外圈加一圈
def assignBCs(elevGrid):
    ny, nx = elevGrid.shape
    Zbc = np.zeros((ny + 2, nx + 2))
    Zbc[1:-1, 1:-1] = elevGrid
 
    Zbc[0, 1:-1] = elevGrid[0, :]
    Zbc[-1, 1:-1] = elevGrid[-1, :]
    Zbc[1:-1, 0] = elevGrid[:, 0]
    Zbc[1:-1, -1] = elevGrid[:, -1]
 
    Zbc[0, 0] = elevGrid[0, 0]
    Zbc[0, -1] = elevGrid[0, -1]
    Zbc[-1, 0] = elevGrid[-1, 0]
    Zbc[-1, -1] = elevGrid[-1, 0]
 
    return Zbc
 
 
# 計算dx,dy
def calcFiniteSlopes(elevGrid, dx):
    Zbc = assignBCs(elevGrid)
 
    Sx = (Zbc[1:-1, :-2] - Zbc[1:-1, 2:]) / (2 * dx)  # WE方向
    Sy = (Zbc[2:, 1:-1] - Zbc[:-2, 1:-1]) / (2 * dx)  # NS方向
 
    return Sx, Sy
 
 
# 投影轉換
def convertProjection(data, filename):
    landsatData = gdal.Open(filename, GA_ReadOnly)
 
    oldRef = osr.SpatialReference()
    oldRef.ImportFromWkt(data.GetProjectionRef())
 
    newRef = osr.SpatialReference()
    newRef.ImportFromWkt(landsatData.GetProjectionRef())
 
    transform = osr.CoordinateTransformation(oldRef, newRef)
 
    tVect = data.GetGeoTransform()
    nx, ny = data.RasterXSize, data.RasterYSize
    (ulx, uly, ulz) = transform.TransformPoint(tVect[0], tVect[3])
    (lrx, lry, lrz) = transform.TransformPoint(tVect[0] + tVect[1] * nx, tVect[3] + tVect[5] * ny)
 
    memDrv = gdal.GetDriverByName('MEM')
 
    dataOut = memDrv.Create('name', int((lrx - ulx) / dx), int((uly - lry) / dx), 1, gdal.GDT_Float32)
 
    newtVect = (ulx, dx, tVect[2], uly, tVect[4], -dx)
 
    dataOut.SetGeoTransform(newtVect)
    dataOut.SetProjection(newRef.ExportToWkt())
 
    # Perform the projection/resampling
    res = gdal.ReprojectImage(data, dataOut, oldRef.ExportToWkt(), newRef.ExportToWkt(), gdal.GRA_Cubic)
 
    return dataOut
 
 
if __name__ == '__main__':
    DEMFilename = 'E:/LandsatDEM/dem/DEM.tif'
    LandsatFilename = 'E:/LandsatDEM/clip/L7-B1.tif'
    slopeFilename = 'E:/LandsatDEM/result/slope_prj.tif'
    aspectFilename = 'E:/LandsatDEM/result/aspect_prj.tif'
 
    gdal.AllRegister()
 
    data = gdal.Open(DEMFilename, GA_ReadOnly)
    if data is None:
        print('Cannot open this file:' + DEMFilename)
        sys.exit(1)
 
    dx = 30  # 分辨率
 
    # 投影變換
    projData = convertProjection(data, LandsatFilename)
 
    gridNew = projData.ReadAsArray().astype(np.float)
 
    Sx, Sy = calcFiniteSlopes(gridNew, dx)
    # 坡度計算
    slope = np.arctan(np.sqrt(Sx ** 2 + Sy ** 2)) * 57.29578
 
    # 坡向計算
    aspect = np.ones([Sx.shape[0], Sx.shape[1]]).astype(np.float32)
    for i in range(Sx.shape[0]):
        for j in range(Sy.shape[1]):
            sx = float(Sx[i, j])
            sy = float(Sy[i, j])
            if (sx == 0.0) & (sy == 0.0):
                aspect[i, j] = -1
            elif sx == 0.0:
                if sy > 0.0:
                    aspect[i, j] = 0.0
                else:
                    aspect[i, j] = 180.0
            elif sy == 0.0:
                if sx > 0.0:
                    aspect[i, j] = 90.0
                else:
                    aspect[i, j] = 270.0
            else:
                aspect[i, j] = float(math.atan2(sy, sx) * 57.29578)
                if aspect[i, j] < 0.0:
                    aspect[i, j] = 90.0 - aspect[i, j]
                elif aspect[i, j] > 90.0:
                    aspect[i, j] = 360.0 - aspect[i, j] + 90.0
                else:
                    aspect[i, j] = 90.0 - aspect[i, j]
 
    # 輸出坡度坡向文件
    driver = gdal.GetDriverByName('GTiff')
    if os.path.exists(slopeFilename):
        os.remove(slopeFilename)
    if os.path.exists(aspectFilename):
        os.remove(aspectFilename)
 
    ds1 = driver.Create(slopeFilename, slope.shape[1], slope.shape[0], 1, GDT_Float32)
    band = ds1.GetRasterBand(1)
    band.WriteArray(slope, 0, 0)
 
    ds2 = driver.Create(aspectFilename, aspect.shape[1], aspect.shape[0], 1, GDT_Float32)
    band = ds2.GetRasterBand(1)
    band.WriteArray(aspect, 0, 0)
 
    del ds1
    del ds2
    data = None
    projData = None

 

 

 

 

 

 

 

 

/**
* @brief 坡度算法數據結構體
*/
typedef struct
{
         /*! 南北方向分辨率 */
         double nsres;
         /*! 東西方向分辨率 */
         double ewres;
         /*! 縮放比例 */
         double scale;
         /*! 坡度方式 */
         int    slopeFormat;
} GDALSlopeAlgData;

接下來是坡度算法的回調函數，該函數就是計算坡度的函數，按照上面的公式寫就好，最后需要根據指定的坡度表達方式進行轉換即可：

float GDALSlopeAlg (float* afWin, float fDstNoDataValue,void* pData)
{
         const doubleradiansToDegrees = 180.0 / M_PI;
         GDALSlopeAlgData*psData = (GDALSlopeAlgData*)pData;
         double dx,dy, key;
 
         dx =((afWin[0] + afWin[3] + afWin[3] + afWin[6]) -
                   (afWin[2]+ afWin[5] + afWin[5] + afWin[8])) / psData->ewres;
 
         dy =((afWin[6] + afWin[7] + afWin[7] + afWin[8]) -
                   (afWin[0]+ afWin[1] + afWin[1] + afWin[2])) / psData->nsres;
 
         key = (dx *dx + dy * dy);
 
         if(psData->slopeFormat == 1)
                   return(float)(atan(sqrt(key) / (8*psData->scale)) * radiansToDegrees);
         else
                   return(float)(100*(sqrt(key) / (8*psData->scale)));
}

最后是創建坡度結構體的函數：

void* GDALCreateSlopeData(double* adfGeoTransform,       double scale, int slopeFormat)
{
         GDALSlopeAlgData*pData =
                   (GDALSlopeAlgData*)CPLMalloc(sizeof(GDALSlopeAlgData));
 
         pData->nsres= adfGeoTransform[5];
         pData->ewres= adfGeoTransform[1];
         pData->scale= scale;
         pData->slopeFormat= slopeFormat;
         return pData;
}

 

2、坡向

首先是定義坡向算法的結構體，坡向的參數很簡單，就是一個是否使用方位角來表示，意思就是如果這個值設置為TRUE，坡度的計算結果按照圖2中的角度進行表示，如果是FALSE，計算的結果是多少就是多少：

/**
* @brief 坡向算法數據結構體
*/
typedef struct
{
         /*! 方位角 */
         intbAngleAsAzimuth;
} GDALAspectAlgData;

接下來是坡向算法的回調函數，按照上面的公式寫的，沒啥難度。需要注意的是如果指定了bAngleAsAzimuth是TRUE的話，需要把計算的角度做一個轉換，轉換的結果就是0表示東，90表示北，180表示西，270表示南：

float GDALAspectAlg (float* afWin, float fDstNoDataValue,void* pData)
{
         const doubledegreesToRadians = M_PI / 180.0;
         GDALAspectAlgData*psData = (GDALAspectAlgData*)pData;
         double dx,dy;
         float aspect;
 
         dx =((afWin[2] + afWin[5] + afWin[5] + afWin[8]) -
                   (afWin[0]+ afWin[3] + afWin[3] + afWin[6]));
 
         dy =((afWin[6] + afWin[7] + afWin[7] + afWin[8]) -
                   (afWin[0]+ afWin[1] + afWin[1] + afWin[2]));
 
         aspect =(float)(atan2(dy, -dx) / degreesToRadians);
 
         if (dx == 0&& dy == 0)
         {
                   /*Flat area */
                   aspect= fDstNoDataValue;//或者這里直接是-1
         }
         else if (psData->bAngleAsAzimuth )
         {
                   if(aspect > 90.0)
                            aspect= 450.0f - aspect;
                   else
                            aspect= 90.0f - aspect;
         }
         else
         {
                   if(aspect < 0)
                            aspect+= 360.0;
         }
 
         if (aspect ==360.0)
                   aspect= 0.0;
 
         returnaspect;
}

最后是創建坡向結構體的函數：

void* GDALCreateAspectData(int bAngleAsAzimuth)
{
         GDALAspectAlgData*pData =
                   (GDALAspectAlgData*)CPLMalloc(sizeof(GDALAspectAlgData));
 
         pData->bAngleAsAzimuth= bAngleAsAzimuth;
         return pData;
}