Slope and aspect calculation method applied to geographic coordinate digital elevation model

Abstract

The present invention discloses a slope and aspect calculation method applied to a digital elevation model of geographic coordinates, wherein a local coordinate system of a surface model is established with DEM data under geographic coordinates, and the elevation change rate is calculated according to the elevation value of each pixel in the coordinate system and its coordinate value, and then the slope and aspect value of each pixel is calculated according to the elevation change rate. The slope and aspect calculation method of the present invention solves the problem that the calculation accuracy of the slope and aspect in the prior art is insufficient, and the need for high-precision terrain analysis of DEM data with a large spatial range based on geographic coordinates cannot be met.

Description

Technical Field

[0001] The invention belongs to the technical field of digital terrain analysis methods, and in particular relates to a slope and aspect calculation method applied to a geographic coordinate digital elevation model. Background Art

[0002] Slope is the angle between the normal direction of a point on the ground and the vertical direction, that is, the inclination of the slope (surface) or slope profile (line) relative to the horizontal ground; slope aspect is the angle between the terrain slope and the north direction in the clockwise direction. Slope and aspect are important terrain parameters, which play a key role in the movement of surface materials under the action of gravity, microclimate characteristics, and slope vegetation coverage characteristics. They are also important factors affecting landscape characteristics and are widely used in soil and water conservation planning, soil erosion calculation, land use planning and other fields. Therefore, obtaining accurate surface slope and aspect data is a prerequisite for high-quality geological analysis.

[0003] Traditional slope and aspect algorithms for plane coordinates include projection conversion method and fixed pixel size method. The projection conversion method converts the digital elevation model (DEM) in geographic coordinates into plane coordinates, calculates the slope and aspect using the plane coordinate slope and aspect algorithm, and then converts it back to geographic coordinates. This method will change the elevation value and spatial position of the original DEM pixel during the projection process, affecting the accuracy of the slope and aspect finally obtained; in addition, the two projection processes also seriously reduce the calculation efficiency, especially when a large amount of data needs to be processed. The fixed pixel size method assumes that the size of each pixel of the DEM data in geographic coordinates is equal in length units, and then uses the plane coordinate slope and aspect algorithm to calculate the slope and aspect of the DEM data in geographic coordinates. The premise for using this method is that the latitudinal span of the spatial range of the DEM data cannot be too large, the pixel size cannot be too large, and the DEM cannot be located in a high-latitude area. Otherwise, the size of each pixel in the DEM data will differ too much in length units, which will affect the accuracy of the calculated slope and aspect. In addition, there is currently no ready-made software tool that integrates this method, and users need to implement it through self-programming.

[0004] Therefore, the calculation accuracy and efficiency of slope and aspect in geographic coordinates under the existing technology are insufficient, and cannot meet the needs of high-precision terrain analysis based on large spatial range DEM data based on geographic coordinates. It is urgent to explore a slope and aspect algorithm that can be directly applied to geographic coordinate DEM to provide important data and technical support for the high-quality development of soil and water conservation and related industries. Summary of the invention

[0005] The purpose of the present invention is to provide a slope and aspect calculation method applied to a geographic coordinate digital elevation model, so as to solve the problem that the accuracy of slope and aspect calculation in the prior art is insufficient and cannot meet the needs of high-precision terrain analysis based on a large spatial range DEM of geographic coordinates.

[0006] The technical scheme adopted by the present invention is: a slope and aspect calculation method applied to a geographic coordinate digital elevation model, firstly taking any pixel of DEM data under geographic coordinates as the center, obtaining a local pixel window covering 3 rows×3 columns (simplified as 3×3 below), then projecting the local pixel window onto a plane to establish a local plane rectangular coordinate system, then calculating the elevation change rate of the central pixel according to the elevation value and coordinate value of each pixel, and incorporating it into a slope and aspect formula to calculate the slope and aspect of the central pixel; repeating the above steps for all pixels of the DEM data to complete the calculation of the slope and aspect of the entire DEM data.

[0007] The technical solution adopted by the present invention is also characterized in that:

[0008] Furthermore, the slope and aspect calculation method applied to the geographic coordinate digital elevation model is implemented in the following steps:

[0009] Step 1: Take a pixel in the DEM data under geographic coordinates as the center and obtain a 3×3 local pixel window;

[0010] Step 2: Project the local pixel window onto the plane, and establish a local plane rectangular coordinate system with the central pixel as the origin;

[0011] Step 3, calculate the elevation change rate of the central pixel according to the elevation and coordinate values ​​of each pixel in the local window, and substitute them into the slope aspect formula to obtain the slope aspect value of the central pixel;

[0012] Step 4: Perform the above steps for each pixel of the DEM data under the geographic coordinates in turn to complete the slope and aspect calculation of the entire DEM data.

[0013] Furthermore, the local pixel window of step 1 is located on the earth ellipsoid, wherein the boundaries of each pixel are adjacent to each other; the upper boundary and the lower boundary of each pixel in the same row are respectively located on the same latitude circle; the left boundary and the right boundary of each pixel in the same column are respectively located on the same meridian circle; each pixel covers a certain surface range on the surface of the earth ellipsoid, and the size of the range is represented by the size of the pixel boundary; the size of the pixel boundary under the geographic coordinate is expressed in the angular unit "degree, minute, second", and the sizes of the four boundaries are equal.

[0014] Furthermore, step 2 is specifically implemented according to the following steps:

[0015] Step 2.1: Use the Evans surface model to simulate the terrain surface of a 3×3 local pixel window. The Evans surface formula is shown in (1):

[0016]

[0017] Among them, x and y represent the horizontal and vertical coordinates of each pixel in the 3×3 window respectively, z represents the elevation value of the pixel, r is the quadratic term coefficient of x, t is the quadratic term coefficient of y, s is the quadratic term coefficient of xy, p is the linear term coefficient of x, q is the linear term coefficient of y, and f is the constant term. These coefficients and constant terms are used to control the shape of the surface, among which p and q have the same meaning as p and q in the slope and aspect formula;

[0018] Step 2.2, since the Evans surface model can only be used for plane coordinates, it is necessary to project the local pixel window of the surface on the ellipsoid in step 1 onto the plane. The specific projection method is as follows: make chords of the arc formed by the outer boundaries of the 8 pixels at the non-center position of the 3×3 local pixel window. These chords intersect each other at the outer vertices of the 4 pixels on the diagonal. Since the chords on each latitude circle are parallel, these chords must be located on the same plane. This plane is used as the projection plane of the local surface, and each chord is used as the outer boundary of each pixel; at the same time, the intersection of the latitude circle plane where the other upper and lower boundaries of each pixel are located and the meridian circle plane where the left and right boundaries are located with the above-mentioned projection plane is used as the other boundaries of each pixel on the projection plane;

[0019] Step 2.3: In the 3×3 local pixel window on the projection plane created in step 2.2, a plane rectangular coordinate system is established; wherein the origin of the coordinate system is located at the center of the central pixel, the row direction (east-west direction) is the x-axis, and the column direction (north-south direction) is the y-axis.

[0020] Furthermore, step 3 is specifically implemented according to the following steps:

[0021] Step 3.1: The coordinate values ​​of each pixel, that is, their distances to the x and y axes, are represented by w1~w 12 It can be expressed as follows: par Formula (2) and the chord length C of the warp coil mer Formula (3) approximately expresses:

[0022]

[0023] Among them, a and b are the major and minor semi-axes of the earth ellipsoid, ω is the pixel size of the DEM in geographic coordinates, θ is the eccentricity angle between two adjacent pixels on the latitude circle, θ1 and θ2 are the eccentricity angles between two adjacent pixels on the meridian circle, and the eccentricity angle and latitude are The relationship between can be expressed as formula (4):

[0024]

[0025] In a 3×3 local pixel window, the latitude difference between two adjacent pixels on the meridian circle is a pixel size ω, so the latitude of the top row of pixels in the window is and the latitude of the bottom row The latitude of the middle row of pixels can be It is expressed as formula (5):

[0026]

[0027] In the projection plane rectangular coordinate system, the x-coordinate of each pixel on the left side of the y-axis is negative, and the x-coordinate of each pixel on the right side is positive; the y-coordinate of each pixel above the x-axis is positive, and the y-coordinate of each pixel below is negative;

[0028] Substituting formula (5) and (4) into formula (2) and (3), we can obtain w1~w 12 Formulas (6) to (10):

[0029]

[0030] in, is the latitude of the central pixel of the 3×3 local pixel window; given the DEM data in geographic coordinates, a, b, and ω are all known constants, and It changes with the location of the pixel and is easy to obtain based on the pixel size and spatial range attributes of the DEM data;

[0031] The elevation and coordinate values ​​of each pixel are expressed as follows: (z1, -w1, w7), (z2, 0, w8), (z3, w2, w9), (z4, -w3, 0), (z5, 0, 0), (z6, w4, 0), (z7, -w5, -w10), (z8, 0, -w11), (z9, w6, w12);

[0032] Step 3.2: Substitute the coordinates and elevation values ​​of each pixel in step 3.1 into the Evans surface formula, and let w1=w2=g1, w3=w4=g2, w5=w6=g3, w7=w8=w9=f1, w 10 =w 11 =w 12 =f2, and the calculation formulas of p and q in the geographic coordinate system are obtained based on least squares regression as shown in (11) and (12):

[0033]

[0034] Step 3.3: Substitute the calculated results of p and q in step 3.2 into the slope calculation formula (13) and the aspect formula (14) to obtain the slope and aspect value of the central pixel:

[0035]

[0036] The present invention provides a slope and aspect calculation method applied to a geographic coordinate digital elevation model. Compared with the traditional method, the present invention has the following beneficial effects:

[0037] 1. The method of the present invention avoids the changes in pixel elevation, slope aspect value and pixel position caused by the two projection processes in the traditional projection method. These changes reduce the accuracy of the slope aspect. At the same time, since no additional projection conversion process is required, the calculation efficiency is improved.

[0038] 2. The method of the present invention avoids the defect that the accuracy difference varies greatly when the traditional fixed pixel size method is applied to DEMs of different spatial ranges, different pixel sizes, and different latitudes.

[0039] 3. The method of the present invention can perform calculations with other geographic information raster images in the geographic coordinate system without losing accuracy (avoiding unnecessary projection conversion). BRIEF DESCRIPTION OF THE DRAWINGS

[0040] Figure 1 is an operation flow chart of the method of the present invention;

[0041] Figure 2 It is a 3×3 local pixel window composed of a certain pixel of DEM data and its neighboring pixels under the geographic coordinates of the method of the present invention and its projection plane rectangular coordinate system diagram;

[0042] Figure 3a is a side view of the Gaussian surface in Example 1 of the present invention;

[0043] Figure 3b is a Gaussian surface contour map in Example 1 of the present invention;

[0044] Figure 4a It is a sinusoidal curve projection diagram with the central meridian of the northern hemisphere being 108°E in Example 1 of the present invention;

[0045] Figure 4b is a sinusoidal curve projection diagram near the equator in Example 1 of the present invention;

[0046] Figure 4c It is a sinusoidal projection diagram near the pole in Example 1 of the present invention;

[0047] Figure 5 It is a schematic diagram of the true value, estimated value and difference of the slope of the Gaussian surface near the equator and the pole in geographic coordinates in Example 1 of the present invention;

[0048] Figure 6 It is a schematic diagram of the true value, estimated value and difference of the aspect of the Gaussian surface near the equator and the pole in geographic coordinates in Example 1 of the present invention;

[0049] Figure 7Schematic diagram of three experimental areas with different elevation fluctuations in Example 2 of the present invention;

[0050] Figure 8 It is a comparison diagram of SRTM slope and aspect data in low undulating hilly areas obtained by three methods in Example 2 of the present invention and Comparative Examples 1 and 2;

[0051] Fig. 9 It is a comparison diagram of SRTM slope and aspect data of the rolling hills and gullies area obtained by three methods in Example 2 of the present invention and Comparative Examples 1 and 2;

[0052] Fig.10 is a comparison diagram of SRTM slope and aspect data in steep mountain areas obtained using three methods in Example 2 of the present invention and Comparative Examples 1 and 2;

[0053] Fig.11 is the slope map of the Muli mining area obtained using the method of the present invention in Example 3 of the present invention;

[0054] Fig.12 is the slope aspect map of the Muli mining area obtained using the method of the present invention in Example 3 of the present invention;

[0055] Fig.13 This is the distribution map of hydraulic and freeze-thaw erosion intensity in the Muli mining area obtained using the method of the present invention in Example 3 of the present invention. DETAILED DESCRIPTION

[0056] In order to make the purpose, technical solution and advantages of the present invention clearer and more understandable, the technical solution of the present invention will be clearly and completely described below in conjunction with the accompanying drawings.

[0057] The algorithm for calculating the slope and aspect of the DEM under geographic coordinates in the present invention has been written into an ArcGIS Toolbox tool using Python language and can be used under the ArcGIS platform.

[0058] The present invention discloses a slope and aspect calculation method applied to a geographic coordinate digital elevation model. The local coordinate system of the surface model is established with DEM data under geographic coordinates. The elevation change rate is calculated according to the elevation value of each pixel in the coordinate system and its coordinate value. Then, the slope and aspect value of each pixel is calculated according to the elevation change rate. The specific operation process is as follows. Figure 1 shown.

[0059] The formula for calculating the slope of a point on the surface is: slope = arctan[(p 2 +q 2 ) 1 / 2], the slope aspect calculation formula is: aspect=180–arctan(q / p)+90p / |p|; From the formula, we can see that the value of slope aspect depends on two key parameters p and q, and p and q represent the elevation change rate of the point on the surface along the east-west and north-south directions respectively; therefore, whether p and q are accurately calculated directly affects the accuracy of slope aspect data.

[0060] The traditional p and q algorithms are aimed at DEM data under plane projection coordinates. The units of the pixel coordinate parameters contained in their formulas are length units such as "meters", and the distance between pixels is a straight line on the plane.

[0061] Geographic coordinates are based on the earth ellipsoid and use longitude and latitude to divide the surface of the ellipsoid. In the DEM data under this coordinate system, each pixel represents a local surface of the earth ellipsoid, and the distance between two adjacent pixels is the length of the arc on the ellipsoid surface, expressed in "degrees, minutes, seconds" angular units. Therefore, the slope and aspect formula based on plane coordinates cannot be used for DEM data under geographic coordinates.

[0062] To calculate the slope and aspect of a pixel in the DEM data under geographic coordinates, we first need to establish a local coordinate system containing 3×3 pixels with the pixel as the center, and calculate p and q based on the elevation value and coordinate position of each pixel in the coordinate system.

[0063] The local coordinate system of a pixel and its neighboring pixels in DEM data under geographic coordinates is as follows: Figure 2 As shown, the gray solid curve represents the position and shape of each pixel on the ellipsoid, and the boundary of each pixel is located on different longitude and latitude circles.

[0064] To calculate p and q of the central pixel z5, the Evans surface model is used to simulate the local terrain surface of the 3×3 pixel. The Evans surface is shown in (1):

[0065]

[0066] Among them, x and y represent the horizontal and vertical coordinates of each pixel in the 3×3 window respectively, z represents the elevation value of the pixel, r is the quadratic term coefficient of x, t is the quadratic term coefficient of y, s is the quadratic term coefficient of xy, p is the linear term coefficient of x, q is the linear term coefficient of y, and f is the constant term. These coefficients and constant terms are used to control the shape of the surface, among which p and q have the same meaning as p and q in the slope and aspect formula;

[0067] Substitute the elevation values ​​and coordinate values ​​of the 9 pixels in the 3×3 window into the Evans surface formula, and then use the least squares method to express p and q as functions of z, x, and y;

[0068] z1~z9 are the elevation values ​​of each pixel, which can be directly obtained from the DEM pixel value, while x and y need to be calculated through the major and minor semi-axis parameters of the earth ellipsoid, the pixel size of the DEM and the latitude of the pixel location.

[0069] Since the Evans surface equation is only applicable to plane coordinates (because the units of x, y, and z are the same), it is necessary to project the local surface coordinate system on the ellipsoid onto the plane coordinates. The projection process is as follows: first, the upper boundary latitude circle arc segment of the pixel corresponding to the 3×3 local pixel window z1, z2, and z3 on the ellipsoid is projected onto its chord, and the lower boundary latitude circle of the pixel corresponding to z7, z8, and z9 is projected onto its chord. The left boundary of the pixel corresponding to z1, z4, and z9 is projected onto its chord through the coil arc segment, and the right boundary of the pixel corresponding to z3, z6, and z9 is projected onto its chord through the coil arc segment. These chords are connected end to end and are located on the same plane. Then the other boundaries of each pixel are also projected onto the plane. The projection results are shown as follows: Figure 2 As shown by the grey dashed line in .

[0070] According to the principles of geodesy, in an area of ​​100 km 2 Whether measuring horizontal distance or horizontal angle within the range, the influence of the earth curvature can be ignored. Therefore, when the boundary length of the 3×3 local window does not exceed 10km, that is, the pixel size does not exceed about 10 / 3km, the error caused by the local coordinate projection can be ignored; this pixel size corresponds to the minimum angular interval of about 1.8′ on the earth ellipsoid, and in actual terrain analysis and other applications, the commonly used DEM data pixel size is much higher than this value (such as the pixel size of SRTM elevation data is 1″ and 3″, and the pixel size of ASTER GDEM data is 1″), so the above projection method is feasible.

[0071] On the plane formed by the gray dotted line, a local plane coordinate system is established with the pixel center corresponding to z5. The x-axis is the east-west direction, and the y-axis is the north-south direction. The coordinate value of each pixel, that is, their distance to the x-axis and the y-axis, is represented by w1~w 12 It can be expressed as follows: par Formula (2) and the chord length C of the warp coil mer Formula (3) calculates:

[0072]

[0073] Among them, a and b are the major and minor semi-axes of the ellipsoid (length units), ω is the pixel size of the DEM in geographic coordinates (angle unit), θ is the eccentricity angle between two adjacent pixels on the latitude circle, θ1 and θ2 are the eccentricity angles between two adjacent pixels on the meridian circle, and the eccentricity angle is related to latitude. The relationship between can be expressed as formula (4):

[0074]

[0075] In a 3×3 window, the latitude difference between two adjacent pixels on the meridian circle is a pixel size ω, so the latitude of the top row of pixels in the window is and the latitude of the bottom row The latitude of the middle row of pixels can be It is expressed as formula (5):

[0076]

[0077] In the projection plane coordinates of the 3×3 window, the x-coordinates of the pixels on the left side of the y-axis are negative, and the x-coordinates of the pixels on the right side are positive; the y-coordinates of the pixels above the x-axis are positive, and the y-coordinates of the pixels below are negative;

[0078] Combining formulas (2), (3), (4), and (5), we can get w1~w 12 Formulas (6) to (10):

[0079]

[0080] in, is the latitude (angular unit) of the central pixel of the 3×3 local window; given the DEM data in geographic coordinates, a, b, and ω are all known constants, and It changes with the location of the pixel and is easy to obtain based on the pixel size and spatial range attributes of the DEM;

[0081] The elevation and coordinate values ​​of each pixel are expressed as follows: (z1, -w1, w7), (z2, 0, w8), (z3, w2, w9), (z4, -w3, 0), (z5, 0, 0), (z6, w4, 0), (z7, -w5, -w10), (z8, 0, -w11), (z9, w6, w12);

[0082] At this time, w1~w 11 is a function of φ, let w1=w2=g1, w3=w4=g2, w5=w6=g3, w7=w8=w9=f1, w 10 =w 11 =w 12 =f2, and substitute them and z1~z9 into the Evans surface formula, and obtain the calculation formulas of p and q in the geographic coordinate system according to the least squares regression as shown in (11) and (12):

[0083]

[0084] Substituting the above p and q formulas into the slope and aspect formulas (13) and (14) can obtain the slope and aspect value of the central pixel;

[0085]

[0086] The above process is carried out for each pixel of the DEM data under the geographic coordinates in turn to complete the slope and aspect calculation of the entire DEM data.

[0087] The present invention is further described below in conjunction with examples and drawings, but the embodiments of the present invention are not limited thereto.

[0088] Example 1

[0089] Calculates the slope and aspect of a mathematically simulated terrain surface.

[0090] Step 1: Use Gaussian synthetic surface to simulate the terrain surface;

[0091] The Gaussian surface formula is:

[0092]

[0093] Where x and y represent the position of any point on the surface on the horizontal plane, z represents the elevation of the point, and A, B, C, m, and n are parameters that control the shape of the Gaussian surface;

[0094] Let A=90, B=300, C=10, m=750, n=750, and control the value range of x and y to [-1500, 1500]. The Gaussian surface obtained is as follows Figure 3a and Figure 3b As shown;

[0095] It should be noted that the Gaussian surface, like the Evans surface model, can only be established on plane coordinates.

[0096] Step 2: Place the Gaussian surface in the equatorial and polar regions;

[0097] The unit of DEM pixel size under geographic coordinates is "degree", and the pixel boundary is parallel to the meridian and latitude circles; when the same angle interval on the meridian and latitude circles is converted into the corresponding length unit, their length will change with latitude; near the equator, the lengths corresponding to the same angle interval on the meridian and latitude circles are basically equal, and as the latitude increases toward the North and South Poles, the difference between the two will gradually increase;

[0098] In addition, the length corresponding to the same angular interval on the meridian circle changes little from the equator to the pole, while the length corresponding to the same angular interval on the latitude circle gradually decreases from the equator to the pole; therefore, when calculating the slope and aspect of the Gaussian surface in geographic coordinates, two boundary conditions are considered to place the Gaussian surface: one is near the equator, and the other is near the pole;

[0099] The center point of the Gaussian surface near the equator is placed at (108°E, 0°N), and the center point of the Gaussian surface near the pole is placed at (108°E, 83°N);

[0100] From the p and q formulas, it can be seen that the choice of the longitude of the center point does not affect the calculation of the slope and aspect, so it can be set randomly.

[0101] Step 3: Select the sinusoidal curve projection (belonging to plane coordinates) as the slope and direction of the bridge Gaussian surface in geographic coordinates;

[0102] The sinusoidal projection has lower geometric distortion at both high and low latitudes. Its geometric structure is closer to the geographic coordinates before projection than other common projections. Specifically, the latitude circles after projection are perpendicular to the central meridian, which helps to eliminate the influence of projection distortion on the accuracy of slope and aspect calculation.

[0103] The geometry of the sinusoidal projection is given by Figure 4a , Figure 4b , Figure 4c As shown; along the central meridian, the projected longitude and latitude grid near the equator is approximately a square, and the projected longitude and latitude grid near the pole is approximately a rectangle.

[0104] Step 4, calculate the slope and aspect of the Gaussian surface near the equator and the poles in geographic coordinates;

[0105] Set the units of the x, y, and z values ​​of the Gaussian surface to "meters", and then locate its origin (0, 0) near the equator (108°E, 0°N), with the x-axis coplanar with the equator to the east and the y-axis coplanar with the central meridian to the north;

[0106] Within the Gaussian surface near the equator, based on the geographic coordinates of the WGS1984 ellipsoid, a sample point set with an interval of 1″ (consisting of 99×99 sample points) is arranged, and then the equivalent position distribution of the sample point set in the sinusoidal projection is obtained through projection transformation; the first-order partial derivatives of the Gaussian surface at the sample point positions in the x and y directions (i.e., the elevation change rate of the Gaussian surface at the sample point positions along the x and y directions) are calculated, and substituted into the slope aspect formula to calculate the true value of the slope aspect of these sample points;

[0107] Substitute the x and y coordinates of these sample points in the sinusoidal projection into the Gaussian surface formula to obtain the elevation value set of these sample point locations; then use these elevation point sets and the geographic coordinates corresponding to their locations to create a Gaussian surface DEM with a pixel size of 1″, and combine the p and q formulas under geographic coordinates to calculate the slope and aspect estimation value of the DEM;

[0108] The origin (0, 0) of the Gaussian surface near the equator is moved to the vicinity of the pole (108°E, 93°N), and the x-axis is kept coplanar with the latitude circle to the east, and the y-axis is kept coplanar with the central meridian to the north; based on the geographic coordinates of the WGS1984 ellipsoid, a sample point set with an interval of 1″ is arranged, but this time the sample point set consists of 99×809 points; the calculation method of the true value of the slope and aspect of these sample point locations and the estimated value of the slope and aspect of the Gaussian surface DEM under the corresponding geographic coordinates is consistent with that near the equator.

[0109] Step 5: Compare the slope and aspect of the Gaussian surface near the equator and the poles in geographic coordinates;

[0110] The estimated and true values ​​of the slope and aspect of the Gaussian surface near the equator and the poles in geographic coordinates and their differences are shown in Figure 5 , Figure 6 As shown, it can be seen that the difference between the slope and aspect of the Gaussian surface DEM calculated by the algorithm in the present invention in geographic coordinates and the true value is very small, wherein near the equator, the root mean square error (RMSE) between the estimated slope and the true value is 0.0332, and the RMSE of the aspect is 0.1989; near the poles, the RMSE of the slope is 0.0220, and the RMSE of the aspect is 0.1399, indicating that the algorithm is very accurate.

[0111] Example 2

[0112] Calculate slope and aspect from real-world SRTM data.

[0113] Determination of experimental area and data: The pixel size of SRTM elevation data is 1″, which is widely used in the field of terrain analysis. This case uses it to verify the geographic coordinate slope aspect algorithm proposed in this invention;

[0114] Considering that elevation change is one of the key factors affecting the accuracy of slope and aspect algorithms under geographic coordinates, three different terrain experimental areas, namely low-undulating hilly areas, medium-undulating hilly gully areas, and steep mountainous areas, were selected to calculate the slope and aspect. Figure 7 As shown, they are located in the Northeast China Plain, the Loess Plateau and the Hengduan Mountains in Southwest China, respectively, and each area covers a spatial range of 1°×1°.

[0115] The method of the present invention is used to calculate the slope and direction of the above three different terrains.

[0116] Comparative Example 1

[0117] For the slope and aspect of the three different terrains in Example 2, the projection transformation method is used to calculate the slope and aspect in geographic coordinates; before the algorithm for directly calculating the slope and aspect in geographic coordinates appears, the projection transformation method is one of the commonly used methods for calculating the slope and aspect in geographic coordinates;

[0118] The process is as follows: first, convert the DEM in geographic coordinates into projection coordinates, then use the traditional slope and aspect algorithm for plane coordinates to calculate the DEM slope and aspect, and finally convert it back to geographic coordinates;

[0119] This method undergoes two projection conversion processes, and the elevation value and spatial position of the original DEM pixel are changed, which will affect the accuracy of the slope and aspect finally obtained. In addition, the two projection conversion processes also seriously reduce the calculation efficiency, especially when a large amount of data needs to be processed. This case uses this method to calculate the slope and aspect data of 1″ pixel size SRTM data in three experimental areas.

[0120] Comparative Example 2

[0121] For the slope and aspect of the three different terrains in Example 2, the fixed pixel size method is used to calculate the slope and aspect in geographic coordinates; the fixed pixel size method is another commonly used method for calculating the slope and aspect in geographic coordinates in the past;

[0122] The process is as follows: first, assume that the pixel size (in length units) of the DEM data in geographic coordinates does not change with the change of latitude, take the size (in length units) of the central pixel of the entire DEM data as the size (in length units) of all pixels in the DEM data, and then use the traditional slope and aspect algorithm for plane coordinates to calculate the DEM slope and aspect;

[0123] The premise of using this method is that the latitude span of the DEM data space range cannot be too large, the pixel size cannot be too large, and the DEM cannot be located in a high-latitude area. Otherwise, the DEM pixel size (in length units) of the pixels at the north and south ends will be too different, which will affect the accuracy of the calculated slope and aspect. In addition, there is currently no ready-made software tool that integrates this method, and users need to implement it through programming.

[0124] This case also uses this method to calculate the slope and aspect data of 1″ pixel size SRTM data in three experimental areas.

[0125] Compare the algorithm in Example 2 of the present invention to calculate the slope and aspect under geographic coordinates with the two commonly used methods in Comparative Example 1 and Comparative Example 2;

[0126] like Figure 8 , Fig. 9 , Fig.10The slope and aspect data of the three experimental areas obtained by the algorithm of the present invention and the other two methods are shown respectively; it can be seen that the macroscopic spatial distribution texture of the slope and aspect data obtained by the three methods in the three experimental areas is similar; in addition, from the statistical data: in the low-undulating hilly area, the RMSE of the slope between the algorithm of the present invention and the fixed pixel size method is 0.0100, and the RMSE of the aspect is 0.1259, and the RMSE of the slope between the algorithm of the present invention and the projection transformation method is 0.5360, and the RMSE of the aspect is 45.4278; in the medium-undulating hilly area In the hilly and gully area, the RMSE of slope between the algorithm of the present invention and the fixed pixel size method is 0.0513, and the RMSE of aspect is 0.1167, while the RMSE of slope between the algorithm of the present invention and the projection conversion method is 1.9558, and the RMSE of aspect is 30.8595; in the steep mountainous area, the RMSE of slope between the algorithm of the present invention and the fixed pixel size method is 0.0641, and the RMSE of aspect is 0.1205, while the RMSE of slope between the algorithm of the present invention and the projection conversion method is 2.3175; the RMSE of aspect is 23.8447;

[0127] In summary, there is a difference between the slope and aspect under geographic coordinates obtained by the traditional projection transformation method and the fixed pixel size method and the slope and aspect obtained by the algorithm in the present invention, wherein the difference between the slopes increases with the increase of terrain undulation, while the difference between the aspects decreases with the increase of terrain undulation; in addition, the difference in the aspects is greater than the difference in the slope; compared with the projection transformation method, the advantage of the algorithm in the present invention is more obvious.

[0128] Example 3

[0129] Application of slope and aspect data in soil erosion calculation in Muli mining area, Qinghai Province.

[0130] Step 1: The Muli mining area in Qinghai Province is a composite erosion area of ​​hydraulic erosion and freeze-thaw erosion. Therefore, the slope and aspect data are respectively introduced into the hydraulic erosion and freeze-thaw erosion formulas to calculate the intensity of the two types of erosion. First, the slope and aspect algorithm derived by the present invention is applied to the DEM data of the high-cold and high-altitude Muli mining area in the geographic coordinate system to obtain the regional slope and aspect raster data, and determine the unique slope and aspect values ​​of different grids. The slope and aspect of the Muli mining area are as follows: Fig.11 and Fig.12 As shown;

[0131] Step 2: Classify the slope data of the Muli mining area, and then calculate the slope factor, with thresholds of 5°, 10° and 30°; in the interval (0, 5°), the slope factor is calculated using the formula S=10.8sinθ+0.03; in the interval [5°, 10°), the slope factor is calculated using the formula S=16.8sinθ-0.05; in the interval ≥10°, the slope factor is calculated using the formula S=21.91sinθ-0.96, where when the slope of the land use (including forest land and grassland) plot is greater than 30°, 30° is used to calculate the slope factor; in addition to the above provisions, forest land and grassland are calculated using the formula S=10.8sinθ+0.03, and the slope factor S value is extracted by classification;

[0132] Step 3: The slope factor S is obtained and the Chinese Soil Loss Equation (CSLE) is used to calculate hydraulic erosion by superposition of the seven-factor method. The slope factor S is multiplied with the rainfall erosivity factor R, the soil erodibility factor K, the slope length factor L, the vegetation coverage and biological measures factor B, the engineering measures factor E, and the tillage measures factor T. Combined with the human-disturbed plots in the area, field investigations are conducted one by one to investigate the erosion type and intensity of the human-disturbed plots. Based on the restoration and management of soil and water loss, the surface coverage, forest and grass plant measures, the implementation area of ​​engineering measures, and the reduction of soil and water loss area or intensity of the human-disturbed plots are analyzed to obtain the hydraulic erosion data of the Muli mining area.

[0133] Step 4: The calculated slope and aspect data are respectively brought into the freeze-thaw erosion calculation, and the six-factor superposition calculation and analysis is carried out with the annual average freeze-thaw cycle days, daily average freeze-thaw phase change water volume, annual average precipitation and vegetation coverage data to obtain the freeze-thaw erosion data of Muli mining area, such as Fig.13 As shown;

[0134] Step 5: Determine the dominant erosion type in the overlapping area of ​​hydraulic erosion and freeze-thaw erosion; when hydraulic erosion and freeze-thaw erosion exist simultaneously on the same plot, the dominant erosion type shall be designated as the erosion type of the plot.

[0135] Comparative Example 3

[0136] The traditional projection transformation method was used to obtain the slope and aspect data of the Muli mining area in Qinghai Province, and to calculate soil erosion.

[0137] Comparative analysis of the accuracy of the algorithm of the present invention and the traditional projection conversion method in the process of calculating erosion intensity shows that: the slope and aspect data directly obtained from the geographic coordinates by the method of the present invention have significantly improved accuracy than those obtained by the traditional projection conversion method. According to calculation and analysis, the average error of the slope under the traditional method is 0.033, and the average error of the aspect is 0.094; the new method for calculating the slope and aspect using geographic coordinates obtains an average error of 0.025 for the slope, and an average error of 0.002 for the aspect; the new method improves the slope accuracy by more than 24.2% compared with the traditional method, and improves the aspect accuracy by more than 97.9% compared with the traditional method; in addition, compared with the projection conversion method, the slope and aspect calculated from geographic coordinates saves two projection conversion steps, and the efficiency is improved by more than 200%. This advantage is particularly prominent in large-area data processing.

[0138] On the basis of the results of multi-factor quantitative evaluation of regional soil erosion, combined with land use types, soil erosion plots, and combined with field investigation and verification, the soil erosion map and its intensity spatial distribution of the Muli mining area were formed; through field investigation and verification results, it was found that the slope and aspect data calculated using the present invention can effectively improve the accuracy of monitoring results.

Claims

1. A slope and aspect calculation method for a geographic coordinate digital elevation model, characterized in that: Please follow the steps below to implement: Step 1: Take a pixel in the DEM data under geographic coordinates as the center and obtain a 3×3 local pixel window; The local pixel window of step 1 is based on the earth ellipsoid, and the surface of the ellipsoid is divided using longitude and latitude. The boundaries of each pixel are located on different meridian circles and latitude circles, respectively. Each pixel represents a local curved surface on the surface of the earth ellipsoid, and the distance between two adjacent pixels is the length of the arc on the surface of the ellipsoid, expressed in the angle unit "degree, minute, second"; Step 2: Project the local pixel window onto the plane, and establish a local plane rectangular coordinate system with the central pixel as the origin; The step 2 is specifically implemented according to the following steps: Step 2.1: Use the Evans surface model to simulate the local terrain surface of 3×3 pixels. The Evans surface formula is shown in (1): Among them, x and y represent the horizontal and vertical coordinates of each pixel in the 3×3 window respectively, z represents the elevation value of the pixel, r is the quadratic term coefficient of x, t is the quadratic term coefficient of y, s is the quadratic term coefficient of xy, p is the linear term coefficient of x, q is the linear term coefficient of y, and f is the constant term. These coefficients and constant terms are used to control the shape of the surface, among which p and q have the same meaning as p and q in the slope and aspect formula; Step 2.2: Project the local surface in step 1 onto the plane coordinates. The specific projection method is as follows: Draw the chords of the arc formed by the outer boundaries of the 8 pixels at the non-central position of the 3×3 local window. These chords intersect at the outer vertices of the pixels on the four diagonals. Since the chords on the latitude circle are parallel, these chords must be located on the same plane. This plane is used as the projection surface of the local surface. Step 2.3, on the plane constructed in step 2.2, establish a local plane coordinate system at the center of the pixel, with the row direction as the x-axis and the column direction as the y-axis; Step 3, calculate the elevation change rate of the central pixel according to the elevation and coordinate values ​​of each pixel in the local window, and substitute them into the slope aspect formula to obtain the slope aspect value of the central pixel; The step 3 is specifically implemented according to the following steps: Step 3.1: The coordinate values ​​of each pixel, that is, their distances to the x and y axes, are represented by w1~w 12 It can be expressed according to the chord length C of the weft loop. par Formula (2) and the chord length C of the warp coil mer Formula (3) calculates: Among them, a and b are the major and minor semi-axes of the earth ellipsoid, ω is the pixel size of the DEM in geographic coordinates, θ is the eccentricity angle between two adjacent pixels on the latitude circle, θ1 and θ2 are the eccentricity angles between two adjacent pixels on the meridian circle, and the eccentricity angle and latitude are The relationship between can be expressed as formula (4): In a 3×3 local pixel window, the latitude difference between two adjacent pixels on the meridian circle is a pixel size ω, so the latitude of the top row of pixels in the window is and the latitude of the bottom row The latitude of the middle row of pixels can be It is expressed as formula (5): In the projection plane rectangular coordinate system, the x-coordinate of each pixel on the left side of the y-axis is negative, and the x-coordinate of each pixel on the right side is positive; the y-coordinate of each pixel above the x-axis is positive, and the y-coordinate of each pixel below is negative; Substituting formula (5) and (4) into formula (2) and (3), we can obtain w1 ~ w 12 Formulas (6) to (10): in, is the latitude of the central pixel of the 3×3 local pixel window; given the DEM data in geographic coordinates, a, b, and ω are all known constants, and It changes with the location of the pixel and is easy to obtain based on the pixel size and spatial range attributes of the DEM data; The elevation and coordinate values ​​of each pixel are expressed as follows: (z1, -w1, w7), (z2, 0, w8), (z3, w2, w9), (z4, -w3, 0), (z5, 0, 0), (z6, w4, 0), (z7, -w5, -w10), (z8, 0, -w11), (z9, w6, w12); Step 3.2: Substitute the coordinates and elevation values ​​of each pixel in step 3.1 into the Evans surface formula, and let w1=w2=g1, w3=w4=g2, w5=w6=g3, w7=w8=w9=f1, w 10 =w 11 =w 12 =f2, and the calculation formulas of p and q in the geographic coordinate system are obtained based on least squares regression as shown in (11) and (12): Step 3.3: Substitute the calculated results of p and q in step 3.2 into the slope calculation formula (13) and the slope aspect formula (14): The slope and aspect value of the central pixel can be obtained; Step 4: Perform the above steps for each pixel of the DEM data under the geographic coordinates in turn to complete the slope and aspect calculation of the entire DEM data.

Images