GDAL使用DEM資料計算山體陰影（Hillshade）

零、        前言

說起Hillshade，其實就是類比太陽光照射地形所引起的明暗對比，然後來對地形圖進行渲染，使之看起來具有立體效果的一種方式，常用於地圖的渲染，如表1所示，具體的可以參考文獻[1]，表1中的圖均來自參考文獻[1]。

表1 DEM、山體陰影以及應用對比

DEM映像（使用顏色渲染）	從左圖的DEM映像中計算的山體陰影圖
Paper Map Without Hillshade	Paper Map With Hillshade
Satellite Image Without Hillshade	Satellite Image With Hillshade

從表1中的圖可以看出，使用山體陰影對地圖和衛星映像進行渲染後，視覺效果得到了很大的提升，具有很明顯的立體感。下面就對如何?山體陰影進行一個說明。

一、        山體陰影原理簡介

山體陰影通過為柵格中的每個像元指定太陽高度等資訊，來計算表面的假定亮度值。通過設定假定光源的位置和計算與相鄰像元相關的每個像元的亮值，即可得出假定亮度度。進行分析或圖形顯示時，特別是使用透明度後，“山體陰影”可大大增強表面的可視化。

預設情況下，陰影和光線是與介於 0 和 255 之間的整數相關的灰階梯度（從黑色漸層為白色）。【摘自參考文獻[4,5]】。

山體陰影計算的時候，需要指定太陽在天空中的位置，這個位置可以用下面兩個參數來進行描述。

1、太陽方位角

太陽方向角指的是太陽的角度方向，是以北為基準方向，在0度到360度範圍內按順時針進行測量的。90º的方位角為東。預設方向角為 315º (NW，西北方向)。1所示，圖中黃色圓圈代表太陽。

2、太陽高度角

太陽高度角指的是太陽高出地平線的角度或坡度。高度的單位為度，範圍為0度（位於地平線上）到 90度（位於頭頂上）之間。預設值為45度。2所示，途中黃色圓圈代表太陽。

               

圖1 太陽方位角                圖2 太陽高度角

二、        計算方法

在計算山體陰影的公式如式（1）所示，詳情參考文獻[6]。式（1）中的Zenith表示太陽天頂角，Azimuth表示太陽方位角，Slope和Aspect分別表示坡度和坡向。尾碼_rad表示所有的角度都是以弧度為單位的。

 Hillshade = 255.0 * ((cos(Zenith_rad) * cos(Slope_rad)) + (sin(Zenith_rad) * sin(Slope_rad) * cos(Azimuth_rad - Aspect_rad)))（1）

通過上式計算山體陰影時，計算的結果可能是小於0的值，此時應該將該值設定為0。下面對計算山體陰影的過程進行拆分說明，主要有下面幾個步驟：

1、計算太陽入射角度

指定太陽的高度角必須大於地平線的角度（也就是0度）。但是，用於計算山體陰影值的公式要求以弧度為單位表示角度並且要求是從垂直方向偏轉的角度。將垂直於地表面的方向（頭頂正上方）標註為天頂。天頂角是從天頂點到太陽的方向測之間的角度，也就是太陽高度角的餘角（即90度減去太陽高度角）。要計算太陽入射角度，第一步是將太陽高度角轉換為天頂角。第二步是將天頂角轉換為弧度。

將太陽高度角轉換為天頂角：

Zenith_deg = 90 - Altitude   （2）

轉換為弧度：

Zenith_rad = Zenith * pi / 180.0  （3）

2、計算太陽方位角

山體陰影公式要求方位角採用弧度作為單位。首先，將天頂角從地理單位（羅盤方向）轉換為數學單位（直角）。然後再將方位角轉換為弧度。具體見公式（4）（5）（6）。

轉換太陽方位角的方法：

Azimuth_math = 360.0 - Azimuth + 90  (4)

請注意，如果，Azimuth_math >= 360.0 則：

Azimuth_math = Azimuth_math - 360.0  (5)

轉換為弧度：

   Azimuth_rad = Azimuth_math * pi / 180.0   (6)

3、計算坡度坡向

關於坡度坡向的計算方法這裡就不說了，可以參考前面兩篇部落格，具體參考這裡《GDAL使用DEM資料計算坡度坡向》，參考文獻[7]。需要注意的是，要將計算的坡度和坡向轉為弧度單位。

至此，所有的變數都計算出來了，然後帶入公式（1）即可求得山體陰影的值。圖3是一個山體陰影計算的，太陽方位角為99度，太陽高度角為33度的計算結果。

圖3 山體陰影計算樣本

三、        演算法編寫

由於計算山體陰影還是一個3×3的視窗，所以我們繼續用之前的那個3×3的演算法作為基礎，這裡只實現演算法部分即可。在計算的過程中需要考慮到像元的大小和高程的單位，如果兩者的單位不一致，需要做一個轉換，此外還有的時候需要對高程進行一個誇大，目的是為了更突出地形資訊。

首先是山體陰影的演算法結構體，在演算法結構體中對輸入的參數提前進行了一系列的處理，比如角度直接計算出正弦和餘弦值，這樣就會在後面的演算法中直接使用，而不用再次計算，提高演算法的效率，當然都是些小的細節問題。

/*** @brief 山體陰影演算法結構體*/typedef struct{         /*! 南北方向解析度 */         double nsres;         /*! 東西方向解析度 */         double ewres;         /*! 方位角的正弦值 */         doublesin_altRadians;         /*! 方位角的餘弦值乘以Z縮放比 */         doublecos_altRadians_mul_z_scale_factor;         /*! 太陽高度角 */         doubleazRadians;         /*! z縮放比平方 */         doublesquare_z_scale_factor;} GDALHillshadeAlgData;

接下來是演算法實現函數，這個函數基本就是按照上面的公式來寫的，沒啥難度。

float GDALHillshadeAlg (float* afWin, floatfDstNoDataValue, void* pData){         GDALHillshadeAlgData*psData = (GDALHillshadeAlgData*)pData;         double x, y,aspect, xx_plus_yy, cang;          // 計算坡度         x =((afWin[0] + afWin[3] + afWin[3] + afWin[6]) -                   (afWin[2]+ afWin[5] + afWin[5] + afWin[8])) / psData->ewres;          y =((afWin[6] + afWin[7] + afWin[7] + afWin[8]) -                   (afWin[0]+ afWin[1] + afWin[1] + afWin[2])) / psData->nsres;          xx_plus_yy =x * x + y * y;          // 計算坡向         aspect =atan2(y, x);          // 計算山體陰影值         cang =(psData->sin_altRadians -                   psData->cos_altRadians_mul_z_scale_factor* sqrt(xx_plus_yy) *                   sin(aspect- psData->azRadians)) /                   sqrt(1+ psData->square_z_scale_factor * xx_plus_yy);          if (cang<= 0.0)                   cang= 1.0;         else                   cang= 1.0 + (254.0 * cang);          return(float)cang;}

最後還需要一個建立演算法所需參數的一個函數，這個函數輸入的有映像的六參數，用於擷取映像的解析度；z是高程的縮放比例，就是有的時候地形比較平坦，出來的效果不太好，可以使用z值來進行一個誇大；sacle是水平像素單位和高程單位的一個換算係數，有的時候像素單位是按照經緯度的，而高程單位一般是米，這時就需要設定sacle了；最後兩個參數就是上面說的方位角和高度角。

void* GDALCreateHillshadeData(double* adfGeoTransform,                                                                    double z,                                                                    double scale,                                                                    double alt,                                                                    double az){         GDALHillshadeAlgData*pData =                   (GDALHillshadeAlgData*)CPLMalloc(sizeof(GDALHillshadeAlgData));          const doubledegreesToRadians = M_PI / 180.0;         pData->nsres= adfGeoTransform[5];         pData->ewres= adfGeoTransform[1];         pData->sin_altRadians= sin(alt * degreesToRadians);         pData->azRadians= az * degreesToRadians;         doublez_scale_factor = z / (8 * scale);         pData->cos_altRadians_mul_z_scale_factor=                   cos(alt* degreesToRadians) * z_scale_factor;         pData->square_z_scale_factor= z_scale_factor * z_scale_factor;         return pData;}

最後在網上有一個使用Matlab編寫的計算山體陰影的m檔案，具體參考文獻[2]。有興趣的可以下載下來用Matlab試試。

四、        處理結果

最後沒啥好說的，貼張處理後的圖給大家看看（方位角為315度，高度角為45度）。希望對大家有用。

圖4 DEM資料

圖5 DEM渲染後的資料

圖5 高程縮放為1時的效果

圖5 高程縮放為2時的效果

圖4 高程縮放為10時的效果

五、        參考文獻

[1] http://www.gichd.org/fileadmin/pdf/IMSMA/fact-sheets/GICHDFactSheet-Hillshade%20Imagery.PDF(Hillshade的一個說明)

[2] http://www.mathworks.ch/matlabcentral/fileexchange/14863-hillshade/content/hillshade.m(MatLab的版本)

[3]http://download.osgeo.org/ossim/tutorials/pdfs/HillShade.pdf（OSSIM軟體中的使用說明）

[4]http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/na/009z000000v0000000/

[5]http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/How_Hillshade_works/009z000000z2000000/

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

[7] http://blog.csdn.net/liminlu0314/article/details/8498985（GDAL計算坡度坡向）

[8] http://www.aubreyrhea.net/gis/index.php/tag/hillshade/（ArcMap軟體中的使用說明）