GDAL使用DEM資料計算坡度坡向

零、        前言

之前寫過一個3×3的通用模板運算元函數的部落格《基於GDAL的一個通用的3×3模板函數》，網址：http://blog.csdn.net/liminlu0314/article/details/8316156。當時說是要基於這個函數寫一個計算坡度坡向的函數。由於這段時間一直忙於別的事情，這件事情就拖著了，今天給大家補上。

一、        簡介

坡度（slope）是地表單元陡緩的程度，通常把坡面的垂直高度h和水平距離l的比叫做坡度（或叫做坡比）用字母i表示。【即坡角的正切值（可寫作：i=tan坡角）】。坡度的表示方法有百分比法、度數法、密位法和分數法四種，其中以百分比法和度數法較為常用（1所示）。

圖1 坡度的兩種表示方法

1、百分比法

表示坡度最為常用的方法，即兩點的高程差與其水平距離的百分比，其計算公式如下：

坡度 = (高程差/水平距離)x100%

使用百分比表示時，即：i=h/l×100%

例如：坡度3% 是指水平距離每100米,垂直方向上升(下降)3米；1%是指水平距離每100米，垂直方向上升(下降)1米。以次類推！

2、度數法

用度數來表示坡度，利用反三角Function Compute而得，其公式如下：

                                                                tanα(坡度)= 高程差/水平距離

                                                                所以α(坡度)= arctan(高程差/水平距離)

如果將高程增量百分比視為高程增量除以水平增量後再乘以 100，就可以更好地理解高程增量百分比。請考慮下面的三角形 B。當角度為 45 度時，高程增量等於水平增量，所以高程增量百分比為 100%。如三角形 C 所示，當坡度角接近直角（90 度）時，高程增量百分比開始接近無窮大。

坡向(aspect) 是指地形坡面的朝向。坡向用於識別出從每個像元到其相鄰像元方向上值的變動率最大的下坡方向。坡向可以被視為坡度方向。坡向是一個角度，將按照順時針方向進行測量，角度範圍介於 0（正東）到 360（仍是正東）之間，即完整的圓。不具有下坡方向的平坦地區將賦值為-1。2所示。

圖2 坡向的取值

二、        坡度坡向計算方法

坡度計算一般採用擬合曲面法。擬合曲面一般採用二次曲面，即3×3的視窗，3所示。每個視窗的中心為一個高程點。圖3中中心點e的坡度和坡向的計算公式如下：

上式中：Slope為坡度，Aspect為坡向，Slopewe為X方向的坡度，Slopesn為Y方向的坡度。

圖3 一個3×3的視窗

關於Slopewe、Slopesn的計算可以採用一下幾種常用的方法：

l  演算法1

l  演算法2

l  演算法3

   

l  演算法4

   

上式中的Cellsize為格網DEM的間隔長度。演算法1的精度最高，計算效率也最高，其次是演算法2。ERDAS Imagine中採用的是演算法4，ArcMap採用的是演算法2。本次也採用演算法2進行實現。

三、        坡度坡向演算法的編寫

這裡編寫計算坡度和坡向的函數依賴於之前的部落格《基於GDAL的一個通用的3×3模板函數》中的函數，所以這裡唯寫出計算坡度和坡向的演算法代碼，具體調用參考部落格《基於GDAL的一個通用的3×3模板函數》。通過部落格《基於GDAL的一個通用的3×3模板函數》我們知道，要實現一個演算法，需要寫三部分的內容，分別是：演算法參數的結構體；一個演算法實現的回呼函數和建立演算法參數結構體指標的函數。

1、坡度

首先是定義坡度演算法的結構體，計算坡度需要的參數有DEM格網的大小，也就是東西方向和南北方向的解析度；高程的縮放比例；以及坡度的表示方式，即百分比還是度：

/*** @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;}

四、        結果展示

原始的DEM資料4所示，圖5是計算的坡度值，圖6是計算的坡向值。以下三個映像均是展開後的顯示效果，由於DEM資料是16bit資料，計算的坡度和坡向資料是一個32位浮點型的資料。映像也沒有進行分級渲染，用ArcMap處理的坡度和坡向資料做了一個彩色的渲染，效果比較好，當然這個就不再這裡討論了。

圖4 澳洲新南威爾士州的一塊DEM（SRTM資料）

圖5 計算的坡度結果

圖6 計算的坡向結果

五、        參考文獻：

[1]Xiangling LIU, Glenn Otto (July 2013). The Determinants of Supply Elasticity in Sydney Housing Market,Working Paper,University of New South Wales

[2] http://help.arcgis.com/zh-cn/arcgisdesktop/10.0/help/index.html#//009z000000v2000000

[3] http://baike.baidu.com/view/1353583.htm

[4] http://help.arcgis.com/zh-cn/arcgisdesktop/10.0/help/index.html#//009z000000tr000000

[5] Burrough, P. A., and McDonell, R. A.,1998. Principles of Geographical Information Systems (OxfordUniversity Press, New York), 190 pp.

[6]http://baike.baidu.com/view/762866.htm