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