Calculation of slope and Slope Using DEM Data in gdal

0. Preface

I have previously written a blog titled gdal-based common 3x3 Template operator function. The website is http://blog.csdn.net/liminlu0314/article/details/8%156. At that time, it was necessary to write a function to calculate the slope and slope based on this function. As this time has been busy with other things, this is a drag-and-go task. I will make it up for you today.

I. Introduction

The slope is the steep degree of the surface unit. Generally, the ratio of the vertical height h and the horizontal distance L of the slope is called the slope ratio, which is expressed by the letter I. [That is, the tangent of the slope angle (writeable: I = tan slope angle )]. There are four methods to indicate the slope: percentage, degree, password, and score. The percentage and degree methods are commonly used (1 ).

Figure 1 two slope Representation Methods

1. Percentage Method

It indicates the most common method of slope, that is, the percentage of the distance between two points and its horizontal distance. The formula is as follows:

Slope = (Height Difference/horizontal distance) x100 %

When the percentage is used, that is, I = H/L × 100%

For example, if the slope is 3%, the horizontal distance is 100 meters, and the vertical direction is increased (decreased) by 3 meters. If the slope is 1%, the horizontal distance is 100 meters, and the vertical direction is increased (decreased) by 1 meter. And so on!

2. Level Method

The degree is used to represent the slope, and the inverse trigonometric function is used for calculation. The formula is as follows:

Tan α (slope) = gradient/horizontal distance

Therefore, α (slope) = arctan (Height Difference/horizontal distance)

If the elevation increment percentage is treated as the elevation increment divided by the horizontal increment and multiplied by 100, the elevation increment percentage can be better understood. Consider the triangle belowB. When the angle is 45 degrees, the elevation increment is equal to the horizontal increment, so the elevation increment percentage is 100%. TriangleCAs shown in, when the slope angle is close to the right angle (90 degrees), the elevation increment percentage begins to approach infinity.

Aspect refers to the orientation of the terrain slope. The slope direction is used to identify the downhill direction with the highest change rate from each pixel to its adjacent pixel. The slope direction can be considered as the slope direction. The aspect direction is an angle, which will be measured clockwise. The angle ranges from 0 (East) to 360 (still east), that is, the complete circle. Flat areas that do not have a downhill direction are assigned a value of-1. 2.

Figure 2 value of slope direction

Ii. Slope and slope Calculation Method

Generally, the method of fitting a surface is used for slope calculation. A fitting surface generally uses a quadratic surface, that is, a window of 3 × 3, as shown in 3. The center of each window is a path. In Figure 3, the calculation formula for the slope and slope of center e is as follows:

In the above formula, slope is the slope, aspect is the slope of the slope, slopewe is the slope of the X direction, and slopesn is the slope of the Y direction.

Figure 3 a 3 × 3 window

Slopewe and slopesn can be calculated using the following common methods:

L algorithm 1

L algorithm 2

L algorithm 3

L algorithm 4

The cellsize in the preceding formula is the Interval Length of the Grid DEM. Algorithm 1 has the highest accuracy and computing efficiency, followed by algorithm 2. ERDAS Imagine uses algorithm 4 and ArcMap uses algorithm 2. Algorithm 2 is also used this time.

Iii. Compiling of slope and slope Algorithms

The function for calculating the slope and slope is dependent on the function in the previous blog gdal-based general 3x3 template function, therefore, only the algorithm code for calculating the slope and slope is provided here. For details about calling this code, refer to gdal-based general 3x3 template function. Through the blog "gdal-based common 3 × 3 Template Function", we know that to implement an algorithm, we need to write three parts: Structure of algorithm parameters; a callback function implemented by an algorithm and a function that creates a struct pointer of an algorithm parameter.

1. Slope

The first is to define the structure of the slope algorithm. The parameters required for calculating the slope include the DEM grid size, that is, the resolution between East and West, and the North and South, the scale of the elevation, and the representation of the slope, that is, the percentage or degree:

/*** @ Brief slope algorithm Data Structure */typedef struct {/*! North-South resolution */Double nsres ;/*! East-West resolution */Double ewres ;/*! Scaling ratio */double scale ;/*! Slope mode */INT slopeformat;} gdalslopealgdata;

Next is the callback function of the slope algorithm. This function is used to calculate the slope. You can write it according to the above formula. Finally, you need to convert it based on the specified expression of the slope:

float GDALSlopeAlg (float* afWin, float fDstNoDataValue,void* pData){         const doubleradiansToDegrees = 180.0 / M_PI;         GDALSlopeAlgData*psData = (GDALSlopeAlgData*)pData;         double dx,dy, key;          dx =((afWin[0] + afWin[3] + afWin[3] + afWin[6]) -                   (afWin[2]+ afWin[5] + afWin[5] + afWin[8])) / psData->ewres;          dy =((afWin[6] + afWin[7] + afWin[7] + afWin[8]) -                   (afWin[0]+ afWin[1] + afWin[1] + afWin[2])) / psData->nsres;          key = (dx *dx + dy * dy);          if(psData->slopeFormat == 1)                   return(float)(atan(sqrt(key) / (8*psData->scale)) * radiansToDegrees);         else                   return(float)(100*(sqrt(key) / (8*psData->scale)));}

Finally, the function for creating a slope struct is as follows:

void* GDALCreateSlopeData(double* adfGeoTransform,       double scale, int slopeFormat){         GDALSlopeAlgData*pData =                   (GDALSlopeAlgData*)CPLMalloc(sizeof(GDALSlopeAlgData));          pData->nsres= adfGeoTransform[5];         pData->ewres= adfGeoTransform[1];         pData->scale= scale;         pData->slopeFormat= slopeFormat;         return pData;}

2. aspect

The first is to define the structure of the slope direction algorithm. The parameter of the slope direction is very simple, that is, whether the azimuth is used for representation, that is, if this value is set to true, the slope calculation result is expressed according to the angle in Figure 2. If it is false, the calculation result is as much as possible:

/*** @ Brief slope aspect algorithm Data Structure */typedef struct {/*! Azimuth */intbangleasazimuth;} gdalaspectalgdata;

The next step is the callback function of the slope direction algorithm. According to the above formula, there is no difficulty. Note that if the value of bangleasazimuth is true, you need to convert the calculation angle. The conversion result is 0 indicating east, 90 indicating North, and 180 indicating west, 270 indicates South:

Float gdalaspectalg (float * afwin, float metadata, void * pdata) {const metadata = m_pi/180.0; gdalaspectalgdata * psdata = (gdalaspectalgdata *) pdata; double dx, Dy; float aspect; DX = (afwin [2] + afwin [5] + afwin [5] + afwin [8]) -(afwin [0] + afwin [3] + afwin [3] + afwin [6]); DY = (afwin [6] + afwin [7] + afwin [7] + afwin [8]) -(afwin [0] + afwin [1] + afwin [1] + afwin [2]); aspect = (float) (atan2 (dy,-dx) /degreestoradians); If (dx = 0 & DY = 0) {/* flat area */aspect = fdstnodatavalue; // or here it is-1} else if (psdata-> bangleasazimuth) {If (aspect> 90.0) aspect = pai0000f-aspect; else aspect = 90.0f-aspect ;} else {If (aspect <0) aspect + = 360.0;} If (aspect = 360.0) aspect = 0.0; returnaspect ;}

Finally, the function for creating the slope direction struct is as follows:

void* GDALCreateAspectData(int bAngleAsAzimuth){         GDALAspectAlgData*pData =                   (GDALAspectAlgData*)CPLMalloc(sizeof(GDALAspectAlgData));          pData->bAngleAsAzimuth= bAngleAsAzimuth;         return pData;}

Iv. Result Display

As shown in original DEM data 4, Figure 5 shows the calculated slope value, and Figure 6 shows the calculated slope value. The following three images are displayed after stretching. Because the DEM data is 16 bit data, the calculated slope and slope data is a 32-bit floating point data. No hierarchical rendering is performed on the image, and a color rendering is performed on the slope and slope data processed by ArcMap. The effect is better. Of course, this will not be discussed here.

Figure 4 a dem (SRTM data) in New South Wales, Australia)

Figure 5 slope Calculation Result

Figure 6 calculation result of slope direction

V. References:

[1] xiangling Liu, Glenn Otto (July 2013). The Determinants of supply elasticity in sysydney housing market, Working Paper, University of New South Wales

[2] http://help.arcgis.com/zh-cn/arcgisdesktop/10.0/help/index.html#//009z000000v2000000

[3] http://baike.baidu.com/view/1353583.htm

[4] http://help.arcgis.com/zh-cn/arcgisdesktop/10.0/help/index.html#//009z000000tr000000

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

[6] http://baike.baidu.com/view/762866.htm