随笔- 183 文章- 0 评论- 18 阅读- 41万

【GDAL】python读取DEM计算坡度与坡向

利用GDAL读入DEM与Landsat影像，由于DEM是WG84坐标系，Landsat是WGS84坐标系UTM投影，因此处理在实际应用中需要将DEM进行投影转换。

大概分为以下几个步骤：

读取DEM，读取Landsat影像
获取Landsat影像的投影信息，将其赋给DEM，并对DEM进行重采样
计算dx和dy
计算坡度和坡向
输出坡度和坡向的影像

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