A LS factor extraction method suitable for large-scale geographic coordinate system raster data

Abstract

The invention relates to an LS factor extraction method applicable to large-scale geographic coordinate system raster data, and is applicable to the design and implementation of the LS factor extraction method for geographic coordinate system raster data, which can improve the LS factor extraction efficiency, and the designed calculation process for large-scale LS factors improves the large-scale LS factor extraction method; the method can quickly and effectively extract and calculate LS factors based on high-resolution global scale SRTM1; the method comprises the following steps: step 1, data block division; step 2, adding a buffer and converting the raster data format into ASCII data; step 3, extracting LS factors with buffer data; step 4, converting the ASCII result data format and extracting each LS factor data without buffer; step 5, extracting the LS factors of the raster of the geographic coordinate system of the required scale by data fusion; the LS factor data obtained in step 4 is fused by a mosaic tool in Arcmap software to obtain the raster LS factors of the geographic coordinate system of the required scale.

Description

Technical Field

[0001] The invention relates to a digital terrain analysis technology intersecting computer science, erosion science and geography, and in particular to an LS factor extraction method suitable for large-scale geographic coordinate system raster data. Background Art

[0002] Topography is an important factor affecting soil erosion. In the US Universal Soil Loss Equation / revised Universal Soil Loss Equation (USLE / REUSLE) and the Chinese Soil Loss Equation (CSLE), the topography (LS) factor is an important indicator for the quantitative calculation of soil loss. At present, the extraction algorithm for regional scale LS factor is mainly based on DEM (Digital Elevation Model) raster data in the projection coordinate system, while there are few reports on the extraction of LS factor based on raster data in geographic coordinate system (such as SRTM, ASTER GDEM, AW3D30 DEM, etc.) in large scale. The data uses raster data with a resolution of 1 arc second (about 30 meters) as the basic unit, adopts the commonly used "WGS_1984" geographic coordinate system, and the file format is hgt, with a resolution of about 30m. Each map of the entire SRTM1 data file covers 1° of longitude and latitude. The global data volume is very large, with a total area of ​​more than 1.19×10 8 km 2 It is necessary to use geographic coordinate system raster data based on SRTM1 to extract large-scale LS factors.

[0003] Most current studies convert the format and projection of SRTM1 and other geographic coordinate system raster data, and then use common GIS tools (such as ArcGIS, ENVI and SAGA software) to gradually complete the calculation of LS factors. For example, Yang Qinke et al. used SRTM1 and 30-meter resolution ASTER GDEM elevation data in the Northeast sample area and the Loess sample area. The SRTM1 data based on the geographic coordinate system uses longitude and latitude as the basic unit. Since the data is projected as spherical distance, the shape of the SRTM1 grid is approximately trapezoidal. How to calculate the terrain factor related parameters (such as slope, slope length, slope length (L) factor, slope (S) factor and terrain (LS) factor, etc.) directly through SRTM1 data without projection conversion, so as to complete the extraction of large-scale or even global LS factors has become a difficult problem in research.

[0004] At present, most LS factor extraction algorithms are based on DEM. They mainly traverse the grid, compare the elevation values ​​of the central grid with the surrounding grids, and determine the corresponding slope, flow direction, unit slope length and catchment area according to the flow direction algorithm idea, so as to determine the L factor and S factor to calculate the LS factor. The flow direction algorithm is divided into single flow direction and multi-flow direction, among which the single flow direction D8 algorithm is the most widely used. Since the DEM grid is a square with equal side lengths and the side lengths are known, it is relatively easy to calculate the terrain parameters. However, the extraction of large-scale LS factors is very time-consuming due to the large amount of data. Some studies have used SRTM to calculate relevant terrain parameters for the European scale soil erosion problem, but the research scope is limited to the European scale; there are also studies on the calculation of LS factors based on the RUSLE model using high-resolution spatial distribution data for global scale soil erosion problems, but the 3 arc second resolution SRTM data used is difficult to guarantee the accuracy of the extraction accuracy. At present, there is also a lack of specific calculation processes and methods for extracting LS factors through large-scale SRTM1 data sets, and the estimation of large-scale or even global LS factors is not fast and efficient.

[0005] How to complete the extraction of LS factors in the geographic coordinate system through the single-flow D8 algorithm and realize the extraction of massive LS factors through large-scale geographic coordinate system raster data is a technical problem to be solved.

[0006] The design and implementation of the LS factor extraction method suitable for raster data in geographic coordinate system effectively avoids the process of coordinate conversion of raster data in geographic coordinate system, which can improve the efficiency of LS factor extraction. In addition, the calculation process for large-scale LS factors is designed to improve the method of large-scale LS factor extraction. This method can quickly and effectively extract and calculate LS factors based on high-resolution global scale SRTM1, and then be applied to global soil erosion mapping and analysis. Summary of the invention

[0007] The purpose of the present invention is to overcome the shortcomings of the prior art and provide an LS factor extraction method suitable for large-scale geographic coordinate system raster data. The LS factor extraction method suitable for geographic coordinate system raster data is designed and implemented, which effectively avoids the process of coordinate conversion of raster data under the geographic coordinate system, can improve the LS factor extraction efficiency, and designs a calculation process for large-scale LS factors, which improves the way of large-scale LS factor extraction; the method can quickly and effectively extract and calculate LS factors based on high-resolution global scale SRTM1, and then be applied to global soil erosion mapping and analysis.

[0008] To achieve the above technical objectives, the technical solution of the present invention is: a LS factor extraction method suitable for large-scale geographic coordinate system raster data, comprising the following steps:

[0009] Step 1: Data segmentation:

[0010] Step 1.1, merge the grid data of the required scale geographic coordinate system into a whole large grid through Arcmap mosaic tool;

[0011] Step 1.2: Divide the large raster data into small blocks according to actual needs using the Arcmap clipping tool rules;

[0012] Step 2: Add a buffer and convert the raster data format:

[0013] Step 2.1, using the mosaic tool in Arcmap software, each single block of geographic coordinate system raster data after division, considering the four directions around the current data block, when raster data exists, add 1° buffer range data for data fusion;

[0014] Step 2.2: Use the raster-to-text tool in Arcmap software to convert the fused raster data into ASCII data format;

[0015] Step 3: LS factor extraction with buffer data:

[0016] Step 3.1, create a log file;

[0017] Step 3.2, read the ASCII data header file and the parameter information in the elevation and LS factor extraction;

[0018] Step 3.3, fill the valueless points and depressions in the geographic coordinate system raster data and update the elevation values;

[0019] Step 3.4, traverse the elevation two-dimensional array to calculate the slope, flow direction and unit slope length; apply for a matrix space to store the slope, flow direction and unit slope length. Each time the elevation matrix is ​​traversed, first determine whether the current grid is a valueless point; if it is a valueless point, directly record the slope as "0" and skip the point to proceed to the next point judgment; if it is not a valueless point, according to the D8 flow algorithm idea, calculate the maximum slope value according to formula (1) as the slope of the grid and record it in the slope matrix: where E c Represents the elevation value of the center grid, E i Represents the elevation value of the current grid, cellsize is the distance between the current grid and the center grid, and the north-south distance of the grid is recorded as h x , the east-west distance of the grid is recorded as h y , the diagonal distance of the grid is recorded as diagcellsize, h x and h y From equations (2) and (3), we can conclude that diagcellsize is based on h x and h yCalculated by the Pythagorean theorem, where θ is the north-south width of the grid pixel in the geographic coordinate system; set the slope in the case of flat land and depression to 0.1, repeat the above steps until the slope calculation of all value points is completed; set the direction of the maximum slope as the flow direction of the grid, record the value in the flow direction matrix according to the corresponding direction flow code, repeat the above steps until the flow direction calculation of all value points is completed; record the celsize value as the unit slope length value in the unit slope length matrix;

[0020] slope=max(deg·arctan((E c -E i ) / cellsize)) (1)

[0021] hx=30.8874791 (2)

[0022] hy=30.8874791·cosθ (3)

[0023] Step 3.5, traverse the corresponding two-dimensional array, calculate the initial catchment area, initial slope length and catchment area:

[0024] Step 1: Calculate the initial catchment area and apply for a matrix space to store the initial catchment area. Each time you traverse the elevation matrix and flow matrix, first determine whether the current grid is a valueless point; if there is no value, proceed to the next grid calculation; if there is a value, initialize the catchment area based on the flow matrix and the cumulative number of grid records; the area of ​​the grid is equal to the product of the length and width, that is, h x ·h y , this area is the initial catchment area value; repeat the above steps until the initial catchment area of ​​all valued points is assigned; Step 2: Calculate the initial slope length, apply for the matrix space to store the initial slope length, each time traversing the elevation matrix and the flow matrix, first determine whether the current grid is a valueless point; if there is no value, proceed to the next grid calculation, if there is a value, initialize the slope length according to the flow matrix and the unit slope length matrix record; the unit slope length value is the initial slope length value; repeat the above steps until the initial slope length of all valued points is assigned; Step 3: Calculate the catchment area, apply for the matrix space to store the catchment area, traverse the flow matrix and the initial catchment area matrix, The sum of the catchment area values ​​flowing to the current grid is recorded as amount; compare amount with the catchment area value of the current grid, and select the larger value as the catchment area value of the current grid; traverse the entire initialized catchment area array forward to calculate the catchment area value of the entire grid data; traverse the entire initialized catchment area array backward to calculate the catchment area value of the entire grid data; if the value of amount is not assigned to the catchment area value of the current grid during the forward and reverse processes, it means that the catchment area extraction is completed and the loop ends, otherwise, the calculation is repeated from the beginning; repeat the above steps until the catchment area assignment of all points is completed;

[0025] Step 3.6, traverse the corresponding two-dimensional array, and calculate the cumulative slope length and the entrance and exit slope length according to the slope cutoff and channel cutoff:

[0026] Step 1: Set slope cutoff and channel cutoff; apply for matrix space to store cutoff values; each time the elevation matrix is ​​traversed, determine whether the current grid is a valueless point, and consider the following cutoff situations:

[0027] If there is no value, the point is set to truncation, otherwise it is set to non-truncation; traverse the slope matrix, take 5% of the slope as the dividing point, less than 5%, the truncation factor is set to 0.7; when it is greater than or equal to 5%, the truncation factor value is set to 0.5; when the product of the grid slope and the truncation factor is greater than the grid slope in the outflow direction, the grid is set to truncation;

[0028] Traverse the catchment area matrix to determine whether the catchment area value of the grid is greater than the set river network threshold. If so, set the grid to be truncated, otherwise, not set to be truncated;

[0029] Repeat the above steps to set the cutoff for each grid;

[0030] Step 2: Calculate the cumulative slope length, apply for the matrix space to store the cumulative slope length, traverse the truncation matrix and the initial slope length matrix, first determine whether the current grid is truncated, if it is truncated, the slope length of the current grid is equal to half of the initial slope length, if it is not truncated, the initial slope length of the current grid remains unchanged; repeat the above steps to set the calculation of each grid to add the initial slope length after truncation to the initial slope length array; then declare that the initial value of the temporary variable total is set to 0; traverse the flow direction matrix, assuming that grid a is a grid adjacent to the current grid c, and the flow direction of grid a points to grid c, if grid a is truncated, total plus half of the slope length of grid a, if not If it is truncated, total is added to the slope length value of grid a; the sum of the slope length values ​​flowing to the current grid is calculated by this method and recorded as total; total is compared with the slope length value of the current grid, and the larger value is selected as the slope length value of the current grid; the entire initialized slope length array is traversed forward to calculate the slope length value of the entire grid data; the entire initialized slope length array is traversed backward to calculate the slope length value of the entire grid data; if the operation of assigning the value of total to the slope length value of the current grid does not occur in the forward and reverse processes, it means that the cumulative slope length extraction is completed and the loop ends, otherwise the calculation is repeated from the beginning; the above steps are repeated until the cumulative slope length assignment of all points is completed;

[0031] Step 3: Calculate the entrance and exit slope lengths, apply for matrix space to store the exit and entrance slope lengths, and each time you traverse the elevation matrix, first determine whether the current grid is a valueless point; if there is no value, proceed to the next grid for calculation; if there is a value, calculate the entrance and exit slope lengths of the current grid based on the flow direction matrix, initial slope length matrix, and cumulative slope length matrix; repeat the above steps until all value points are found and the entrance slope length assignment is completed;

[0032] Step 3.7, traverse the corresponding two-dimensional array and calculate the S factor and L factor:

[0033] Step 1: Calculate the S factor and apply for a matrix space to store the slope factor (S); each time the elevation matrix is ​​traversed, first determine whether the current grid is a valueless point; if there is no value, proceed to the next grid for calculation; if there is a value, calculate the S factor value according to the slope matrix and formula (4); repeat the above steps to complete the S factor calculation for each grid; where θ is the slope;

[0034]

[0035] Step 2: Calculate the L factor and apply for a matrix space to store the slope length factor (L); each time the elevation matrix is ​​traversed, first determine whether the current grid is a valueless point; if there is no value, proceed to the next grid for calculation; if there is a value, calculate the segmented slope L value according to the entrance and exit slope length matrix according to formulas (5) and (6); when the entrance slope length is less than the exit slope length value, use formula (6); otherwise use formula (6); repeat the above steps to complete the L factor calculation of each grid; where λ is the slope length, m is the slope length index, and λout and λin are the slope lengths of the grid exit and entrance, respectively (m);

[0036] L = (λ / 22.13) m (5)

[0037]

[0038] Step 3.8, traverse the corresponding two-dimensional array and calculate the LS factor: apply for a matrix space to store the slope length factor (LS); each time the elevation matrix is ​​traversed, first determine whether the current grid is a valueless point; if there is no value, proceed to the next grid for calculation; if there is a value, calculate the LS factor value according to the L factor and S factor matrix according to formula (7); repeat the above steps to complete the LS factor calculation for each grid;

[0039] LS=L·S (7)

[0040] Step 4: Convert the result data format into ASCII and extract the LS factor data without buffer:

[0041] Step 4.1, convert the LS factor data of each block with buffer back to raster data through the ASCII to raster tool in Arcmap software;

[0042] Step 4.2, using the clipping tool in Arcmap software, clip the LS factor raster data obtained above using each block raster data without buffer;

[0043] Step 5: Data fusion to extract the LS factor of the geographic coordinate system grid of the required scale; the LS data obtained in step 4 can be fused through the mosaic tool in Arcmap software to obtain the LS factor of the grid in the geographic coordinate system of the required scale.

[0044] Preferably, in step 1, a block method is used to process data in a large scale range.

[0045] Preferably, before step 3, add 1° buffer range data and convert the raster format to ASCII text; this method uses Arcmap software to complete buffer addition and raster data conversion through mosaicking and raster to ASCII conversion.

[0046] Preferably, in step 3.2, the process of reading the ASCII file of the geographic coordinate system raster data and the parameter information in the LS factor extraction is:

[0047] Step 1: Create a structure named DemData to store the ASCII header information and the parameter information set in the LS factor extraction, and apply for a two-dimensional array to save the elevation value;

[0048] Step 2: Open the raster data text file. If the opening fails, write the log and stop the execution.

[0049] Step 3: First read the contents of the ASCII file line by line, and record them in the format of "name-space-value" in the file header; then store each line of data read into each string, and then split the string with spaces, convert the obtained value into the type of the value and save it to the corresponding attribute of the created structure, and repeat this process until the header file is read; then record the parameter information such as flow direction coding, truncation factor and river network threshold set in the LS factor extraction in the same form in the corresponding attribute of the structure.

[0050] Preferably, in step 3.4, the grid encoding method of the D8 flow direction algorithm takes the central grid as an example, and determines the flow directions in eight directions around the central grid according to the elevation values ​​corresponding to different grids. The directions from east, southeast, south, southwest, west, northwest, north to northeast are recorded as 1, 2, 4, 8, 16, 32, 64 and 128 respectively.

[0051] Preferably, in Step 2 of Step 3.6, the initial slope length value is updated by re-assigning the unit slope length according to whether the grid is truncated: for the base grid, the initial slope length is the original unit slope length value; for the truncated grid, the initial slope length is half of the original unit slope length value.

[0052] Preferably, before step 5, the calculated ASCII result data is converted back to raster data and the buffer is removed. This method uses Arcmap software to complete the raster data conversion and the construction of the LS factor of the non-overlapping area through ASCII to raster conversion and clipping.

[0053] It can be seen from the above description that the beneficial effects of the present invention are:

[0054] 1. According to the characteristics of geographic coordinate system grid data, a geographic coordinate system grid model based on SRTM1 is established. Based on the idea of ​​single flow direction D8 algorithm, the LS factor is calculated and extracted step by step according to the grid characteristics. The LS factor extracted by the algorithm is in line with the actual research and has high extraction efficiency.

[0055] 2. In view of the large amount of large-scale raster data, a large-scale LS factor extraction process is designed. Based on the idea of ​​"block calculation and then merging", a buffer is set when calculating each block of data to ensure the accuracy of the calculation. The extraction of LS factors in a larger range is completed through ARCGIS and other geographic information systems. The process is operational and provides a solution for the extraction of large-scale and even global LS factors;

[0056] 3. According to the final results of the method, it can be seen that: the LS factor results calculated by this method and the LS factor results obtained by DEM transformation data under the projection coordinate system are unified in the same coordinate system, and 99% of the LS differences are concentrated between ±1; the global LS factor extraction results completed by this algorithm and process show that the global LS factor point map has a smooth and seamless transition, which conforms to the value range of LS factor calculation, and the calculation results of the European region are similar to the existing European regional LS factor extraction research results in distribution and the terrain features are more obvious. The calculation results of this method are reasonable and highly feasible. BRIEF DESCRIPTION OF THE DRAWINGS

[0057] Figure 1 This is a flow chart of the LS factor extraction method applicable to large-scale geographic coordinate system grids (taking SRTM1 as an example);

[0058] Figure 2 This is a flow chart of the LS factor extraction algorithm based on the geographic coordinate system grid (taking SRTM1 as an example);

[0059] Figure 3 It is a SRTM1 grid data longitude and latitude surface map

[0060] Figure 4SRTM1 grid flow direction and encoding method diagram

[0061] Figure 5 This is the LS factor result map of Xiannangou based on DEM data

[0062] Figure 6 This is the LS factor result diagram of Xiannangou based on SRTM1 data

[0063] Figure 7 It is the frequency statistics chart of the LS factor difference between SRTM1 and DEM. DETAILED DESCRIPTION

[0064] The present invention will be further described below in conjunction with specific embodiments.

[0065] Soil erosion has been identified as one of the major soil threats. In order to cope with the complex challenges of soil erosion, it is necessary to conduct a global soil erosion assessment and produce a global soil erosion map to address the risks brought by soil erosion. The topographic (LS) factor is an important factor in commonly used soil erosion models (USLE / RUSLE / CSLE). The present invention designs and implements an LS factor extraction method suitable for large-scale geographic coordinate system raster data, providing technical support and solutions for calculating the global LS factor.

[0066] The LS factor extraction method based on GIS requires first transforming the coordinates of the grid in the geographic coordinate system, and then using common hydrological analysis tools such as slope, flow direction, and grid calculator to gradually extract the LS factor. In addition, when extracting large-scale LS factors, the increase in the amount of data will make the operation very complicated, requiring a high degree of professionalism and rigor. In view of the defects or deficiencies in the above-mentioned existing research technologies, the design and implementation of the LS factor extraction method and large-scale LS factor extraction process suitable for geographic coordinate system raster data will greatly facilitate user use and provide technical support and solutions for the extraction of large-scale and even global LS factors.

[0067] See also Figure 1 and Figure 2 Taking the raster data in the global geographic coordinate system (SRTM1 as the main data, and 30-meter resolution ASTERGDEM as the fusion data of the hole area as the supplementary data) as an example, the specific process of the LS factor extraction method applicable to the large-scale geographic coordinate system raster data is divided into the following steps:

[0068] Step 1: Divide the geographic coordinate system raster data calculated at the global scale into blocks;

[0069] Using the clipping tool in Arcmap software, the large blocks of SRTM1 and ASTER GDEM fusion data that are spliced ​​globally are divided according to actual needs. The structure of each data block after division is consistent, all of which are grid rasters. The global grid data is divided into 228 large blocks.

[0070] Step 2: Add a 1° buffer for each chunk of data and convert the chunk into ASCII text;

[0071] Use the fusion tool in Arcmap software to add each single block of divided raster data, and add 1° buffer range data in four directions around the raster data (when raster data exists) to reduce the data boundary effect. In this way, the data is fused into larger blocks of raster data; then use the raster to text tool in Arcmap software to convert the fused raster data blocks into ASCII data. Each block of data contains header information and elevation values.

[0072] Step 3: Complete the LS factor extraction of each block of raster data with buffer based on the LS factor extraction algorithm;

[0073] (1) Create a log file

[0074] Considering that problems may occur during algorithm execution, especially when the amount of raster data is large and the algorithm runs for a long time, a log file is created at the output location when the algorithm is executed and a log file class is created to output the necessary log information.

[0075] (2) Reading data

[0076] i. Input raster fusion data with buffer and read parameter information in file header and LS factor extraction: first apply to create a structure to store the header file information of geographic coordinate system raster data; then open the raster data text file. If the opening fails, write to the log and stop execution; then read the content in the file line by line, and the format of the file header is recorded in the form of "name-space-value"; then store each line of data read into each string, and then split the string with spaces, convert the obtained value into the type of the value and save it to the corresponding attribute of the created structure, and repeat the process until the header file is read; finally, the parameter information such as flow direction coding, truncation factor and river network threshold set in LS factor extraction is recorded in the corresponding attribute of the structure in the same form.

[0077] ii. Read the data elevation value: Apply to create a two-dimensional array of elevation. Since each data in each row of elevation data is separated by spaces. During the algorithm execution, the elevation value is read row by row, and each row of data is read in the form of a string and separated by spaces. Then, each string of separated data is converted into a float type and stored in the elevation data matrix. For the convenience of calculation, the elevation matrix is ​​expanded twice, so reading the data will place the data in a matrix that is expanded twice.

[0078] (3) Filling of points with no value and depressions.

[0079] i. Fill the valueless points and depressions in the elevation data: Valueless points are erroneous values ​​in the elevation data due to errors in the measurement method itself. In terrain data, depressions are manifested as the phenomenon that the surrounding grid data are higher than the central grid data. When the algorithm is executed, these grids need to be filled with the minimum elevation value of the surrounding eight grids.

[0080] (4) Traverse the two-dimensional elevation array to calculate the slope, flow direction and unit slope length.

[0081] i. Apply for matrix space to store slope, flow direction and unit slope length. Each time you traverse the elevation matrix, first determine whether the current grid is a valueless point. If it is a valueless point, directly record the slope as "0" and skip the point to proceed to the next point judgment; if it is not a valueless point, calculate the grid cell side length cellsize as follows:

[0082] The result obtained from formula (2) and (3) is used to set the grid width (distance) h in the longitude and latitude directions of the geographic coordinate system grid data. x With h y , for h y There is formula (8):

[0083] cosθ=cos((yllcorner+((nrows-1)-(i-2))·cellsize) / deg (8)

[0084] Among them, yllcorner is the coordinate information of the space where the data is located obtained from the grid header file information; nrows is the number of rows after two circles are expanded based on the original number of rows of the grid data; i represents the offset of the current traversed grid relative to the starting grid in the longitude direction. For the pixel distance diagcellsize in the diagonal direction of the grid, it is calculated according to the Pythagorean theorem using formula (9):

[0085]

[0086] ii. Calculate the slope, flow direction and unit slope length according to the D8 flow direction algorithm:

[0087] ①: Compare the elevation values ​​of the current grid with those of the surrounding eight grids in turn, and calculate the angle between the surrounding grids and the central grid. The calculation method is shown in formula (10):

[0088] angle = deg arctan ((E c -E i ) / cellsize) (10)

[0089] Where E c Represents the elevation value of the center grid, E i Represents the elevation value of the current grid. Cellsize is the distance between the current grid and the center grid, where the east-west direction of the grid is h y , the north-south direction of the grid is h x , the diagonal direction of the grid is diagcellsize.

[0090] ②: Select the maximum value of the surrounding angles as the slope value of the central grid. The calculation method is shown in formula (11):

[0091] slope=max(deg·arctan((E c -E i ) / cellsize)) (11)

[0092] ③: Select the direction with the maximum surrounding angle as the flow direction of the central grid, and assign a value to it according to the flow direction code.

[0093] ④: According to the grid flow direction, the grid cell side length value celsize is assigned as the cell slope length value.

[0094] iii. Setting the slope in the case of flat land and depression: Since flat land and depression have no flow direction, this algorithm simplifies this process. To ensure the connection between the grid and the river network, the slope is set to the minimum value of 0.1. Repeat the above steps until all grids are traversed, and output the calculated slope results to the slope array. (5) Traverse the corresponding two-dimensional array and calculate the initial catchment area, initial slope length and catchment area.

[0095] i. Extraction of initial catchment area: Apply for matrix space to store the initial catchment area. Each time you traverse the elevation matrix and flow matrix, first determine whether the current grid is a valueless point. If there is no value, proceed to the next grid calculation. If there is a value, initialize the catchment area based on the flow matrix and the cumulative number of grid records. The area of ​​the grid is equal to the product of the length and width, that is, h x ·h y, this area, i.e. the initial catchment area value, can be obtained by formula (2)(3)(8). Repeat the above steps until the initial catchment area values ​​of all valued points are assigned, and the calculated initial catchment area values ​​are output to the initial catchment area array;

[0096] ii. Extraction of initial slope length: Apply for matrix space to store initial slope length. Each time you traverse the elevation matrix and flow matrix, first determine whether the current grid is a valueless point. If there is no value, proceed to the next grid calculation. If there is a value, initialize the slope length according to the flow matrix and unit slope length matrix records. Assign the unit slope length value to the initial slope length. Repeat the above steps until the initial slope length assignment of all value points is completed, and output the calculated initial slope length to the initial slope length array.

[0097] iii. Extraction of catchment area: Apply for matrix space to store catchment area, traverse the flow direction matrix and the initial catchment area matrix, calculate the sum of the catchment area values ​​flowing to the current grid and record it as amount; compare amount with the catchment area value of the current grid, and select the larger value as the catchment area value of the current grid; traverse the entire initialized catchment area array forward to calculate the catchment area value of the entire grid data; traverse the entire initialized catchment area array backward to calculate the catchment area value of the entire grid data; if the operation of assigning the value of amount to the catchment area value of the current grid does not occur in the forward and reverse processes, it means that the catchment area extraction is completed and the loop ends, otherwise the calculation is repeated from the beginning. Repeat the above steps until the catchment area assignment of all points is completed, and finally output the calculated catchment area to the catchment area array.

[0098] (6) Traverse the corresponding two-dimensional array and calculate the cumulative slope length and the entrance and exit slope length according to the slope cutoff and channel cutoff.

[0099] i. Apply for matrix space to store truncation values. Each time the elevation matrix is ​​traversed, determine whether the current grid is a valueless point. If there is no value, set the point to truncation, otherwise set it to non-truncation; traverse the slope matrix to determine whether the slope of the grid is less than 5% (about 2.861°). Since the truncation factor is 0.7 at this time, when the product of the grid slope and the truncation factor is greater than the grid slope in the outflow direction, the grid is set to truncation; if the slope is greater than or equal to 5% (about 2.861°), the truncation factor is 0.5 at this time. When the product of the grid slope and the truncation factor is greater than the grid slope in the outflow direction, the grid is also set to truncation; then traverse the catchment area matrix to determine whether the catchment area value of the grid is greater than the set river network threshold. If so, set the grid to truncation, otherwise not set to truncation; repeat the above steps to set the truncation of each grid and output the calculated truncation setting to the truncation array.

[0100] ii. Cumulative slope length extraction: Apply for matrix space to store cumulative slope length, traverse the truncation matrix and the initial slope length matrix, first determine whether the current grid is truncated, if it is truncated, the slope length of the current grid is equal to half of the initial slope length, if it is not truncated, the initial slope length of the current grid remains unchanged; repeat the above steps to set the calculation of each grid and add the initial slope length after truncation to the initial slope length array.

[0101] Again, declare that the initial value of the temporary variable total is set to 0; traverse the flow direction matrix, assuming that grid a is a grid adjacent to the current grid c, and the flow direction of grid a points to grid c. If grid a is truncated, total plus half of the slope length of grid a, if it is not truncated, total plus the slope length of grid a. In this way, the sum of the slope length values ​​flowing to the current grid is calculated and recorded as total; compare total with the slope length value of the current grid, and select the larger value as the slope length value of the current grid; traverse the entire initialized slope length array forward to calculate the slope length value of the entire grid data; traverse the entire initialized slope length array backward to calculate the slope length value of the entire grid data; if the value of total is not assigned to the slope length value of the current grid during the forward and reverse processes, it means that the cumulative slope length extraction is completed and the loop ends, otherwise the calculation is repeated from the beginning. Repeat the above steps until the cumulative slope length assignment of all points is completed, and finally output the calculated cumulative slope length to the cumulative slope length array;

[0102] iii. Extraction of entrance and exit slope lengths: Apply for matrix space to store exit and entrance slope lengths. Each time the elevation matrix is ​​traversed, first determine whether the current grid is a valueless point. If there is no value, proceed to the next grid for calculation. If there is a value, calculate the entrance and exit slope lengths of the current grid based on the flow direction matrix, initial slope length matrix, and cumulative slope length matrix; repeat the above steps until all exit and entrance slope lengths with values ​​are assigned, and output the calculated exit and entrance slope lengths to the exit and entrance slope length arrays.

[0103] (7) Traverse the corresponding two-dimensional array and calculate the S factor and L factor.

[0104] iS factor extraction: Apply for a matrix space to store the slope factor (S). Each time you traverse the elevation matrix, first determine whether the current grid is a valueless point. If there is no value, proceed to the next grid for calculation. If there is a value, calculate the S factor value based on the slope matrix and formula (4); repeat the above steps to complete the S factor calculation for each grid;

[0105] ii. L factor extraction: Apply for a matrix space to store the slope length factor (L). Each time the elevation matrix is ​​traversed, first determine whether the current grid is a valueless point. If there is no value, proceed to the next grid for calculation. If there is a value, calculate the L value of the segmented slope according to the entrance and exit slope length matrix according to formulas (5) and (6). When the entrance slope length is less than the exit slope length value, use formula (6), otherwise use formula (5); repeat the above steps to complete the L factor calculation of each grid.

[0106] (8) Traverse the corresponding two-dimensional array and calculate the LS factor.

[0107] Apply for the matrix space to store the slope length factor (LS). Each time you traverse the elevation matrix, first determine whether the current grid is a valueless point. If there is no value, proceed to the next grid for calculation. If there is a value, traverse the L factor and S factor matrices and calculate the LS factor value according to formula (7); repeat the above steps to complete the LS factor extraction of each grid.

[0108] Step 4: Convert the LS factor data in ASCII format obtained by the LS factor extraction algorithm into raster data, and use the terrain data of each block without buffer to clip the obtained LS factor data blocks, and extract the data of each LS factor block without buffer;

[0109] The calculation results of each local block completed by the LS factor extraction algorithm of the present invention are ASCII data, and each block of LS factor data is converted into raster data by using the ASCII to raster tool in the Arcmap software. Then, the clipping tool in the Arcmap software is used to clip the LS factor raster data blocks obtained above through each block raster data, thereby obtaining the local block LS factor raster data without buffer (no overlap between blocks).

[0110] Step 5: Fuse the LS factor grid data obtained in step 4 to finally obtain the global LS factor with a resolution of 1 arc second.

[0111] Use the fusion tool in Arcmap software to fuse the local block LS factors obtained in step 4 to obtain seamless LS factor raster data with a global resolution of 1 arc second.

[0112] Experimental part:

[0113] The LS factor extraction algorithm of the geographic coordinate system raster data is compared with the LS factor extraction algorithm of the traditional projection coordinate system raster data, and the global LS factor extraction is completed using the calculation process and algorithm, and compared with the existing European regional research results:

[0114] Experimental background:

[0115] The most direct and effective way to verify the correctness and feasibility of the calculation results obtained from raster data in the geographic coordinate system is to compare and analyze the calculation results with the LS factor extraction algorithm in the traditional projection coordinate system, and to obtain large-scale or even global LS factor results by implementing this method, and to analyze their rationality by referring to existing research.

[0116] Experimental area:

[0117] The SRTM1 and 30-meter resolution DEM data of the Nangou Basin in Suide County, northern Shaanxi, and the global geographic coordinate system raster data (SRTM1 is the main data, and 30-meter resolution ASTERGDEM is the fusion data of the supplementary data in the hole area).

[0118] Experimental methods:

[0119] 1. Project the 1 arc second resolution Xiannangou SRTM1 data obtained by the above LS factor extraction algorithm, and place them in the same coordinate system with the 30 m resolution Xiannangou DEM data calculated by the LS factor extraction algorithm in the traditional projection coordinate system to compare the differences between the two results. Figure 5 It is the LS value of Xiannangou based on the LS factor extraction algorithm in the projection coordinate system of DEM. Figure 6 It is the LS factor value of Xiannangou based on the LS factor extraction algorithm in the geographic coordinate system of SRTM1. Figure 7 This is the frequency statistics of the difference of LS factors of Xiannangou obtained by expanding the difference of LS factors of SRTM1 and DEM by 100 times;

[0120] 2. Use the above process and algorithm to calculate the global geographic coordinate system raster data to complete the extraction of the global LS factor with a resolution of 1 arc second, check whether its distribution and value range are reasonable, and compare the existing European regional LS factor data to see the distribution of the results.

[0121] Result analysis:

[0122] DEM data is a 30-meter digital elevation model digital map, and SRTM1 is a remote sensing elevation map with a 1 arc second resolution (i.e. 30-meter resolution SRTM data). Figure 4 and Figure 5 It can be obtained that the overall spatial distribution map of the LS factor based on the two data is very similar to the range of the LS factor value, and the trends of the LS factors extracted from the two different data maps are consistent. Since the LS factor is calculated from the S factor and the L factor, where the S factor is calculated from the slope and the L factor is calculated from the slope length, and the LS factor is mainly affected by the slope, the calculated value of the LS factor will be affected to a certain extent.

[0123] The STRM1 data value is the span of the grid in longitude. The span of the SRTM data in 1 arc second calculated by the radius of the earth and the geometric formula is approximately equal to 30 meters. Since the span length of the SRTM1 data grid in longitude is fixed, its span length in latitude decreases with the increase of latitude, so the grid size of the SRTM1 data not at the equator is slightly smaller than the DEM grid with a resolution of 30m. The LS factor results show that since the elevation difference between the DEM data and the STRM1 data at the same location is equal, when the flow direction of the grid is not north-south, the grid spacing of SRTM1 will be smaller than that of DEM, so the LS factor value calculated by SRTM1 will be greater than the LS factor value calculated by DEM; when the grid flow direction is north-south, the grid spacing of SRTM1 will be slightly larger than that of DEM, so the LS factor value extracted by STRM1 is sometimes greater than the LS factor value extracted by DEM. From the LS factor difference frequency statistics, we can see that 99% of the differences are concentrated between ±1, and the results are highly consistent with the results of the traditional DEM method. From the global LS factor extraction result map, we can see that the LS factor point map transitions smoothly and seamlessly, conforms to the value range of LS factor calculation, and the value distribution is reasonable. Compared with the existing European regional extraction results, the terrain distribution characteristics of the LS factor results of this method are more obvious, and the stratification phenomenon is more prominent in areas with larger and smaller LS factor values, making it easy to find local LS factor high and low points. Although the LS factor value distribution in some areas is biased, this is based on different LS factor calculation formulas (the actual slope and its corresponding LS factor calculation formula have different parameter settings), and the error is within a reasonable range.

[0124] The LS factor results calculated by SRTM1 in the Nangou area of ​​the county are slightly different from those calculated by DEM, but the results produced by the two data are highly similar, which verifies the accuracy of this LS factor extraction algorithm. The global LS factor data obtained based on the global geographic coordinate system raster data is relatively "smooth", with a correct range of values ​​and a reasonable spatial distribution. It is similar to the existing European regional LS factor results and has more obvious terrain features. The resulting 1 arc second resolution global seamless map can effectively solve the shortcoming of low resolution of the global LS factor point map, and can also better reflect the impact of terrain on slope erosion. This method can quickly, efficiently and accurately extract large-scale and even global LS factors.

[0125] The above are preferred implementation modes of the present invention. Those skilled in the art to which the present invention belongs can also change and modify the above implementation modes. Therefore, the present invention is not limited to the above specific implementation modes. Any obvious improvements, substitutions or modifications made by those skilled in the art on the basis of the present invention belong to the protection scope of the present invention.

Claims

1. A LS factor extraction method suitable for large-scale geographic coordinate system raster data, It is characterized in that The steps include: Step 1: Data segmentation: Step 1.1, merge the grid data of the required scale geographic coordinate system into a whole large grid through Arcmap mosaic tool; Step 1.2: Divide the large raster data into small blocks according to actual needs using the Arcmap clipping tool rules; Step 2: Add a buffer and convert the raster data format: Step 2.1, using the mosaic tool in Arcmap software, each single block of geographic coordinate system raster data after division, considering the four directions around the current data block, when raster data exists, add 1° buffer range data for data fusion; Step 2.2: Use the raster-to-text tool in Arcmap software to convert the fused raster data into ASCII data format; Step 3: LS factor extraction with buffer data: Step 3.1, create a log file; Step 3.2, read the ASCII data header file and the parameter information in the elevation and LS factor extraction; Step 3.3, fill the valueless points and depressions in the geographic coordinate system raster data and update the elevation values; Step 3.4, traverse the elevation two-dimensional array to calculate the slope, flow direction and unit slope length; apply for a matrix space to store the slope, flow direction and unit slope length. Each time the elevation matrix is ​​traversed, first determine whether the current grid is a valueless point; if it is a valueless point, directly record the slope as "0" and skip the point to proceed to the next point judgment; if it is not a valueless point, according to the D8 flow algorithm idea, calculate the maximum slope value according to formula (1) as the slope of the grid and record it in the slope matrix: where E c Represents the elevation value of the center grid, E i Represents the elevation value of the current grid, cellsize is the distance between the current grid and the center grid, and the north-south distance of the grid is recorded as h x , the east-west distance of the grid is recorded as h y , the diagonal distance of the grid is recorded as diagcellsize, h x and h y From equations (2) and (3), we can conclude that diagcellsize is based on h x and h y Calculated by the Pythagorean theorem, where θ is the north-south width of the grid pixel in the geographic coordinate system; set the slope in the case of flat land and depression to 0.1, repeat the above steps until the slope calculation of all value points is completed; set the direction of the maximum slope as the flow direction of the grid, record the value in the flow direction matrix according to the corresponding direction flow direction code, repeat the above steps until the flow direction calculation of all value points is completed; record the cellsize value as the unit slope length value in the unit slope length matrix; slope=max(deg·arctan((E c -E i ) / cellsize)) (1) hx=30.8874791 (2) hy=30.8874791·cosθ (3) Step 3.5, traverse the corresponding two-dimensional array, calculate the initial catchment area, initial slope length and catchment area: Step 1: Calculate the initial catchment area and apply for a matrix space to store the initial catchment area. Each time you traverse the elevation matrix and flow matrix, first determine whether the current grid is a valueless point; if there is no value, proceed to the next grid calculation; if there is a value, initialize the catchment area based on the flow matrix and the cumulative number of grid records; the area of ​​the grid is equal to the product of the length and width, that is, h x ·h y , this area is the initial catchment area value; repeat the above steps until the initial catchment area of ​​all value points is assigned; Step 2: Calculate the initial slope length, apply for the matrix space to store the initial slope length, and each time the elevation matrix and flow direction matrix are traversed, first determine whether the current grid is a valueless point; If there is no value, proceed to the next grid calculation. If there is a value, initialize the slope length according to the flow direction matrix and the unit slope length matrix records; the unit slope length value is the initial slope length value; repeat the above steps until the initial slope length assignment of all valued points is completed; Step 3: Calculate the catchment area, apply for the matrix space to store the catchment area, traverse the flow direction matrix and the initial catchment area matrix, and record the sum of the catchment area values ​​flowing to the current grid as amount; compare amount with the catchment area value of the current grid, and select the larger value as the catchment area value of the current grid; traverse the entire initialized catchment area array forward to calculate the catchment area value of the entire grid data; traverse the entire initialized catchment area array backward to calculate the catchment area value of the entire grid data; if the operation of assigning the value of amount to the catchment area value of the current grid does not occur in the forward and reverse processes, it means that the catchment area extraction is completed and the loop ends, otherwise the calculation is repeated from the beginning; repeat the above steps until the catchment area assignment of all points is completed; Step 3.6, traverse the corresponding two-dimensional array, and calculate the cumulative slope length and the entrance and exit slope length according to the slope cutoff and channel cutoff: Step 1: Set slope cutoff and channel cutoff; apply for matrix space to store cutoff values; each time the elevation matrix is ​​traversed, determine whether the current grid is a valueless point, and consider the following cutoff situations: If there is no value, the point is set to truncation, otherwise it is set to non-truncation; traverse the slope matrix, take 5% of the slope as the dividing point, less than 5%, the truncation factor is set to 0.7; when it is greater than or equal to 5%, the truncation factor value is set to 0.5; when the product of the grid slope and the truncation factor is greater than the grid slope in the outflow direction, the grid is set to truncation; Traverse the catchment area matrix to determine whether the catchment area value of the grid is greater than the set river network threshold. If so, set the grid to be truncated, otherwise, not set to be truncated; Repeat the above steps to set the cutoff for each grid; Step 2: Calculate the cumulative slope length, apply for the matrix space to store the cumulative slope length, traverse the truncation matrix and the initial slope length matrix, first determine whether the current grid is truncated, if it is truncated, the slope length of the current grid is equal to half of the initial slope length, if it is not truncated, the initial slope length of the current grid remains unchanged; Repeat the above steps to set the calculation of the initial slope length after truncation of each grid and add it to the initial slope length array; then declare the initial value of the temporary variable total to be 0; traverse the flow direction matrix, assuming that grid a is a grid adjacent to the current grid c, and the flow direction of grid a points to grid c, if grid a is truncated, total plus half of the slope length of grid a, if not truncated, total plus the slope length value of grid a; use this method to calculate the sum of the slope length values ​​flowing to the current grid and record it as total; compare total with the slope length value of the current grid, and select the larger value as the slope length value of the current grid; traverse the entire initialized slope length array forward to calculate the slope length value of the entire grid data; traverse the entire initialized slope length array backward to calculate the slope length value of the entire grid data; if the operation of assigning the value of total to the slope length value of the current grid does not occur in the forward and reverse processes, it means that the cumulative slope length extraction is completed and the loop ends, otherwise the loop calculation is repeated from the beginning; repeat the above steps until the cumulative slope length assignment of all points is completed; Step 3: Calculate the entrance and exit slope lengths, apply for matrix space to store the exit and entrance slope lengths, and each time you traverse the elevation matrix, first determine whether the current grid is a valueless point; If there is no value, proceed to the next grid calculation. If there is a value, calculate the entrance slope length and exit slope length of the current grid according to the flow direction matrix, initial slope length matrix and cumulative slope length matrix. Repeat the above steps until all value points are found and the entrance slope length assignment is completed. Step 3.7, traverse the corresponding two-dimensional array and calculate the S factor and L factor: Step 1: Calculate the S factor and apply for a matrix space to store the slope factor (S); each time the elevation matrix is ​​traversed, first determine whether the current grid is a valueless point; if there is no value, proceed to the next grid calculation; if there is a value, calculate the S factor value based on the slope matrix and formula (4); Repeat the above steps to complete the S factor calculation for each grid; where θ is the slope; Step 2: Calculate the L factor and apply for a matrix space to store the slope length factor (L); each time the elevation matrix is ​​traversed, first determine whether the current grid is a valueless point; If there is no value, proceed to the next grid calculation. If there is a value, calculate the segment slope L value according to the entrance and exit slope length matrix according to formulas (5) and (6). When the entrance slope length is less than the exit slope length, use formula (6), otherwise use formula (6). Repeat the above steps to complete the L factor calculation of each grid. In the formula, λ is the slope length, m is the slope length index, and λout and λin are the slope lengths of the grid exit and entrance, respectively (m). L=(λ / 22.13) m (5) Step 3.8, traverse the corresponding two-dimensional array and calculate the LS factor: apply for a matrix space to store the slope length factor (LS); each time the elevation matrix is ​​traversed, first determine whether the current grid is a valueless point; If there is no value, proceed to the next grid calculation. If there is a value, calculate the LS factor value according to the L factor and S factor matrix according to formula (7); Repeat the above steps to complete the LS factor calculation for each grid; LS=L·S (7) Step 4: Convert the result data format into ASCII and extract the LS factor data without buffer: Step 4.1, convert the LS factor data of each block with buffer back to raster data through the ASCII to raster tool in Arcmap software; Step 4.2, using the clipping tool in Arcmap software, clip the LS factor raster data obtained above using each block raster data without buffer; Step 5: Data fusion to extract the LS factor of the geographic coordinate system grid of the required scale; the LS data obtained in step 4 can be fused through the mosaic tool in Arcmap software to obtain the LS factor of the grid in the geographic coordinate system of the required scale.

2. The LS factor extraction method applicable to large-scale geographic coordinate system raster data as claimed in claim 1, It is characterized in that In step 1, the data is processed in blocks for large-scale data.

3. The LS factor extraction method applicable to large-scale geographic coordinate system raster data as claimed in claim 1, It is characterized in that Before step 3, add the 1° buffer range data and convert the raster format to ASCII text; this method uses Arcmap software to complete the buffer addition and raster data conversion through mosaicking and raster to ASCII conversion.

4. The LS factor extraction method applicable to large-scale geographic coordinate system raster data as claimed in claim 1, It is characterized in that In step 3.2, the process of reading the ASCII file of the geographic coordinate system raster data and the parameter information in the LS factor extraction is as follows: Step 1: Create a structure named DemData to store the ASCII header information and the parameter information set in the LS factor extraction, and apply for a two-dimensional array to save the elevation value; Step 2: Open the raster data text file. If the opening fails, write the log and stop the execution. Step 3: First read the content in the ASCII file line by line, and record it in the format of "name-space-value" in the file header; then store each line of data read into each string, and then split the string with spaces, convert the obtained value into the type of the value and save it to the corresponding attribute of the created structure, and repeat this process until the header file is read; then record the parameter information such as flow direction coding, truncation factor and river network threshold set in LS factor extraction in the same form in the corresponding attribute of the structure.

5. The LS factor extraction method applicable to large-scale geographic coordinate system raster data as claimed in claim 1, It is characterized in that In step 3.4, the grid encoding method of the D8 flow direction algorithm takes the central grid as an example. The flow directions in eight directions around the central grid are determined according to the elevation values ​​corresponding to different grids. The directions from east, southeast, south, southwest, west, northwest, north to northeast are recorded as 1, 2, 4, 8, 16, 32, 64 and 128 respectively.

6. The LS factor extraction method applicable to large-scale geographic coordinate system raster data as claimed in claim 1, It is characterized in that In Step 2 of Step 3.6, the initial slope length value is updated by re-assigning the unit slope length according to whether the grid is truncated: for the base grid, the initial slope length is the original unit slope length value; for the truncated grid, the initial slope length is half of the original unit slope length value.

7. The LS factor extraction method applicable to large-scale geographic coordinate system raster data as claimed in claim 1, It is characterized in that Before step 5, the calculated ASCII result data is converted back to raster data and the buffer is removed. This method uses Arcmap software to complete the raster data conversion and the construction of the LS factor of the non-overlapping area through ASCII to raster conversion and clipping.

Images