当前位置:   article > 正文

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实现的代码:

  1. __kernel void hillshade_kernel( __global const float *pSrcData,
  2. __global float *pDestData,const int nWidth,const int nHeight
  3. , struct HillshadeOption hillOption)
  4. {
  5. int j = get_global_id(0);
  6. int i = get_global_id(1);
  7. if (j >= nWidth || i >= nHeight)
  8. return;
  9. int nTopTmp = i-1;
  10. int nBottomTmp = i+1;
  11. int nLeftTep = j-1;
  12. int nRightTmp = j+1;
  13. //处理边界情况
  14. if (0 == i)
  15. {
  16. nTopTmp = i;
  17. }
  18. if (0 == j)
  19. {
  20. nLeftTep = j;
  21. }
  22. if (i == nHeight-1)
  23. {
  24. nBottomTmp = i;
  25. }
  26. if (j == nWidth-1)
  27. {
  28. nRightTmp = j;
  29. }
  30. __local float afRectData[9];
  31. afRectData[0] = pSrcData[nTopTmp*nWidth+nLeftTep];
  32. afRectData[1] = pSrcData[nTopTmp*nWidth+j];
  33. afRectData[2] = pSrcData[nTopTmp*nWidth+nRightTmp];
  34. afRectData[3] = pSrcData[i*nWidth+nLeftTep];
  35. afRectData[4] = pSrcData[i*nWidth+j];
  36. afRectData[5] = pSrcData[i*nWidth+nRightTmp];
  37. afRectData[6] = pSrcData[nBottomTmp*nWidth+nLeftTep];
  38. afRectData[7] = pSrcData[nBottomTmp*nWidth+j];
  39. afRectData[8] = pSrcData[nBottomTmp*nWidth+nRightTmp];
  40. const float degreesToRadians = (M_PI_F / 180);
  41. float dx = ((afRectData[2]+ afRectData[5]*2 + afRectData[8]) -
  42. (afRectData[0] + afRectData[3]*2 + afRectData[6])) / (8 * hillOption.dbEwres);
  43. float dy = ((afRectData[6] + afRectData[7]*2 + afRectData[8]) -
  44. (afRectData[0]+ afRectData[1]*2 + afRectData[2])) / (8 * hillOption.dbNsres);
  45. //计算坡度(弧度)
  46. float key = sqrt(dx *dx + dy * dy);
  47. float dfSlope = atan( hillOption.dfzScaleFactor * key);
  48. //计算坡向(弧度)
  49. float dfAspect = 0;
  50. if (dx != 0)
  51. {
  52. dfAspect = atan2(dy,-dx);
  53. if (dfAspect < 0)
  54. {
  55. dfAspect += 2* M_PI_F;
  56. }
  57. }
  58. if (dx == 0)
  59. {
  60. if (dy > 0.0f)
  61. {
  62. dfAspect = M_PI_F / 2;
  63. }
  64. else if (dy < 0.0f)
  65. dfAspect = M_PI_F + M_PI_F / 2;
  66. }
  67. //将太阳高度和太阳光线角度转换为要求的格式
  68. float dfZenithDeg = hillOption.dfAltitude;
  69. float dfAzimuthRad = hillOption.dfAzimuth;
  70. //最后计算山体阴影值
  71. float dfHillshade = 255 * (cos(dfZenithDeg)*cos(dfSlope) +
  72. sin(dfZenithDeg)*sin(dfSlope)*cos(dfAzimuthRad-dfAspect));
  73. if (dfHillshade < 0)
  74. {
  75. dfHillshade = 0;
  76. }
  77. if (dfHillshade >= 255)
  78. {
  79. dfHillshade = 255;
  80. }
  81. pDestData[i*nWidth+j] = (int)(dfHillshade+1/2);
  82. }


声明:本文内容由网友自发贡献,不代表【wpsshop博客】立场,版权归原作者所有,本站不承担相应法律责任。如您发现有侵权的内容,请联系我们。转载请注明出处:https://www.wpsshop.cn/w/寸_铁/article/detail/768917
推荐阅读
相关标签
  

闽ICP备14008679号