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

 

 

 

 

posted @   羊大葱  阅读(2454)  评论(2编辑  收藏  举报
相关博文:
阅读排行:
· Manus重磅发布:全球首款通用AI代理技术深度解析与实战指南
· 被坑几百块钱后,我竟然真的恢复了删除的微信聊天记录!
· 没有Manus邀请码?试试免邀请码的MGX或者开源的OpenManus吧
· 园子的第一款AI主题卫衣上架——"HELLO! HOW CAN I ASSIST YOU TODAY
· 【自荐】一款简洁、开源的在线白板工具 Drawnix
点击右上角即可分享
微信分享提示