【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;
}