一、概述

山體陰影工具通過為柵格中的每個像元確定照明度，來獲取表面的假定照明度。 通過設置假定光源的位置并計算每個像元相對于相鄰像元的照明度值來實現此目的。 它可以顯著增強用于分析或圖形顯示的表面的可視化效果，尤其是在使用透明度時。

默認情況下，陰影和光是與從 0 到 255（從黑色到白色遞增）的整數關聯的灰色陰影。

二、參數

1、太陽方位角

方向角指的是太陽的角度方向，是以北為基準方向在 0 到 360 度范圍內按順時針進行測量的。 90 度的方位角指向為東。 默認方位角為 315 度 (NW)。

2、太陽高度角

高度指的是照明源高出地平線的角度或坡度。 高度的單位為度，范圍為 0（位于地平線上）到 90（位于頭上）之間。 默認值為 45 度。

3、高程縮放因子

一個表面 z 單位中地面 x,y 單位的數量。

z 單位與輸入表面的 x,y 單位不同時，可使用 z 因子調整 z 單位的測量單位。 計算最終輸出表面時，將用 z 因子乘以輸入表面的 z 值。

如果 x,y 單位和 z 單位采用相同的測量單位，則 z 因子為 1。 這是默認設置。

如果 x,y 單位和 z 單位采用不同的測量單位，則必須將 z 因子設置為適當的因子，否則會得到錯誤的結果。 例如，如果 z 單位是英尺，而 x,y 單位是米，則應使用 z 因子 0.3048 將 z 單位從英尺轉換為米（1 英尺 = 0.3048 米）。

三、示例

下面的山體陰影示例的方位角為 315 度，高度角為 45 度，z因子為1.0。

原TIF為：

四、gdal代碼實現

def calculateSlopeAndAspect(demData: Array[Float],width: Int,height: Int,demx: Double,   demy: Double, zFactor: Double): (Array[Float], Array[Float]) = {val slopeData = new Array[Float](width * height)val aspectData = new Array[Float](width * height)// 遍歷每個像素（跳過邊緣1像素，避免越界）for (y <- 0 until height - 1; x <- 0 until width - 1) {// 獲取3x3窗口的DEM值val z1 = demData((y - 1) * width + (x - 1))val z2 = demData((y - 1) * width + x)val z3 = demData((y - 1) * width + (x + 1))val z4 = demData(y * width + (x - 1))val z6 = demData(y * width + (x + 1))val z7 = demData((y + 1) * width + (x - 1))val z8 = demData((y + 1) * width + x)val z9 = demData((y + 1) * width + (x + 1))// x方向一階偏導數val dx = ((z3 + 2 * z6 + z9) - (z1 + 2 * z4 + z7)) / (8.0 * demx)// y方向一階偏導數val dy = ((z7 + 2 * z8 + z9) - (z1 + 2 * z2 + z3)) / (8.0 * demy)// 計算坡度// 公式：slope = arctan(zFactor * sqrt(dx2 + dy2))val slopeRad = atan(zFactor * sqrt(dx * dx + dy * dy))slopeData(y * width + x) = slopeRad.toFloat  // 弧度制坡度// 計算坡向var aspectRad = atan2(dy, -dx)  if (aspectRad < 0) {aspectRad += 2 * Pi  // 轉換為0~2π范圍}aspectData(y * width + x) = aspectRad.toFloat  // 弧度制坡向}(slopeData, aspectData)
}
def calculateHillShade(slopeData: Array[Float], aspectData: Array[Float], azimuth: Double, altitude: Double): Array[Float] = {//將角度值轉為弧度val azimuthRad = toRadians(azimuth)val altitudeRad = toRadians(altitude)val hillShadeData = new Array[Float](slopeData.length)for (i <- slopeData.indices) {val slope = slopeData(i)val aspect = aspectData(i)val hillShade = cos(altitudeRad) * cos(slope) + sin(altitudeRad) * sin(slope) * cos(azimuthRad - aspect)hillShadeData(i) = (255 * hillShade).toFloat}hillShadeData
}

五、geotrellis-raster代碼

def calculateHillshade(dem: Tile,cellSize: CellSize,solarAzimuth: Double,solarAltitude: Double,zFactor: Double): Tile = {val (cols, rows) = (dem.cols, dem.rows)val hillshade = ArrayTile.empty(DoubleConstantNoDataCellType, cols, rows)// 太陽參數轉換為弧度val azimuthRad = toRadians(solarAzimuth)val altitudeRad = toRadians(solarAltitude)val zenithRad = toRadians(90 - solarAltitude)def getZ(col: Int, row: Int): Double = {if (col < 0 || col >= cols || row < 0 || row >= rows) Double.NaNelse {val v = dem.getDouble(col, row)if (isNoData(v)) Double.NaN else v}}// 遍歷每個像素for (row <- 0 until rows; col <- 0 until cols) {val (z1, z2, z3) = (getZ(col-1, row-1), getZ(col, row-1), getZ(col+1, row-1))val (z4, z5, z6) = (getZ(col-1, row),   getZ(col, row),   getZ(col+1, row))val (z7, z8, z9) = (getZ(col-1, row+1), getZ(col, row+1), getZ(col+1, row+1))if (!e.isNaN && isData(e)) {// 計算x、y方向導數val dzdx = ((z3 + 2*z6 + z9) - (z1 + 2*z4 + z7)) / (8 * cellSize.width)val dzdy = ((z7 + 2*z8 + z9) - (z1 + 2*z2 + z3)) / (8 * cellSize.height)// 計算坡度（弧度）val slopeRad = atan(zFactor * math.sqrt(dzdx*dzdx + dzdy*dzdy))// 計算坡向（度）val aspectDeg = {if (dzdx == 0 && dzdy == 0) 0.0else {var aspect = math.toDegrees(atan2(dzdy, -dzdx))if (aspect < 0) aspect + 360 else aspect}}// 計算山體陰影val cosSlope = cos(slopeRad)val sinSlope = sin(slopeRad)val cosZenith = cos(zenithRad)val sinZenith = sin(zenithRad)val deltaAzimuth = azimuthRad - toRadians(aspectDeg)val cosDelta = cos(deltaAzimuth)val shade = math.max(0.0, cosZenith * cosSlope + sinZenith * sinSlope * cosDelta)hillshade.setDouble(col, row, (shade * 255.0))} else {hillshade.setDouble(col, row, Double.NaN)}}hillshade
}

六、3x3網格示例