DEM山体阴影原理以及算法详解

山体阴影原理以及算法详解

山体阴影基本原理：

山体阴影是假想一个光源在某个方向和某个太阳高度的模拟下，用过临近像元的计算来生成一副0-255的灰度图。

一、山体阴影的主要参数：

1、  太阳光线的入射角度：这个角度的量算起点是正北方向，按照顺时针的方向，角度的范围是0到360度，如下图所示，默认的角度是315度，西北方向，如下图所示：

2、  太阳高度角：太阳高度角也简称太阳高度。是太阳光线和当地地平面之间的夹角，范围是0-90度，默认的太阳高度是45度，如下图所示：

二、山体阴影计算方法

山体阴影的计算公式如下

(1)  Hillshade = 255.0 * ((cos(Zenith_rad)* cos(Slope_rad)) +

(sin(Zenith_rad) * sin(Slope_rad) * cos(Azimuth_rad- Aspect_rad)))

如果Hillshade < 0, 则设Hillshade=0.

其中，Zenith_rad是太阳天顶角的的弧度数，Slope_rad是某一点的坡度弧度数，Azimuth_rad是指太阳光线方向角的弧度数，Aspect_rad是某一点的坡向弧度数。

计算山体阴影的照明光源的角度默认是太阳高度角，但是真正计算时，需要用到太阳天顶角，太阳天顶角的计算方法是90°-太阳高度角。所以有如下计算公式：

计算太阳天顶角弧度数:

(2)  Zenith_deg = 90 - Altitude

转换为弧度数:

(3)  Zenith_rad = Zenith * pi / 180.0

计算照明的方向：

照明的方向角是指定的角度数，山体阴影的计算公式需要弧度数。首先，需要将地理上的指南针方向转换为数学上的向右的方向，即向右为起算的方向；其次，需要将角度转换为弧度。

转为数学上的方向角:

(4)  Azimuth_math = 360.0 - Azimuth + 90

注意如果 Azimuth_math >=360.0, 那么:

(5)  Azimuth_math = Azimuth_math - 360.0

转换为弧度:

(6)  Azimuth_rad = Azimuth_math * pi / 180.0

计算坡度和坡向

坡度和坡向是利用一个3*3的窗口在输入影像中访问每个像素，9个像素从左到右、从上到下分别用a-i表示，如图所示：

a
b
c


d
e
f


g
h
i

E像元X方向上的变化率采用如下的算法: 

(7)  [dz/dx] = ((c + 2f + i) - (a + 2d + g)) / (8 * cellsize)

E像元Y方向上的变化率采用如下的算法:

 (8)  [dz/dy] = ((g + 2h + i) - (a + 2b + c)) / (8 * cellsize)

坡度的弧度计算公式，考虑了Z因子（协调Z方向的单位和Xy平面上单位的一个系数）：

(9)  Slope_rad = atan (z_factor * sqrt ([dz/dx]2 + [dz/dy]2)) 

坡向通过下面的方法进行计算：

If [dz/dx] is non-zero:

Aspect_rad= atan2 ([dz/dy], -[dz/dx])

if Aspect_rad< 0 then

Aspect_rad= 2 * pi + Aspect_rad

If [dz/dx] iszero:

if [dz/dy] >0 then

Aspect_rad= pi / 2

else if [dz/dy]< 0 then

Aspect_rad= 2 * pi - pi / 2

else

Aspect_rad = Aspect_rad

山体阴影计算示例：出自arcgis10.0帮助文档。

最后，奉上OpenCL实现的代码：

__kernel void hillshade_kernel( __global const float *pSrcData,
						 __global float *pDestData,const int nWidth,const int nHeight
						 , struct HillshadeOption hillOption)
{
	int j = get_global_id(0);
	int i = get_global_id(1);

	if (j >= nWidth || i >= nHeight)
		return;

	int nTopTmp = i-1;
	int nBottomTmp = i+1;
	int nLeftTep = j-1;
	int nRightTmp = j+1;

	//处理边界情况
	if (0 == i)
	{
		nTopTmp = i;
	}
	if (0 == j)
	{
		nLeftTep = j;
	}
	if (i == nHeight-1)
	{
		nBottomTmp = i;
	}
	if (j == nWidth-1)
	{
		nRightTmp = j;
	}
	__local float afRectData[9];
	afRectData[0] = pSrcData[nTopTmp*nWidth+nLeftTep];
	afRectData[1] = pSrcData[nTopTmp*nWidth+j];
	afRectData[2] = pSrcData[nTopTmp*nWidth+nRightTmp];

	afRectData[3] = pSrcData[i*nWidth+nLeftTep];
	afRectData[4] = pSrcData[i*nWidth+j];
	afRectData[5] = pSrcData[i*nWidth+nRightTmp];

	afRectData[6] = pSrcData[nBottomTmp*nWidth+nLeftTep];
	afRectData[7] = pSrcData[nBottomTmp*nWidth+j];
	afRectData[8] = pSrcData[nBottomTmp*nWidth+nRightTmp];

	const float degreesToRadians = (M_PI_F / 180);

	float dx = ((afRectData[2]+ afRectData[5]*2 + afRectData[8]) -
		(afRectData[0] + afRectData[3]*2 + afRectData[6])) / (8 * hillOption.dbEwres);

	float dy = ((afRectData[6] + afRectData[7]*2 + afRectData[8]) -
		(afRectData[0]+ afRectData[1]*2 + afRectData[2])) / (8 * hillOption.dbNsres);

	//计算坡度(弧度)
	float key = sqrt(dx *dx + dy * dy);
	float dfSlope = atan( hillOption.dfzScaleFactor *  key);

	//计算坡向（弧度）
	float dfAspect = 0;
	if (dx != 0)
	{
		dfAspect = atan2(dy,-dx);
		if (dfAspect < 0)
		{
			dfAspect += 2* M_PI_F;
		}
	}

	if (dx == 0)
	{
		if (dy > 0.0f)
		{
			dfAspect = M_PI_F / 2;
		}

		else if (dy < 0.0f)
			dfAspect = M_PI_F + M_PI_F / 2;
	}

	//将太阳高度和太阳光线角度转换为要求的格式
	float dfZenithDeg = hillOption.dfAltitude;

	float dfAzimuthRad = hillOption.dfAzimuth;

	//最后计算山体阴影值
	float dfHillshade = 255 * (cos(dfZenithDeg)*cos(dfSlope) +
		sin(dfZenithDeg)*sin(dfSlope)*cos(dfAzimuthRad-dfAspect));
	if (dfHillshade < 0)
	{
		dfHillshade = 0;
	}

	if (dfHillshade >= 255)
	{
		dfHillshade = 255;
	}

	pDestData[i*nWidth+j] = (int)(dfHillshade+1/2); 

}