A method for automatically extracting a gully network based on accumulation of surface runoff lines

Abstract

The application discloses a kind of based on surface catchment line accumulation valley network automatic extraction method, it is related to the topographic analysis technical field based on digital elevation model, including the following steps: input DEM data, judge whether there is depression in data, if there is, fill the depression as flat ground;In the DEM after completing depression filling, find the grid with maximum elevation value, not belonging to flat ground, and without calculating the number of catchment path as the current processing grid;The application takes DEM data as data source, can directly calculate the catchment line cumulative value of analysis area, then set reasonable threshold value can extract the valley network of analysis area, compared with traditional method, the degree of automation is obviously improved, and the valley network extracted in flat area has no parallel valley and other unreasonable phenomena, the result is more in line with actual situation, therefore, it has important application value in the field of agricultural planning, soil and water conservation planning and geomorphology related scientific research.

Description

Technical Field

[0001] This invention relates to the field of terrain analysis technology based on digital elevation models, specifically to an automatic extraction method for gully networks based on the accumulation of surface runoff. Background Technology

[0002] A gully is a narrow depression on the Earth's surface formed by the impact and erosion of flowing water. In areas with severe water erosion, gullies of different levels are interlocked and form a gully network. In areas with severe water erosion, such as the Loess Plateau in my country, accurately extracting the gully network is of great significance for agricultural production planning, soil and water conservation, sand control, and scientific research on landform development and evolution.

[0003] The current method for extracting gully networks is mainly based on DEM data. It utilizes surface runoff models to simulate the flow of surface runoff to generate a water system, i.e., a gully network. The basic process involves first filling depressions in the DEM data to eliminate local depressions. Then, based on the principle of maximizing elevation difference, the direction of water flow on each grid cell is determined, thus simulating the confluence path of water flow on a surface area. Next, based on this confluence path, the number of catchment cells passed through by each grid cell is calculated, and the total number is called the cumulative runoff of that grid cell. Once the cumulative runoff of a grid cell reaches a certain value range, that grid cell can be considered part of a water system or gully; that is, the location of a gully is determined by judging the cumulative runoff value of a grid cell.

[0004] The above technical process for extracting gully networks is based on the hydrological concept of confluence to determine the flow path and simulate surface runoff. It has a certain model foundation and the algorithm is simple. It can directly generate continuous flow paths and is currently the commonly used method for extracting gully networks. However, this technical method still has some problems in practical applications, such as: (1) The cumulative flow calculated by this technical method is the total number of grids through which all upstream water flows. Therefore, when the analysis area is large, especially when there is a large area of ​​flat land, the calculated cumulative flow value will be very large, and it is easy to produce the phenomenon of pseudo-gullies, that is, non-gully areas are marked as gullies; (2) This technical method calculates the grid flow direction and the cumulative flow independently. In other words, when calculating the flow direction of a grid, the characteristics of a water flow path are not considered. This is not a problem between grids with elevation differences, but in flat areas, the extracted gully morphology and spatial location will be unreasonable, such as parallel gullies.

[0005] It should be noted that many experts and scholars in fields such as geography and hydrology have proposed corresponding solutions to the above problems. For example, integrating image data to eliminate possible false valley phenomena and utilizing the spatial relationship between main valleys and tributary valleys to eliminate parallel valley phenomena can be addressed. However, these methods all involve secondary corrections based on the initial valley extraction results. The operation process is complex, inefficient, and difficult to automate, thus hindering their widespread adoption. Summary of the Invention

[0006] The purpose of this invention is to provide an automatic extraction method for gully networks based on surface runoff accumulation, so as to solve the problems mentioned in the background art.

[0007] To achieve the above objectives, the present invention provides the following technical solution:

[0008] An automatic method for extracting gully networks based on surface runoff accumulation includes the following steps:

[0009] S1: Input DEM data, determine if there are depressions in the data, and fill the depressions with flat land if they exist;

[0010] S2: In the DEM after the depression filling is completed, find the grid with the largest elevation value, which is not flat and has not calculated the number of water runoff paths as the current processing grid. If it is found, use the grid as the current processing grid, increase its cumulative value of water runoff paths by 1, and execute step S3; otherwise, execute step S6.

[0011] S3: Using the current processing raster as the center, use the principle of maximum elevation difference to find the next raster on the water runoff path. If it can be found, the newly found raster is used as the current processing raster, and its water runoff path accumulation value is increased by 1. Then continue to search in the same way. If it cannot be found, it is necessary to make a judgment and execute the corresponding steps according to the judgment result: if the processing raster belongs to the flat area, execute step S4; if the processing raster is located at the data boundary, execute step S5.

[0012] S4: Mark the first grid cell in the current catchment line that reaches the flat area as the flat area inlet grid. Determine the flat area outlet grid in the flat area and its surrounding adjacent grid cells, and generate the catchment path from the flat area inlet grid cell to the flat area outlet grid cell. At the same time, increment the cumulative value of the catchment path of all path grid cells between the flat area inlet grid cell and the outlet grid cell by 1. Then make a judgment. If the cumulative value of the catchment path of the flat area outlet grid cell is equal to 1, then the flat area outlet grid cell is used as the current processing grid cell, and continue to execute step S3. If the cumulative value of the catchment path of the flat area outlet grid cell is greater than 1, it means that the grid cell has been marked as a catchment path grid cell before this catchment path. It is necessary to increment the cumulative value of the catchment path of all grid cells between this grid cell and the outlet of this catchment path cell by 1 to complete the calculation of a catchment path, and continue to execute step S2.

[0013] S5: The water flow path reaches the data boundary, completing the calculation of one water flow path, and continue to execute step S2;

[0014] S6: After all watershed calculations are completed, a certain threshold is set according to the requirements for extracting the gully network. The grids that exceed the threshold are marked, thus obtaining the gully network result.

[0015] Preferably, the specific process of step S4 is as follows:

[0016] S41: Traverse the DEM data and mark the raster cells;

[0017] S42: Determine if there are any grid cells in the current flat area that have already calculated the cumulative value of the confluence path. If there are and the number is greater than 1, proceed to step S43. If there are and the number is 1, proceed to step S44. If not found, proceed to step S45.

[0018] S43: Determine the path from the flatland entrance grid to the grid in the current flatland region where the cumulative value of the confluence path has been calculated, and calculate the path length;

[0019] S44: The only grid cell in the current flat area that has already calculated the cumulative value of the confluence path is used as the flat area exit grid cell. Then, the confluence path from the flat area entrance grid cell to the flat area exit grid cell is generated in the manner described in step S43.

[0020] S45: Calculate the difference between the elevation values ​​of the flatland entrance grid and the adjacent grids in the current flatland area. If there is a grid with a difference greater than 0, take its maximum value. If there is only one maximum value, take the grid corresponding to the maximum value as the flatland exit grid. Then generate the confluence path from the flatland entrance grid to the flatland exit grid in the manner described in step S43. If there are multiple maximum values, generate the confluence path from the flatland entrance grid to the grid with the maximum difference in the manner described in step S43, calculate the total length of the confluence path, select the one with the smallest length as the flatland exit grid, and then take the corresponding path grid as the grid of this confluence path. If the difference between the elevation values ​​of the flatland entrance grid and the adjacent grids in the current flatland area is less than 0, proceed to step S46.

[0021] S46: When the current flat area reaches the data boundary, if there is only one raster that belongs to both the current flat area and the boundary, then this raster is used as the flat area exit raster. Then, the confluence path from the flat area entrance raster to the flat area exit raster is generated in the manner described in step S43. If there are multiple rasteres that belong to both the current flat area and the boundary, then the confluence path from the flat area entrance raster to the boundary raster is generated separately in the manner described in step S43. The total length of the confluence path is calculated, and the one with the shortest length is selected as the flat area exit raster. Then, the corresponding path raster is used as the raster of this confluence path.

[0022] Preferably, the marking method in step S41 is:

[0023] The grid with the same elevation value as the flatland entrance grid and connected to the flatland entrance grid is marked as 1, which is the current flatland area. Then, the grids adjacent to the flatland area are marked as -1, and the other grids are marked as 0.

[0024] Preferably, the specific operation method in step S43 is as follows:

[0025] S431: Set up a stack and push the flat ground entrance grid onto the stack, initialize the distance value to be calculated to 0;

[0026] S432: Determine whether the topmost grid cell in the stack is adjacent to the grid cell whose distance is to be calculated. If so, calculate the distance between the topmost grid cell and the grid cell whose distance is to be calculated, and add it to the distance value to be calculated. If the topmost grid cell in the stack is not adjacent to the grid cell whose distance is to be calculated, select all adjacent grid cells of the topmost grid cell that belong to the flat area, are not in the stack, and are not marked as non-drainage path grid cells. Calculate the distance between them and the grid cell whose distance is to be calculated. Select the adjacent grid cell with the smallest distance value and push it into the stack. Add the smallest distance to the distance value to be calculated, and then continue to execute S432.

[0027] S433: If none of the adjacent grids of the topmost grid in the current stack meet the path condition, then pop the topmost grid from the current stack and mark it as a non-catchment path grid.

[0028] S434: Continuously execute S432 and S433 until the distance between the flatland inlet grid and the marker catchment cumulative value grid is calculated. Select the grid with the smallest calculated catchment path cumulative value as the flatland outlet grid and use the corresponding path as the path from the flatland inlet grid to the flatland outlet grid.

[0029] Compared with the prior art, the beneficial effects of the present invention are:

[0030] 1. This invention uses DEM data as the data source, which can directly calculate the cumulative value of the catchment area in the analysis area. Then, by setting a reasonable threshold, the gully network of the analysis area can be extracted. Compared with traditional methods, the degree of automation is significantly improved.

[0031] 2. The method proposed in this invention, which uses the cumulative value of the water catchment line as the basis for determining whether a valley is a valley, is theoretically more reasonable than the traditional method of calculating the number of grids through which the water flows.

[0032] 3. Unlike traditional methods that independently calculate grid flow direction and runoff accumulation, the method of this invention combines the calculation of grid flow direction with the calculation of runoff paths. The biggest advantage of this approach is in the handling of flat areas. When calculating the grid flow direction in flat areas, the overall flow direction and other characteristics of the currently calculated runoff paths are fully considered. Therefore, the extracted gully network is reasonably distributed in flat areas, eliminating unreasonable phenomena such as parallel gullies. Attached Figure Description

[0033] Figure 1 This is a flowchart illustrating the method of the present invention.

[0034] Figure 2 This is a schematic diagram illustrating the calculation of the elevation difference between adjacent grid cells in the method of this invention.

[0035] Figure 3 This is a schematic diagram of a grid in a flat area in the method of the present invention.

[0036] Figure 4 This is a schematic diagram of the process of finding the entrance grid and the exit grid in the flat area in the method of the present invention.

[0037] Figure 5 This is a schematic diagram of the DEM of the experimental area in an embodiment of the present invention.

[0038] Figure 6 This is a schematic diagram of the cumulative water catchment value calculated in the experimental area in an embodiment of the present invention.

[0039] Figure 7This is a comparison diagram of the valley network extracted by the method of the present invention and the traditional method in the experimental area of ​​the present invention embodiment.

[0040] Figure 8 This is a comparison diagram of local areas of the valley network extracted by the method of the present invention and the traditional method in the experimental area of ​​the present invention embodiment. Detailed Implementation

[0041] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0042] Please see Figure 1-8 The present invention provides the following technical solution:

[0043] An automatic method for extracting gully networks based on surface runoff accumulation includes the following steps:

[0044] S1: Input DEM data, determine if there are depressions in the data, and fill the depressions with flat land if they exist;

[0045] S2: In the DEM after the depression filling is completed, find the grid with the largest elevation value, which is not flat and has not calculated the number of water runoff paths as the current processing grid. If it is found, use the grid as the current processing grid, increase its cumulative value of water runoff paths by 1, and execute step S3; otherwise, execute step S6.

[0046] S3: Using the current processing raster as the center, use the principle of maximum elevation difference to find the next raster on the water runoff path. If it can be found, the newly found raster is used as the current processing raster, and its water runoff path accumulation value is increased by 1. Then continue to search in the same way. If it cannot be found, it is necessary to make a judgment and execute the corresponding steps according to the judgment result: if the processing raster belongs to the flat area, execute step S4; if the processing raster is located at the data boundary, execute step S5.

[0047] like Figure 2 As shown, the method for calculating the elevation difference is to use the right, bottom, left, and top four grids of the currently processed grid according to the formula...

[0048] Formula 1 is used for calculation, and the remaining four grid cells are calculated according to Formula 2:

[0049]

[0050]

[0051] In the formula, Δh is the elevation difference between two grids, h0 is the currently processed grid, and h i For adjacent grid cells, r is the grid size.

[0052] S4: Mark the first grid cell in the current catchment line that reaches the flat area as the flat area inlet grid. Determine the flat area outlet grid in the flat area and its surrounding adjacent grid cells, and generate the catchment path from the flat area inlet grid cell to the flat area outlet grid cell. At the same time, increment the cumulative value of the catchment path of all path grid cells between the flat area inlet grid cell and the outlet grid cell by 1. Then make a judgment. If the cumulative value of the catchment path of the flat area outlet grid cell is equal to 1, then the flat area outlet grid cell is used as the current processing grid cell, and continue to execute step S3. If the cumulative value of the catchment path of the flat area outlet grid cell is greater than 1, it means that the grid cell has been marked as a catchment path grid cell before this catchment path. It is necessary to increment the cumulative value of the catchment path of all grid cells between this grid cell and the outlet of this catchment path cell by 1 to complete the calculation of a catchment path, and continue to execute step S2.

[0053] The specific process of step S4 is as follows:

[0054] S41: Traverse the DEM data and mark the grids; the marking method in step S41 is as follows: mark the grids with the same elevation value as the flat land entrance grid and connected to the flat land entrance grid as 1, which is the current flat land area, and then mark the grids adjacent to the flat land area as -1, and mark the other grids as 0.

[0055] S42: Determine if there are any grid cells in the current flat area that have already calculated the cumulative value of the confluence path. If there are and the number is greater than 1, proceed to step S43. If there are and the number is 1, proceed to step S44. If not found, proceed to step S45.

[0056] S43: Determine the path from the flat land inlet grid to the grid in the current flat land area where the cumulative value of the confluence path has been calculated, and calculate the path length; the specific operation method in step S43 is as follows:

[0057] S431: Set up a stack and push the flat ground entrance grid onto the stack, initialize the distance value to be calculated to 0;

[0058] S432: Determine whether the topmost grid cell in the stack is adjacent to the grid cell whose distance is to be calculated. If so, calculate the distance between the topmost grid cell and the grid cell whose distance is to be calculated, and add it to the distance value to be calculated. If the topmost grid cell in the stack is not adjacent to the grid cell whose distance is to be calculated, select all adjacent grid cells of the topmost grid cell that belong to the flat area, are not in the stack, and are not marked as non-drainage path grid cells. Calculate the distance between them and the grid cell whose distance is to be calculated. Select the adjacent grid cell with the smallest distance value and push it into the stack. Add the smallest distance to the distance value to be calculated, and then continue to execute S432.

[0059] S433: If none of the adjacent grids of the topmost grid in the current stack meet the path condition, then pop the topmost grid from the current stack and mark it as a non-catchment path grid.

[0060] S434: Continuously execute S432 and S433 until the distance between the flatland inlet grid and the marker catchment cumulative value grid is calculated. Select the grid with the smallest calculated catchment path cumulative value as the flatland outlet grid and use the corresponding path as the path from the flatland inlet grid to the flatland outlet grid.

[0061] S44: The only grid cell in the current flat area that has already calculated the cumulative value of the confluence path is used as the flat area exit grid cell. Then, the confluence path from the flat area entrance grid cell to the flat area exit grid cell is generated in the manner described in step S43.

[0062] S45: Calculate the difference between the elevation values ​​of the flatland entrance grid and the adjacent grids in the current flatland area. If there is a grid with a difference greater than 0, take its maximum value. If there is only one maximum value, take the grid corresponding to the maximum value as the flatland exit grid. Then generate the confluence path from the flatland entrance grid to the flatland exit grid in the manner described in step S43. If there are multiple maximum values, generate the confluence path from the flatland entrance grid to the cumulative value of the waterline of the grid with the maximum difference in the manner described in step S43, and calculate the total length of the confluence path. Select the one with the smallest length as the flatland exit grid, and then take the corresponding path grid as the grid of this confluence path. If the difference between the elevation values ​​of the flatland entrance grid and the adjacent grids in the current flatland area is less than 0, proceed to step S46.

[0063] like Figure 4 As shown, the method for calculating the total length of the confluence path is as follows: (1) Set up a stack and push the flat ground entrance grid onto the stack, and initialize the distance value to be calculated to 0; (2) Determine whether the topmost grid in the stack is adjacent to the grid whose distance is to be calculated. If so, calculate the distance between the topmost grid in the stack and the grid whose distance is to be calculated, and add it to the distance value to be calculated. The formula for calculating the distance between two grids is:

[0064]

[0065] In the formula, d represents the distance between two grid cells, r is the grid cell size, and (i1,j1) and (i2,j2) are the row and column numbers of the two grid cells, respectively.

[0066] S46: When the current flat area reaches the data boundary, if there is only one raster that belongs to both the current flat area and the boundary, then this raster is used as the flat area exit raster. Then, the confluence path from the flat area entrance raster to the flat area exit raster is generated in the manner described in step S43. If there are multiple rasteres that belong to both the current flat area and the boundary, then the confluence path from the flat area entrance raster to the boundary raster is generated separately in the manner described in step S43. The total length of the confluence path is calculated, and the one with the shortest length is selected as the flat area exit raster. Then, the corresponding path raster is used as the raster of this confluence path.

[0067] S5: The water flow path reaches the data boundary, completing the calculation of one water flow path, and continue to execute step S2;

[0068] S6: After all watershed calculations are completed, a certain threshold is set according to the requirements for extracting the gully network. The grids that exceed the threshold are marked, thus obtaining the gully network result.

[0069] To verify the applicability of the method of the present invention in extracting gully networks, a small watershed in a hilly area in southeastern my country was selected as the experimental area. This experimental area has a complete gully system, with gentle terrain undulations and relatively wide main gullies, such as... Figure 5 As shown, when using traditional methods to extract valley networks, unreasonable phenomena such as parallel valleys are easily observed. Therefore, this method can be used to verify the effectiveness of the method of the present invention.

[0070] Arrangement 1: For the example area, following the calculation process of the present invention, input DEM data and calculate the cumulative value of the catchment area (see results). Figure 6 );

[0071] Step 2: Referring to the DEM data of the example area, determine a threshold of 300, and mark the watercourses with a cumulative value greater than 300. This gives the result of the valley network extraction in the example area (see results). Figure 7 );

[0072] Step 3: Using traditional methods, including depression filling, flow direction calculation, and runoff accumulation calculation (implemented using ArcGIS software), and employing an experimental method with a runoff accumulation threshold of 3000, the extracted gully network is basically similar to the results extracted by the method of this invention in Step 2 (see results). Figure 7 );

[0073] Step 4: Compare the valley network extracted by the method of this invention with the valley network extracted by the traditional method, wherein... Figure 7 It is a comparison of the complete example area. Figure 8 It is a local area comparison near the outlet of a small watershed, in Figure 7 , Figure 8 In the diagram, the solid black lines represent the valley network extracted by the method of this invention, and the dashed lines represent the valley network extracted by the traditional method.

[0074] Through the above comparative analysis, it can be seen that the process of extracting valley networks based on the method of the present invention is simpler, and the valley networks extracted by the method of the present invention are more reasonable in terms of spatial morphology. Therefore, the method of the present invention has better application value in extracting valley networks based on DEM data.

[0075] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for automatically extracting gully networks based on surface runoff accumulation, characterized in that, Includes the following steps: S1: Input DEM data, determine if there are depressions in the data, and fill the depressions with flat land if they exist; S2: In the DEM after the depression filling is completed, find the grid with the largest elevation value, which is not flat and has not calculated the number of water runoff paths as the current processing grid. If it is found, use the grid as the current processing grid, increase its cumulative value of water runoff paths by 1, and execute step S3; otherwise, execute step S6. S3: Using the current processing raster as the center, use the principle of maximum elevation difference to find the next raster on the water runoff path. If it can be found, the newly found raster is used as the current processing raster, and its water runoff path accumulation value is increased by 1. Then continue to search in the same way. If it cannot be found, it is necessary to make a judgment and execute the corresponding steps according to the judgment result: if the processing raster belongs to the flat area, execute step S4; if the processing raster is located at the data boundary, execute step S5. S4: Mark the first grid cell in the current catchment line that reaches the flat area as the flat area inlet grid. Determine the flat area outlet grid in the flat area and its surrounding adjacent grid cells, and generate the catchment path from the flat area inlet grid cell to the flat area outlet grid cell. At the same time, increment the cumulative value of the catchment path of all path grid cells between the flat area inlet grid cell and the outlet grid cell by 1. Then make a judgment. If the cumulative value of the catchment path of the flat area outlet grid cell is equal to 1, then the flat area outlet grid cell is used as the current processing grid cell, and continue to execute step S3. If the cumulative value of the catchment path of the flat area outlet grid cell is greater than 1, it means that the grid cell has been marked as a catchment path grid cell before this catchment path. It is necessary to increment the cumulative value of the catchment path of all grid cells between this grid cell and the outlet of this catchment path cell by 1 to complete the calculation of a catchment path, and continue to execute step S2. The specific process of step S4 is as follows: S41: Traverse the DEM data and mark the raster cells; S42: Determine if there are any grid cells in the current flat area that have already calculated the cumulative value of the confluence path. If there are and the number is greater than 1, proceed to step S43. If there are and the number is 1, proceed to step S44. If not found, proceed to step S45. S43: Determine the path from the flatland entrance grid to the grid in the current flatland region where the cumulative value of the confluence path has been calculated, and calculate the path length; S44: The only grid cell in the current flat area that has already calculated the cumulative value of the confluence path is used as the flat area exit grid cell. Then, the confluence path from the flat area entrance grid cell to the flat area exit grid cell is generated in the manner described in step S43. S45: Calculate the difference between the elevation values ​​of the flatland entrance grid and the adjacent grids in the current flatland area. If there is a grid with a difference greater than 0, take its maximum value. If there is only one maximum value, take the grid corresponding to the maximum value as the flatland exit grid. Then generate the confluence path from the flatland entrance grid to the flatland exit grid in the manner described in step S43. If there are multiple maximum values, generate the confluence path from the flatland entrance grid to the grid with the maximum difference in the manner described in step S43, calculate the total length of the confluence path, select the one with the smallest length as the flatland exit grid, and then take the corresponding path grid as the grid of this confluence path. If the difference between the elevation values ​​of the flatland entrance grid and the adjacent grids in the current flatland area is less than 0, proceed to step S46. S46: When the current flat area reaches the data boundary, if there is only one raster that belongs to both the current flat area and the boundary, then this raster is used as the flat area exit raster. Then, the confluence path from the flat area entrance raster to the flat area exit raster is generated in the manner described in step S43. If there are multiple rasteres that belong to both the current flat area and the boundary, then the confluence path from the flat area entrance raster to the boundary raster is generated separately in the manner described in step S43. The total length of the confluence path is calculated, and the one with the shortest length is selected as the flat area exit raster. Then, the corresponding path raster is used as the raster of this confluence path. S5: The water flow path reaches the data boundary, completing the calculation of one water flow path, and continue to execute step S2; S6: After all watershed calculations are completed, a certain threshold is set according to the requirements for extracting the gully network. The grids that exceed the threshold are marked, thus obtaining the gully network result.

2. The method for automatic extraction of gully networks based on surface runoff accumulation according to claim 1, characterized in that, The marking method in step S41 is as follows: The grid with the same elevation value as the flat land entrance grid and connected to the flat land entrance grid is marked as 1, which is the current flat land area. Then, the grid adjacent to the flat land area is marked as -1, and the other grids are marked as 0.

3. The method for automatic extraction of gully networks based on surface runoff accumulation according to claim 1, characterized in that, The specific operation method in step S43 is as follows: S431: Set up a stack and push the flat ground entrance grid onto the stack, initialize the distance value to be calculated to 0; S432: Determine whether the topmost grid cell in the stack is adjacent to the grid cell whose distance is to be calculated. If so, calculate the distance between the topmost grid cell and the grid cell whose distance is to be calculated, and add it to the distance value to be calculated. If the topmost grid cell in the stack is not adjacent to the grid cell whose distance is to be calculated, select all adjacent grid cells of the topmost grid cell that belong to the flat area, are not in the stack, and are not marked as non-drainage path grid cells. Calculate the distance between them and the grid cell whose distance is to be calculated. Select the adjacent grid cell with the smallest distance value and push it into the stack. Add the smallest distance to the distance value to be calculated, and then continue to execute S432. S433: If none of the adjacent grids of the topmost grid in the current stack meet the path condition, then pop the topmost grid from the current stack and mark it as a non-catchment path grid. S434: Continuously execute S432 and S433 until the distance between the flatland inlet grid and the marker catchment cumulative value grid is calculated. Select the grid with the smallest calculated catchment path cumulative value as the flatland outlet grid and use the corresponding path as the path from the flatland inlet grid to the flatland outlet grid.

Images