Dem Hillshade principle and algorithm specific explanation

Hillshade principle and algorithm specific explanation

Hillshade fundamentals:

A hillshade is a hypothetical light source that simulates a certain direction and the height of a certain sun. Uses the calculation of neighboring cells to generate a 0-255 grayscale image.

first, The main parameters of the mountain shadow:

1, the incidence angle of the sun light: This angle of the acreage starting point is North direction, in accordance with the clockwise direction, the angle of the range is 0-360 degrees. For example, as seen, the default angle is 315 degrees, Northwest direction, For example as seen in:

2, The Sun height angle: the sun height angle also referred to as the Sun Height. is the angle between the sun light and the local ground plane, the range is 0-90 degrees, The default sun height is 45 degrees, for example to See:

second, the method of calculating the mountain shadow

The calculation formula for the Hillshade is as follows

(1) Hillshade = 255.0 * ((cos (zenith_rad) * cos (slope_rad)) +

(sin (zenith_rad) * sin (slope_rad) * cos (azimuth_rad- Aspect_rad)))

Assuming Hillshade < 0, the hillshade=0is Set.

The Zenith_rad is the number of radians in the Sun's Zenith Angle. Slope_rad is the number of slope radians for a point, and Azimuth_rad is the number of radians in the direction of the Sun's Rays. Aspect_rad is the number of slope radians for a point.

The angle of the lighting source that calculates the hillshade is the sun height angle by default, But when it is actually calculated. The Sun Zenith Angle is required, and the solar Zenith angle is calculated by 90°-the height of the Sun. So there are for example the following calculation Formula:

calculate the number of radians of the sun Zenith Angle :

zenith_deg Altitude

Convert to radians Number:

Zenith_rad Zenith * pi/180.0

To calculate the direction of illumination:

The direction angle of the illumination is the specified number of angles, and the calculation formula of the Hillshade requires Radians.

first, you need to convert the compass direction of the geography to a mathematically right direction. That is, the right is the starting direction; You need to convert the angle to Radians.

To a mathematically oriented angle:

Azimuth_math Azimuth + 90

Note that assuming azimuth_math >=360.0, then:

Azimuth_math Azimuth_math -360.0

Convert to Radians:

Azimuth_rad Azimuth_math * pi/180.0

Calculate slope and aspect

Slope and aspect is the use of a 3*3 form to access each pixel in the input image, 9 pixels from left to right, from top to bottom, respectively, with a-i, see:

 
   
  

  
    
   
 
    
 
     
A
 
     
B
 
     
C
 
    
 
    
 
     
D
 
     
E
 
     
F
 
    
 
    
 
     
G
 
     
H
 
     
I
 
    
 
   
 
   
  

the rate of change in E-cell X-direction is based on an algorithm such as the Following:

(7) [dz/dx] = ((c + 2fi)-(a + 2dgcellsize)

The rate of change in the y direction of e cells is based on an algorithm such as the Following:

(8) [dz/dy] = ((g + 2hi)-(a + 2bccellsize)

The arc of the slope is calculated by considering the Z-factor (a factor that coordinates the units in the z-direction and the units on the XY plane):

Slope_rad = atan (z_factor

The aspect is calculated by the following methods:

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

Hillshade calculation demo Sample: from the arcgis10.0 help Document.

finally, the code implemented in Opencl:

__kernel void Hillshade_kernel (__global const float *psrcdata, __global float *pdestdata,const int nwidth,const int NHeig ht, struct Hillshadeoption hilloption) {int j = get_global_id (0); int i = get_global_id (1); if (j >= nwidth | | I >= nH Eight) return;int ntoptmp = i-1;int nbottomtmp = i+1;int Nlefttep = j-1;int nrighttmp = j+1;//processing boundary condition 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 = ((afrectd ata[6] + afrectdata[7]*2 + afrectdata[8])-(afrectdata[0]+ afrectdata[1]*2 + afrectdata[2])/(8 * hilloption.dbnsres);// Calculate slope (radians) Float key = sqrt (dx *dx + dy * dy), float dfslope = atan (hilloption.dfzscalefactor * key);//calculate aspect (radians) float Dfaspec t = 0;if (dx! = 0) {dfaspect = atan2 (dy,-dx), if (dfaspect < 0) {dfaspect + = 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;} Convert the sun height and the sun light angle to the required format float dfzenithdeg = hilloption.dfaltitude;float Dfazimuthrad = hilloption.dfazimuth;// The last computed Hillshade value 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 Hillshade principle and algorithm specific explanation