Gdal uses DEM data to calculate hillshade)

0. Preface

Speaking of hillshade, it is actually a method that simulates the contrast between light and shade caused by sunlight shining on the terrain, and then renders the topographic map to make it look stereoscopic. It is often used for map rendering, as shown in table 1, for more information, see [1]. The diagram in table 1 is from reference [1].

Table 1 DEM, mountain shadow, and Application Comparison

DEM image (color rendering)	Mountain shadow calculated from DEM on the left
Paper map without hillshade	Paper map with hillshade
Satellite image without hillshade	Satellite Image with hillshade

From the figure in Table 1, we can see that after the mountain shadow is used to render the map and satellite images, the visual effect has been greatly improved, with a very obvious stereoscopic effect. The following describes how to implement mountain shadow.

1. Introduction to mountain shadow principles

The mountain shadow calculates the hypothetical brightness of the surface by specifying the sun height and other information for each pixel in the grid. By setting the location of the hypothetical light source and calculating the brightness value of each pixel related to the adjacent pixel, the degree of brightness is assumed. For analysis or graphic display, especially after transparency, the "mountain shadow" greatly enhances the surface visualization.

By default, shadows and light are gray gradients related to integers between 0 and 255 (from black to white ). [From references [4, 5].

When calculating the mountain shadow, you need to specify the position of the sun in the sky. The position can be described using the following two parameters.

1. Sun Azimuth

The direction angle of the sun refers to the angle direction of the Sun. It is measured in clockwise direction from 0 degrees to 360 degrees. The azimuth angle of 90 ° is east. The default direction is 315 ° (NW, Northwest ). 1. The yellow circle represents the sun.

2. Sun height Angle

The sun height angle refers to the angle or slope of the sun above the horizon. The unit of height is degrees, ranging from 0 degrees (on the horizon) to 90 degrees (on the top of the head. The default value is 45 degrees. As shown in 2, the yellow circle on the way represents the sun.

Figure 1 sun azimuth Chart 2 sun elevation angle

Ii. Calculation Method

The formula for calculating mountain shadow is shown in formula (1). For more information, see [6]. In formula (1), Zenith indicates the top corner of the sun, azimuth indicates the sun azimuth, slope and aspect indicate the slope and slope respectively. The suffix _ rad indicates that all angles are in radians.

 Hillshade = 255.0 * (COS (zenith_rad) * Cos (slope_rad) + (sin (zenith_rad) * sin (slope_rad) * Cos (azimuth_rad-aspect_rad )))(1)

When the mountain shadow is calculated in the preceding formula, the calculated result may be less than 0. In this case, set this value to 0. The following steps describe how to split the shadow computing process:

1. Calculate the sun Incidence Angle

Specify that the height of the sun must be greater than the angle of the horizon (that is, 0 degrees ). However, the formula used to calculate the mountain shadow value requires that the angle be expressed in radians and the angle be deflection from the vertical direction. Mark the direction perpendicular to the surface (top above the head) as the top of the sky. The corner of the sky is the angle from the day vertex to the Sun's direction, that is, the corner of the Sun's height angle (90 degrees minus the sun's height angle ). To calculate the sun incidence angle, the first step is to convert the sun height angle to the day top angle. The second step is to convert the corner to radians.

Convert the sun height to the day top angle:

Zenith_deg = 90-Altitude(2)

Convert to radian:

Zenith_rad = Zenith * PI/180.0(3)

2. Calculate the sun Azimuth

The mountain shadow formula requires that the azimuth angle be measured in radians. First, convert the corner of the sky from the geographic unit (in the direction of the compass) to the mathematical unit (right angle ). Then convert the azimuth to radians. For details, see formula (4) (5) (6 ).

Method to convert the sun Azimuth:

Azimuth_math = 360.0-azimuth + 90(4)

Note that if,Azimuth_math> = 360.0Then:

Azimuth_math = azimuth_math-360.0(5)

Convert to radian:

Azimuth_rad = azimuth_math * PI/180.0(6)

3. calculate slope and slope

For more information about the calculation method of slope and slope, refer to the previous two blogs. For details, refer to gdal's calculation of slope and Slope Using DEM Data. Refer to [7]. Note that the calculated slope and slope direction should be converted into radians.

So far, all the variables are calculated, and then the mountain shadow value can be obtained using the formula (1. Figure 3 is a calculation result of the mountain shadow. The sun's azimuth angle is 99 degrees and the sun's elevation angle is 33 degrees.

Figure 3 calculation example of mountain shadow

3. algorithm writing

Since the calculation of mountain shadow is still a 3x3 window, we continue to use the previous 3x3 algorithm as the basis. Here we only implement the algorithm part. In the calculation process, we need to consider the pixel size and the elevation unit. If the two units are inconsistent, we need to make a conversion. In addition, we need to exaggerate the elevation, the purpose is to highlight terrain information.

The first is the algorithm structure of mountain shadow. In the algorithm structure, a series of input parameters are processed in advance. For example, the angle directly calculates the sine and Cosine values, this will be used directly in the following algorithms, instead of re-calculation, to improve the efficiency of the algorithm, of course, it is a small detail.

/*** @ Brief mountain shadow algorithm struct */typedef struct {/*! North-South resolution */Double nsres ;/*! East-West resolution */Double ewres ;/*! Sine value of the azimuth angle */doublesin_altradians ;/*! Multiply the cosine of the azimuth value by the Z zoom ratio */doublecos_altradians_mul_z_scale_factor ;/*! Sun height angle */doubleazradians ;/*! Z scaling ratio x/doublesquare_z_scale_factor;} gdalhillshadealgdata;

The next step is the algorithm implementation function. This function is basically written according to the above formula, with no difficulty.

Float gdalhillshadealg; // calculate the slope x = (afwin [0] + afwin [3] + afwin [3] + afwin [6]) -(afwin [2] + afwin [5] + afwin [5] + afwin [8])/psdata-> ewres; y = (afwin [6] + afwin [7] + afwin [7] + afwin [8]) -(afwin [0] + afwin [1] + afwin [1] + afwin [2])/psdata-> nsres; xx_plus_yy = x * x + y * Y; // calculate the slope aspect = atan2 (Y, x); // calculate the mountain shadow value Cang = (psdata-> sin_altradians-psdata-> cos_altradians_mul_z_scale_factor * SQRT (xx_plus_yy) * sin (Aspect-psdata-> azradians)/SQRT (1 + psdata-> square_z_scale_factor * xx_plus_yy); If (Cang <= 0.0) Cang = 1.0; else Cang = 1.0 + (254.0 * Cang); Return (float) Cang ;}

Finally, a function is required to create parameters required for the algorithm. This function is input with six parameters of the image to obtain the image resolution. Z is the zooming ratio of the elevation, sometimes the terrain is relatively flat, and the effect is not very good. You can use the Z value to exaggerate it. The sacle is a conversion coefficient between the horizontal pixel unit and the elevation unit, sometimes, the pixel unit is based on the longitude and latitude, and the elevation unit is usually meter. In this case, you need to set sacle. The last two parameters are the azimuth and elevation angles mentioned above.

void* GDALCreateHillshadeData(double* adfGeoTransform,                                                                    double z,                                                                    double scale,                                                                    double alt,                                                                    double az){         GDALHillshadeAlgData*pData =                   (GDALHillshadeAlgData*)CPLMalloc(sizeof(GDALHillshadeAlgData));          const doubledegreesToRadians = M_PI / 180.0;         pData->nsres= adfGeoTransform[5];         pData->ewres= adfGeoTransform[1];         pData->sin_altRadians= sin(alt * degreesToRadians);         pData->azRadians= az * degreesToRadians;         doublez_scale_factor = z / (8 * scale);         pData->cos_altRadians_mul_z_scale_factor=                   cos(alt* degreesToRadians) * z_scale_factor;         pData->square_z_scale_factor= z_scale_factor * z_scale_factor;         return pData;}

Finally, there is an M file on the Internet that uses MATLAB to calculate the mountain shadow. For more information, see [2]. If you are interested, download and use MATLAB.

Iv. processing result

Finally, I have nothing to say. I will paste the processed image to show you (the azimuth angle is 315 degrees, and the height angle is 45 degrees ). I hope it will be useful to you.

Figure 4 DEM Data

Figure 5 data after DEM Rendering

Figure 5 Effect of scaling an elevation to 1

Figure 5 effect when the elevation scale is 2

Figure 4 effect of zooming to 10

V. References

[1] http://www.gichd.org/fileadmin/pdf/IMSMA/fact-sheets/GICHDFactSheet-Hillshade%20Imagery.PDF (A description of hillshade)

[2] http://www.mathworks.ch/matlabcentral/fileexchange/14863-hillshade/content/hillshade.m (MATLAB Version)

[3] http://download.osgeo.org/ossim/tutorials/pdfs/HillShade.pdf (use instructions in ossim software)

[4] http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/na/009z000000v0000000/

[5] http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/How_Hillshade_works/009z000000z2000000/

[6] burrough, p. A. and McDonell, R. A., 1998. Principles of Geographical Information Systems (oxforduniversity Press, New York), 190 pp.

[7] http://blog.csdn.net/liminlu0314/article/details/8498985 (gdal calculation slope aspect)

[8] http://www.aubreyrhea.net/gis/index.php/tag/hillshade/ (use instructions In ArcMap software)