A New Algorithm for Calculating the Flow Path Curvature ( C ) from the Square-Grid Digital Elevation Model (DEM)

: This paper proposes a ﬂow-path-network-based (FPN-based) algorithm, constructed from a square-grid digital elevation model (DEM) to improve the simulation of the ﬂow path curvature ( C ). First, the ﬂow-path network model was utilized to obtain an FPN. Then, a ﬂow-path-network-ﬂow-path-curvature (FPN-C) algorithm was proposed to estimate C from the FPN. The experiments consisted of two sections: (1) quantitatively evaluating the accuracy using 5 m DEMs generated from the mathematical ellipsoid and Gauss models, and (2) qualitatively assessing the accuracy using a 30 m DEM of a real-world complex region. The three algorithms proposed by Evans (1980), Zevenbergen and Throne (1987), and Shary (1995) were used to validate the accuracy of the new algorithm. The results demonstrate that the C value of the proposed algorithm was generally closer to the theoretical C value derived from two mathematical surfaces. The root mean standard error (RMSE) and mean absolute error (MAE) of the new method are 0.0014 and 0.0002 m, reduced by 42% and 82% of that of the third algorithm on the ellipsoid surface, respectively. The RMSE and MAE of the presented method are 0.0043 and 0.0025 m at best, reduced by up to 35% and 14% of that of the former two algorithms on the Gauss surface, respectively. The proposed algorithm generally produces better spatial distributions of C on di ﬀ erent terrain surfaces.

For a long time, a variety of curvatures developed by researchers have been used to meet the requirements of practical applications [25][26][27][28][29][30], such as the mean curvature, maximum curvature, minimum curvature, plan curvature, profile curvature, tangential curvature, flow path curvature, and so on. Less acknowledged is the flow path curvature (C), which is also known as the rotor curvature or streamline curvature [31]. It measures the rate of the change of the flow paths along the horizontal direction, and describes the degree of twisting of the flow lines [4,[32][33][34]. Although the flow path curvature was not considered in the complete classification system of curvatures constructed by Shary [35] at first, it has been widely utilized in the theory of the electromagnetic field [36]. For example, f yy )/ f x 2 + f y 2 3/2 (units: m −1 ) by Shary [26]. It is measuring the twist of flow lines. When it is larger than zero, flow lines rotate clockwise. When it is smaller than zero, flow lines rotate counterclockwise, otherwise the flow line does not swing along the straight line [3,42,43]. The flow path curvature is often derived from the square-grid digital elevation model (DEM), which is valued for its simple structure and continuity in the representation of topographic surfaces.
The commonly-used algorithms use the center grid cell and its eight surrounding grid cells based on a moving window of 3 × 3 to calculate the C of the center grid cell. The elevation value of the nine grid cells is approximated by the differentiate operation or local fitting curve. For example, the method proposed by Evans [44] uses the six-parameter second-order polynomial to represent the terrain surface, and derives different partial derivatives by the least-squares fitting algorithm to calculate C. It has a high accuracy when considering the smoothing of the high-frequency noise of the DEM [24,45]. Under the above method, Zevenbergen and Thorne [46] proposed a method which utilizes the partial nine-parameter fourth-order polynomial to describe the surface, and so the fitting curve can be passed through the nine grid cells and obtain the solely different partial derivatives for the C calculation. Its aim was to enhance the accuracy of the different partial derivatives, but it has not achieved the desired results because it lacks the smoothing and denoising effect of DEM [43]. The high-order polynomial interpolation method may result in incorrect topographic features [24].
Because the general second-order surface does not pass through all of the nine grid cells, the method presented by Shary [35] employs the constrained five-parameter second-order polynomial to fit the surface, and is also based on the least-squares fitting algorithm to derive the different partial derivatives for estimating C. It is similar to the Evans algorithm, except for the different averaging processes. Considering the equidistant distribution characteristics of a regular grid, Moore, et al. [47] proposed a difference method using the numerical differentiation to calculate partial derivatives for the C estimation. It directly uses the elevation of the center grid cell and eight neighboring grid cells to derive the different partial derivatives for estimating C. It is similar to the algorithm proposed by Zevenbergen ISPRS Int. J. Geo-Inf. 2020, 9, 510 3 of 25 and Thorne [46], but they calculated the second-order partial derivatives using different methods. In order to improve accuracy and adaptability to the different terrain surfaces, Shary, Sharaya and Mitusov [30] proposed the modified Evans-Young method, which is based on the Evans algorithm, to calculate the curvature after using filters to handle the center grid cell. The above mentioned methods utilize the local terrain surface representation to derive the partial derivatives by the different interpolation algorithm and a moving window. They are more helpful to consistently extract the local curvature, but are not suitable for complex topography regions for higher-scale problems in terrain analysis [48]. Moreover, the accuracy of these algorithms is affected by interpolation errors, and there are difficulties in selecting suitable algorithms for the applications.
To overcome the aforementioned disadvantages, we present a flow-path-network-based (FPN-based) algorithm to derive the flow path curvature (C) based on the vector-based approach. It uses the flow-path network model [49] to generate the one-dimensional flow path network (FPN). Then, a new flow-path-network-flow-path-curvature (FPN-C) algorithm is presented to calculate the C from the FPN. It aims to improve the calculation accuracy of the C and avoid the interpolation error and the choice of the calculation algorithm in practical applications. The experiments consisted of two sections: (1) quantitatively evaluating the accuracy of the new algorithm on the 5 m DEMs generated from the mathematical ellipsoid and Gauss surface models, and (2) qualitatively assessing the accuracy of the proposed method by using a real-world DEM of a hilly plateau and mountainous region.
The structure of the paper is arranged as follows. The methods of the FPN-based approach are presented in Section 2. Section 3 describes the experiments. The experiment results are shown in Section 4. The accuracy of the proposed approach is discussed in Section 5. Section 6 concludes the paper and illustrates directions for further research.

Methods
The methods consist of two sections in this paper: (1) obtaining a flow path network (FPN) using the flow-path network model; (2) proposing an FPN-C algorithm to calculate the flow path curvature (C) from the FPN. The detailed process of the FPN-based algorithm is shown in Figure 1. methods. In order to improve accuracy and adaptability to the different terrain surfaces, Shary, Sharaya and Mitusov [30] proposed the modified Evans-Young method, which is based on the Evans algorithm, to calculate the curvature after using filters to handle the center grid cell. The above mentioned methods utilize the local terrain surface representation to derive the partial derivatives by the different interpolation algorithm and a moving window. They are more helpful to consistently extract the local curvature, but are not suitable for complex topography regions for higher-scale problems in terrain analysis [48]. Moreover, the accuracy of these algorithms is affected by interpolation errors, and there are difficulties in selecting suitable algorithms for the applications.
To overcome the aforementioned disadvantages, we present a flow-path-network-based (FPNbased) algorithm to derive the flow path curvature (C) based on the vector-based approach. It uses the flow-path network model [49] to generate the one-dimensional flow path network (FPN). Then, a new flow-path-network-flow-path-curvature (FPN-C) algorithm is presented to calculate the C from the FPN. It aims to improve the calculation accuracy of the C and avoid the interpolation error and the choice of the calculation algorithm in practical applications. The experiments consisted of two sections: (1) quantitatively evaluating the accuracy of the new algorithm on the 5 m DEMs generated from the mathematical ellipsoid and Gauss surface models, and (2) qualitatively assessing the accuracy of the proposed method by using a real-world DEM of a hilly plateau and mountainous region.
The structure of the paper is arranged as follows. The methods of the FPN-based approach are presented in Section 2. Section 3 describes the experiments. The experiment results are shown in Section 4. The accuracy of the proposed approach is discussed in Section 5. Section 6 concludes the paper and illustrates directions for further research.

Methods
The methods consist of two sections in this paper: (1) obtaining a flow path network (FPN) using the flow-path network model; (2) proposing an FPN-C algorithm to calculate the flow path curvature (C) from the FPN. The detailed process of the FPN-based algorithm is shown in Figure 1.

Obtaining a Flow Path Network (FPN) Using the Flow-Path Network Model
In this paper, the flow path network (FPN) is tracked by the flow-path network model [49], and its detailed steps are shown in Figure 1. First, a no-depression DEM was acquired by filling the sinks and local pits of the original DEM. Second, the triangular facet network algorithm [50] was used to construct the triangular facet network (TFN). Third, the flow direction of the triangular facets over the TFN was determined by its aspect and slope as shown in Figure 2, which were calculated by the equations presented by Zhou, Pilesjö and Chen [50]. When the coordinate values of the three vertices of a triangular facet were assumed as p 1 (x 1 , y 1 , z 1 ), p 2 (x 2 , y 2 , z 2 ), and p 3 (x 3 , y 3 , z 3 ), the plane equation of the facet could be specified as z = a * x + b * y + c. Here, a, b, and c could be derived from Equation (1). The aspect (α) and slope (β) could be derived from Equation (2). Thus, the flow direction over the triangular facet could be represented by a vector whose direction and length were determined by the aspect and slope, respectively. The process of estimating the flow direction was different from that of Terrain Analysis Using Digital Elevation Models (TauDEM). This is because the latter utilizes the multiple flow direction (D∞) algorithm to estimate the flow direction which is represented as the direction of the steepest downward slope over the eight triangular facets centered on a grid cell [51]. The downward slope of each triangular facet is symbolized by a vector whose direction and length are determined by the ratio of elevation change to length in the x direction and y direction, respectively [52]. It obtained the flow direction using grid-based DEM, and the method in this paper estimated the flow direction based on the three-dimensional vector facet. Fourth, the random resampling algorithm [53] was used to obtain the flow source points from the original DEM. Finally, combining the flow direction of the triangular facets over the TFN with the flow source points, an FPN was tracked based on the flow-path network model; its detailed steps have already been described in this paper [49].
ISPRS Int. J. Geo-Inf. 2020, 9, x FOR PEER REVIEW 4 of 25 In this paper, the flow path network (FPN) is tracked by the flow-path network model [49], and its detailed steps are shown in Figure 1. First, a no-depression DEM was acquired by filling the sinks and local pits of the original DEM. Second, the triangular facet network algorithm [50] was used to construct the triangular facet network (TFN). Third, the flow direction of the triangular facets over the TFN was determined by its aspect and slope as shown in Figure 2, which were calculated by the equations presented by Zhou, Pilesjö and Chen [50]. When the coordinate values of the three vertices of a triangular facet were assumed as 1 ( 1 , 1 , 1 ) , 2 ( 2 , 2 , 2 ), and 3 ( 3 , 3 , 3 ) , the plane equation of the facet could be specified as = * + * + . Here, , , and could be derived from Equation (1). The aspect ( ) and slope ( ) could be derived from Equation (2). Thus, the flow direction over the triangular facet could be represented by a vector whose direction and length were determined by the aspect and slope, respectively. The process of estimating the flow direction was different from that of Terrain Analysis Using Digital Elevation Models (TauDEM). This is because the latter utilizes the multiple flow direction (D∞) algorithm to estimate the flow direction which is represented as the direction of the steepest downward slope over the eight triangular facets centered on a grid cell [51]. The downward slope of each triangular facet is symbolized by a vector whose direction and length are determined by the ratio of elevation change to length in the direction and direction, respectively [52]. It obtained the flow direction using grid-based DEM, and the method in this paper estimated the flow direction based on the three-dimensional vector facet. Fourth, the random resampling algorithm [53] was used to obtain the flow source points from the original DEM. Finally, combining the flow direction of the triangular facets over the TFN with the flow source points, an FPN was tracked based on the flow-path network model; its detailed steps have already been described in this paper [49].
(1) The flow path curvature (C) of a grid cell was assumed to be derived from the flow line passing the grid cell. There may be numerous flow lines through a grid cell, and they are parallel for the same flow direction of each grid cell. Therefore, it was necessary to select the suitable flow line from the FPN.
The selection method is as follows: (1) Judge whether the number of flow lines passing through the calculated grid cell is equal to 1. If it is equal to 1, there will be only one flow line passing through the grid cell, and the flow line will be regarded as the suitable flow line, otherwise; (2) find all of the flow lines passing through the calculated grid cell, namely, the pass-lines; (3) select the flow lines throughout the calculated grid cell from the pass-lines, namely, the through-lines. The standard throughout a grid cell depends on the length of the flow line within the grid cell and whether it is larger than the length of the grid cell; (4) regard the longest line among the through-lines as the suitable flow line.

Smoothing the Flow Line by the B-Spline Interpolation Method
The most suitable flow line within a grid cell consists of several break points; each broken line may have significantly different curvature, and it is not reasonable to expect that any of the broken lines can be utilized to estimate C. Thus, we smoothed the flow line so that the whole line in a grid cell could be used to derive an accurate flow path curvature (C).
The spline interpolation method is commonly used to obtain smooth curves in the mathematical [54][55][56], physical [57][58][59][60], and other fields [61][62][63][64]. In this study, the flow line was smoothed by the B-spine interpolation method (B-spline method), which is suitable for handling the multivalued functions that may appear on the flow line. To keep the overall smoothness and a continuous slope and curvature at the break points, the cubic B-spine interpolation algorithm was utilized to smooth the flow line by using these break points within the calculated grid cell. Moreover, the flow line was cut up by a threshold to reduce the likelihood of overfitting. To prevent the lines in the grid cell being too short and the number of the break points not being enough, the threshold should not be too small.
According to the principle of the cubic B-spline interpolation algorithm, we could obtain a B-spline curve (P(u)), which is a piecewise function, as shown in Equation (3). If there were n + 1 points and a node vector (U = u 0 , u 1 , . . . , u n+k+1 , (k = 3)) used to smooth the curve, there would be an n basic function. Each function (N i,k (x), (i = 0, 1, . . . , n)) could be defined as Equation (4), and the operational relationship between the basic functions is shown in Figure 3. Under the rules of interpolation continuity and differential continuity, we could obtain Equation (5). Combing the equation with Equations (3) and (4), we could calculate the P 0 , P 1 , P 2 , . . . , P n to obtain the P(u). Figure 4 illustrates a smoothing flow line by the cubic B-spline method.

Fitting the Circle by the Least Square Algorithm
The premise of fitting a circle is selecting a series of points from the smoothing flow line in the calculated grid cell. Therefore, it was key to choose the proper points from the line in the grid cell.
In this paper, the points were obtained by equally dividing the smoothing flow line, such as the black points shown in Figure 5. The least square algorithm with the greatest effect was utilized to fit the circle (the assumption of the circle equation is 2 + 2 + * + * + = 0) by the several points in the calculated grid cell (the assumptions of these points are ( , ), = 1,2,3, … , ). Once the value of , , and was determined, the circle was obtained. According to the principle of the least square method, we could obtain the objective function (as shown in Equation (6)). The optimal circle was matched when the objective function reached its minimum, and , , and could be acquired from Equation (7).

Fitting the Circle by the Least Square Algorithm
The premise of fitting a circle is selecting a series of points from the smoothing flow line in the calculated grid cell. Therefore, it was key to choose the proper points from the line in the grid cell.
In this paper, the points were obtained by equally dividing the smoothing flow line, such as the black points shown in Figure 5. The least square algorithm with the greatest effect was utilized to fit the circle (the assumption of the circle equation is 2 + 2 + * + * + = 0) by the several points in the calculated grid cell (the assumptions of these points are ( , ), = 1,2,3, … , ). Once the value of , , and was determined, the circle was obtained. According to the principle of the least square method, we could obtain the objective function (as shown in Equation (6)). The optimal circle was matched when the objective function reached its minimum, and , , and could be acquired from Equation (7).

Fitting the Circle by the Least Square Algorithm
The premise of fitting a circle is selecting a series of points from the smoothing flow line in the calculated grid cell. Therefore, it was key to choose the proper points from the line in the grid cell.
In this paper, the points were obtained by equally dividing the smoothing flow line, such as the black points shown in Figure 5. The least square algorithm with the greatest effect was utilized to fit the circle (the assumption of the circle equation is x 2 + y 2 + a * x + b * y + c = 0) by the several points in the calculated grid cell (the assumptions of these points are (x i , y i ), i = 1, 2, 3, . . . , n). Once the value of a,b, and c was determined, the circle was obtained. According to the principle of the least square method, we could obtain the objective function (as shown in Equation (6)). The optimal circle was matched when the objective function reached its minimum, and a, b, and c could be acquired from Equation (7).
where f (a, b, c) denotes the objective function, i denotes the i th point, and x i and y i are the coordinate values of the i th point in the x direction and y direction, respectively. ( , , ) = ∑ ( 2 + 2 + * + * + ) 2 where ( , , ) denotes the objective function, denotes the ℎ point, and and are the coordinate values of the ℎ point in the x direction and y direction, respectively.
The method can use more than three points to complete the circle fitting, and when more points are applied it is possible to achieve a better performance. Thus, the standard of dividing the line evenly in this paper was ensuring that there were more than ten points within the calculated grid cell. The red circle in Figure 5 was acquired by these points in the black grid cell, based on the above algorithm.

Calculating the C by the Fitting Circle
According to Figure 5, the flow path curvature (C) of each grid cell was computed by the radius of the fitting circle (the assumption of the radius is ) and any three ordered points (the assumptions of the three blue points are 0 ( 0 , 0 ), 1 ( 1 , 1 ), 2 ( 2 , 2 )) in the grid cell. The circle radius could be derived from three parameters of the circle (assumption of the parameters are a, b, c) by the equation ( = 0.5 * √ 2 + 2 − 4 * ). Thus, the absolute value of C was the reciprocal of the radius, namely, 1/ . Its plus-minus sign was determined by the cross product of 0 1 ⃗⃗⃗⃗⃗⃗⃗⃗ and 1 2 ⃗⃗⃗⃗⃗⃗⃗⃗ . If 0 1 ⃗⃗⃗⃗⃗⃗⃗⃗ • 1 2 ⃗⃗⃗⃗⃗⃗⃗⃗ was greater than zero, C would be −1/ , which denoted that the counterclockwise flow rotated in the flow direction, otherwise, C would be 1/ , which denoted that the clockwise flow rotated in the flow direction. In addition, if C was zero, the flow would not have any swing.

Experiments
We conducted experiments using two mathematical models to quantitatively evaluate the accuracy of the new algorithm in this study. Moreover, the new algorithm was applied to a realworld DEM of a hilly plateau and mountainous area, located in the central Ganzi Tibetan Autonomous Prefecture of Sichuan Province to qualitatively assess its accuracy.

Quantitative Experiment
The ellipsoid surface model (Equation (8)) and Gauss surface model (Equation (9)) [65,66] were selected to generate the DEMs with a 5 m resolution. We utilized the DEMs generated from four The method can use more than three points to complete the circle fitting, and when more points are applied it is possible to achieve a better performance. Thus, the standard of dividing the line evenly in this paper was ensuring that there were more than ten points within the calculated grid cell. The red circle in Figure 5 was acquired by these points in the black grid cell, based on the above algorithm.

Calculating the C by the Fitting Circle
According to Figure 5, the flow path curvature (C) of each grid cell was computed by the radius of the fitting circle (the assumption of the radius is r) and any three ordered points (the assumptions of the three blue points are P 0 (x 0 , y 0 ), P 1 (x 1 , y 1 ), P 2 (x 2 , y 2 )) in the grid cell. The circle radius could be derived from three parameters of the circle (assumption of the parameters are a, b, c) by the equation . Thus, the absolute value of C was the reciprocal of the radius, namely, 1/r. Its plus-minus sign was determined by the cross product of → P 0 P 1 and → P 1 P 2 . If → P 0 P 1 · → P 1 P 2 was greater than zero, C would be −1/r, which denoted that the counterclockwise flow rotated in the flow direction, otherwise, C would be 1/r, which denoted that the clockwise flow rotated in the flow direction. In addition, if C was zero, the flow would not have any swing.

Experiments
We conducted experiments using two mathematical models to quantitatively evaluate the accuracy of the new algorithm in this study. Moreover, the new algorithm was applied to a real-world DEM of a hilly plateau and mountainous area, located in the central Ganzi Tibetan Autonomous Prefecture of Sichuan Province to qualitatively assess its accuracy.

Quantitative Experiment
The ellipsoid surface model (Equation (8)) and Gauss surface model (Equation (9)) [65,66] were selected to generate the DEMs with a 5 m resolution. We utilized the DEMs generated from four ellipsoid surface models (namely, E1, E2, E3, and E4) and four Gauss surface models (namely, G1, G2, G3, and G4) with different complexities to validate the accuracy of the new method. Table 1 shows the parameters of the eight mathematical surfaces. The theoretical flow path curvature (C) could be derived from the mathematical formulas of the equations summarized in Table 2.

Surface Types Formulas Describing the Surfaces and Calculating Theoretical the Theoretical C Value
Ellipsoid 1 Theoretical flow path curvature (C) value at the point (x p , y p , y p ).
In this paper, the 5 m DEMs were resampled to obtain 15,031 and 40,000 flow source points at random on the ellipsoid and Gauss surfaces, respectively. The threshold values of 30, 50, 100, and 150 m were used to cut up the flow line over the FPN on the E1 and G1 in the process of choosing the optimal threshold. During the circle fitting by the least square algorithm, we ensured that there were more than ten points within the calculated grid cell. Next, we selected the optimum threshold value to estimate the C value on E2, E3, E4, G2, G3, and G4. Finally, the simulated C value was compared with the theoretical C value by the root mean standard error (RMSE) and mean absolute error (MAE), to validate its accuracy. The RMSE and MAE are expressed as follows: ISPRS Int. J. Geo-Inf. 2020, 9, 510 9 of 25 where C i denotes the theoretical value, C i denotes the simulated value, n denotes the number of grid cells, and i denotes the i th grid cell.

Comparison Algorithms
In this study, three common published methods were compared to the proposed method. These algorithms used the elevation values of the surrounding cells to compute the C of the central cell using a 3 × 3 window, as shown in Figure 6. more than ten points within the calculated grid cell. Next, we selected the optimum threshold value to estimate the C value on E2, E3, E4, G2, G3, and G4. Finally, the simulated C value was compared with the theoretical C value by the root mean standard error (RMSE) and mean absolute error (MAE), to validate its accuracy. The RMSE and MAE are expressed as follows: (11) where denotes the theoretical value, ′ denotes the simulated value, denotes the number of grid cells, and denotes the ℎ grid cell.

Comparison Algorithms
In this study, three common published methods were compared to the proposed method. These algorithms used the elevation values of the surrounding cells to compute the C of the central cell using a 3 × 3 window, as shown in Figure 6. The first is the method proposed by Evans [44]. The expression of the fitting surface and the equation for calculating C are expressed as Equation (12). The first is the method proposed by Evans [44]. The expression of the fitting surface and the equation for calculating C are expressed as Equation (12).
where Z i denotes the elevation of the i th grid cell.
The second is the method proposed by Zevenbergen and Thorne [46]. The expression of the fitting surface and the equation for calculating C are expressed as Equation (13).
where Z i denotes the elevation of the i th grid cell. The third is the method proposed by Shary [35]. The functional expression of the fitting surface and the equation for calculating C are expressed as Equation (14).
where Z i denotes the elevation of the i th grid cell.
The three above-mentioned algorithms are referred to hereafter as the Evans, Zevenbergen, and Shary algorithms, respectively. C++ was used to implement the Evans, Zevenbergen, and Shary algorithms, and the proposed algorithm was implemented using C++ and Python.

Results
The experimental results consist of two sections: (1) the results of the quantitative assessment of the accuracy of the flow path curvature (C) derived from the mathematical models using RMSE and MAE; (2) the results of the qualitative estimation of the accuracy of C derived from the real-world DEM using the error of the spatial distribution.

Results
The experimental results consist of two sections: (1) the results of the quantitative assessment of the accuracy of the flow path curvature (C) derived from the mathematical models using RMSE and MAE; (2) the results of the qualitative estimation of the accuracy of C derived from the real-world DEM using the error of the spatial distribution. Figure 8 illustrates the DEMs generated from four ellipsoid and four Gauss surfaces with various parameters, and Table 3 shows the complexity parameters of E1, E2, E3, E4, G1, G2, G3, and G4. Table 4 illustrates the RMSEs and MAEs for the C calculated by the proposed algorithms under the different threshold values (30,50,100, and 150 m) on the E1 and G1.    From Table 4, we can see that the optimum threshold value for cutting the flow line over the FPN was 100 and 50 m on the ellipsoid and Gauss surfaces, respectively. The three comparison algorithms could not derive C on the boundary of the study area using a moving 3 × 3 window. Thus, Table 5 shows the RMSEs and MAEs for the C values calculated by the comparison algorithms and proposed algorithm under the optimal threshold on the eight surfaces, except for the error on the boundary.

Spatial Distribution of Residuals on the Mathematical Surface Models
From Table 5, we can see that the RMSEs and MAEs for the C values calculated by all four algorithms are constant, which is in accordance with the theoretical derivation. Figure 9 shows the spatial distribution of the residuals of the simulated C values relative to the theoretical C values on E1, E2, E3, and E4. Figures 10-13 illustrate the spatial distribution of the residuals of the simulated C values relative to the theoretical C values on G1, G2, G3, and G4, respectively. A positive value indicated that the simulated C value was greater than the theoretical C value, and a negative value indicated that the simulated C value was less than the theoretical C value. Table 6 shows the number of residual values within different levels of the four methods on E1, E2, E3, E4, G1, G2, G3, and G4. Table 7 shows the proportion of residual values within different levels of all four algorithms on E1, E2, E3, E4, G1, G2, G3, and G4. Table 8 shows the change of the proportion of the large and small errors for the new algorithm on E1, E2, E3, E4, G1, G2, G3, and G4.  Table 8 shows the change of the proportion of the large and small errors for the new algorithm on E1, E2, E3, E4, G1, G2, G3, and G4.    Table 8 shows the change of the proportion of the large and small errors for the new algorithm on E1, E2, E3, E4, G1, G2, G3, and G4.          subregion (number of the grid cells of 100 × 100) delineated by the red box in Figure 7, to support a visual comparison. The TFN and FPN of the subregion are illustrated in Figure 14. The C values of the presented method under different threshold values on the subregion are shown in Figure 15. From the figure, we can see the estimated C value under the cutting-up threshold of 200 m was the optimal value to be compared with other algorithms. The spatial distribution of the C values calculated by all four algorithms for the subregion are shown in Figure 16.

Discussion
Our new approach attempts to calculate the flow path curvature (C) values from the vector-based flow path network (FPN). The presented method converts the three-dimensional terrain into a one-dimensional flow line. The vector flow line can directly reflect the curve projected by the flow path over the horizontal surface and only utilizes the geometric parameter (coordinate of the break points and length of the line) to estimate the flow path curvature. Thus, it can avoid the local interpolation error based on the square-grid DEM and is the optimal choice of calculation algorithm, focusing on how the flow moves over the terrain surface. The experimental results demonstrate that the new approach has advantages over the selected comparison methods on the whole and can obtain a high accuracy with different terrains.

Accuracy Measures
The results of the quantitative test demonstrate that the C simulated by the FPN-based algorithm were generally closer to the theoretical value, and the algorithm was able to achieve a high accuracy for two mathematical surfaces. From Table 5, we can see that the RMSEs of the Evans, Zevenbergen, and Shary and proposed methods on E1, E2, E3, and E4 were 0.0012, 0.0024, 0.0012, and 0.0014 m, respectively. Compared to the Evans and Shary algorithms, the RMSE of the proposed method increased by 17% on E1, E2, E3, and E4. Compared to the Zevenbergen algorithm, the RMSE of the new method reduced by 42% on E1, E2, E3, and E4. The RMSEs of the Evans and Shary algorithms on G1, G2, G3, and G4 were 0.0066, 0.0070, 0.0087, and 0.0111 m, respectively. The RMSE of the Zevenbergen algorithm on G1, G2, G3, and G4 was 0.0040, 0.0050, 0.0073, and 0.0101 m, respectively. The RMSE of the proposed algorithm on G1, G2, G3, and G4 was 0.0043, 0.0052, 0.0076, and 0.0102 m, respectively. Compared to the Evans and Shary algorithms, the RMSE of the new algorithm reduced by 35%, 26%, 13%, and 8% on G1, G2, G3, and G4, respectively. Compared to the Zevenbergen algorithm, the RMSE of the proposed method increased by 8%, 4%, 4%, and 1% on G1, G2, G3, and G4, respectively. Moreover, the MAEs of the Evans, Zevenbergen, and Shary and proposed methods on E1, E2, E3, and E4 were 0.0002, 0.0011, 0.0002, and 0.0002 m, respectively. The MAE of the new approach was the same as that of the Evans and Shary algorithms on E1, E2, E3, and E4. Compared to the Zevenbergen algorithm, the MAE of the new approach reduced by 82% on E1, E2, E3, and E4. The MAEs of the Evans and Shary algorithms on G1, G2, G3, and G4 were 0.0029, 0.0041, 0.0059, and 0.0079 m, respectively. The MAE of the Zevenbergen algorithm on G1, G2, G3, and G4 was 0.0019, 0.0032, 0.0051, and 0.0072 m, respectively. The MAE of the proposed algorithm on G1, G2, G3, and G4 was 0.0025, 0.0038, 0.0058, and 0.0077 m, respectively. Compared to the Evans and Shary algorithms, the MAE of the new approach reduced by 14%, 7%, 2%, and 3% on G1, G2, G3, and G4, respectively. Compared to the Zevenbergen algorithm, the MAE of the new approach increased by 32%, 19%, 14%, and 7% on G1, G2, G3, and G4, respectively. These results show that the new approach can reduce the impact of landscape morphology on different terrain surfaces. Therefore, the new algorithm can generally obtain a relatively good result for the two above-mentioned surfaces.

Reasonable Spatial Distribution
The results of the spatial distributions of the residual values on the mathematical surfaces and of the estimated C values on the real-world terrain show the expected distribution patterns, and discrete patterns and some anomalous distributions. For E1, E2, E3, and E4, the results in Table 6 show that the number of the residual values between −∞ and −0.002 and between 0.002 and +∞ for the Evans and Shary algorithms was 563 and 80, respectively. The number of the residual values between −∞ and −0.002 and 0.002 and between 0.002 and +∞ for the Zevenbergen algorithm was 1288 and 804, respectively. The number of the residual values between −∞ and −0.002 and between 0.002 and +∞ for the new algorithm was 42 and 45, respectively. For the Evans and Shary algorithms, the number of the residual value between −∞ and −0.005 and between 0.005 and +∞ on G1 was 3492 and 3384, respectively. The number of the residual value between −∞ and −0.005 and between 0.005 and +∞ on G2 was 8105 and 4710, respectively. The number of the residual value between −∞ and −0.01 and between 0.01 and +∞ on G3 was 5660 and 2112, respectively. The number of the residual value between −∞ and −0.01 and between 0.01 and +∞ on G4 was 8309 and 4892, respectively.
For the Zevenbergen algorithm, the number of the residual value between −∞ and −0.005 and between 0.005 and +∞ on G1 was 1708 and 1334, respectively. The number of the residual value between −∞ and −0.005 and between 0.005 and +∞ on G2 was 6037 and 2363, respectively. The number of the residual value between −∞ and −0.01 and between 0.01 and +∞ on G3 was 5267 and 1554, respectively. The number of the residual value between −∞ and −0.01 and between 0.01 and +∞ on G4 was 7843 and 4477, respectively. For the new algorithm, the number of the residual value between −∞ and −0.005 and between 0.005 and +∞ on G1 was 1341 and 1471, respectively. The number of the residual value between −∞ and −0.005 and between 0.005 and +∞ on G2 was 2846 and 5603, respectively. The number of the residual value between −∞ and −0.01 and between 0.01 and +∞ on G3 was 4673 and 1908, respectively. The number of the residual value between −∞ and −0.01 and between 0.01 and +∞ on G4 was 7356 and 4551, respectively. These results demonstrate that the number of large errors of the new algorithm is generally lower than that of the other comparison algorithms.
From Table 7, we can see that the proportions of the large errors obtained by the Evans, Zevenbergen, and Shary and the proposed method were 4.28%, 13.92%, 4.28%, and 0.58%, respectively. The results in Table 8 show that the proportion of large errors of the new algorithm reduced by 86.45%, 95.83%, and 86.45% of that of the Evans, Zevenbergen, and Shary algorithms, respectively. For the Evans, Zevenbergen, and Shary and the proposed methods, the proportions of large errors on G1 were 17.19%, 7.61%, 17.19%, and 7.03%, respectively. The proportions of large errors on G2 were 32.04%, 21.00%, 32.04%, and 21.12%, respectively. The proportions of large errors on G3 were 19.43%, 17.05%, 19.43%, and 16.45%, respectively. The proportions of large errors on G4 were 33.00%, 30.80%, 33.00%, and 29.77%, respectively. The results in Table 8 show that the proportion of large errors of the new approach reduced by 59.10%, 34.08%, 15.34%, and 9.79% of that of the Evans and Shary algorithms on G1, G2, G3, and G4, respectively. The proportion of large error of the new approach reduced by 7.62%, −0.57%, 3.52%, and 3.34% of that of the Zevenbergen algorithm on G1, G2, G3, and G4, respectively. These results demonstrate that the proportion of large errors of the new algorithm is generally lower than that of other comparison algorithms.
For E1, E2, E3, and E4, the proportion of residual values of the new approach between −0.002 and 0.002 increased by 3.87% compared to the Evans and Shary approaches. The proportion of residual values of the new approach between −0.002 and 0.002 increased by 15.50% compared to the Zevenbergen algorithm. Combing the results in Figure 9, we can see that the residual value of the new approach was mainly distributed between −0.002 and 0.002, and was smaller than that of the comparison algorithms on E1, E2, E3, and E4. For the abovementioned order of Gauss surfaces, the proportion of residual values of the new approach between −0.001 and 0.001 increased by 21.03%, 16.29%, 16.06%, and 21.39% compared to the Evans and Shary approaches, respectively. Compared to the Zevenbergen algorithm, the proportion of residual values of the new approach between −0.001 and 0.001 reduced by 6.60%, 6.47%, 10.47%, and 9.29%, respectively. Combing the results in Figures 10-13, we can see that the spatial error of the proposed method was much higher than that of the Evans and Shary algorithms, and slightly lower than that of the Zevenbergen algorithm on G1, G2, G3, and G4 on the whole. Therefore, the performance of the new approach is generally better than that of the comparison algorithms on the ellipsoid and Gauss surfaces.
From Figures 14 and 16, we can see that the C value of the proposed algorithm was slightly higher than that of the Evans, Zevenbergen, and Shary algorithms on the ridges of extremely complicated terrain surfaces, which symbolized the area with obvious topographic relief. The new algorithm has a relative advantage over the selected comparison methods in gullies. Generally speaking, the FPN-based algorithm produces plausible outcomes, which conform to the real terrain. Because of the lack of field measurements of the partial derivatives, the quantitative evaluation of the proposed algorithm could not be simulated.

Consistent Simulation Results
From Table 4 Figure 14, the accuracy of the new algorithm was generally invariant with the change of the cutting-up threshold for the ellipsoid, Gauss, and real-world surfaces. Therefore, the proposed algorithm achieved the consistent simulation results.

Conclusions
In this paper, we propose a new approach to simulating the flow path curvature (C) using a vector-based method. This approach utilizes a new FPN-C method to derive C from the flow path network (FPN). The new algorithm aims to enhance the accuracy of simulating C and avoid interpolation errors, as well as being the choice of calculation algorithm for practical applications.
The presented method was implemented on the mathematical ellipsoid and Gauss surface models for a quantitative evaluation. Then, it was applied to a hilly plateau and mountainous region, located in the central Ganzi Tibetan Autonomous Prefecture of Sichuan Province, as a qualitative assessment. The results demonstrated that the new algorithm can obtain a relatively good result on different terrain surfaces. The new approach generally performed better than the comparison algorithms on the two above-mentioned surfaces. It was validated by both the quantitative evaluation (RMSE and MAE) and qualitative assessment (the visual comparison of the spatial distribution of the simulated C values on the mathematical and real-world surface). The RMSE and MAE of the new method were 0.0014 and 0.0002 m, reduced by up to 42% and 82% of that of the comparison algorithms on the ellipsoid surface, respectively. The RMSE and MAE of the presented method were 0.0043 and 0.0025 m at best, reduced by up to 35% and 14% of that of the comparison algorithms on the Gauss surface, respectively. For the ellipsoid surfaces, the residual value of the new approach was mainly distributed between −0.002 and 0.002 which was smaller than that of the comparison algorithms. The proportion of large errors of the new algorithm was 0.58%, reduced by up to 95.83% of that of the comparison algorithms. For the Gauss surfaces, the proportion of large errors of the new algorithm was 7.03% at best, reduced by up to 59.10% of that of the comparison algorithms. Moreover, this new approach can achieve consistent simulation results.
However, the vector-based flow line requires a high computing power in practical applications, especially in large DEMs. Therefore, the optimization and parallelization of this algorithm will be the focus of our future work. In addition, the B-spline interpolation method does not pass through the point used for interpolation. Thus, there may be a better algorithm to further enhance the accuracy when estimating the flow path curvature, which will be addressed in our future work as well. We will also continue to improve the calculation of the flow path curvature for the karstic or glacier type relief surfaces.