首页 > 代码库 > DEM山体阴影原理以及算法详解

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表示,如图所示:
abc
def
ghi
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); 

}


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