Math Geosci (2014) 46:445–481 DOI 10.1007/s11004-013-9508-8

Accuracy Analysis of Digital Elevation Model Relating to Spatial Resolution and Terrain Slope by Bilinear Interpolation Wenzhong Shi · Bin Wang · Yan Tian

Received: 28 March 2013 / Accepted: 7 November 2013 / Published online: 23 January 2014 © International Association for Mathematical Geosciences 2014

Abstract This study investigates the accuracy analysis of the digital elevation model (DEM) with respect to the following two major factors that strongly affect the interpolated accuracy: (1) spatial resolution of a DEM and (2) terrain slope. Unlike existing studies based mainly on a simulation approach, this research first provides an analytical approach in order to build the relationship between the interpolated DEM accuracy and its influencing factors. The bi-linear interpolation model was adopted to produce this analytic model formalized as inequalities. Then, our analytic models were verified and further rectified by means of experimental studies in order to derive a practical formula for estimating the DEM accuracy together with an optimization model for calculating the required resolution when a prescribed upper bound to the DEM accuracy is given. Moreover, this analytic approach can cope with either a grid-based DEM or a randomly scattered scenario whose efficacies have been validated by the experiments using both synthetic and realistic data sets. In particular, these findings first establish the rules for directly correlating the horizontal resolution of DEM data with vertical accuracy. Keywords DEM · Analytic inequality · Accuracy assessment · Terrain complexity · Bilinear interpolation method

W. Shi (B) · B. Wang (B) Department of Land Surveying and Geo-Informatics, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong e-mail: [email protected] B. Wang e-mail: [email protected] Y. Tian Department of Electronic and Information Engineering, Huazhong University of Science and Technology, Wuhan, China e-mail: [email protected]

123

446

Math Geosci (2014) 46:445–481

1 Introduction The digital elevation model (DEM) has been widely used in many applications, including digital terrain representation and geo-spatial analysis. A DEM is also used as source data to derive terrain features such as slopes, aspects, profiles and curvatures. The accuracy of a DEM is especially critical for situations where engineering accuracy is essential, such as hydraulic engineering applications. DEM accuracy can be affected by many factors, such as the interpolation model, the quality of the source data (acquisitions error, native scale, data type, etc.), the resolution (grid size) of the DEM, and degree of terrain complexity. Terrain feature is one of the prominent factors affecting soil erosion, and slope information is subsequently one of the important parameters for predicting soil erosion (Armstrong and Martz 2003; Wu et al. 2013). The scientific and effective representation of the relationship between terrain feature and hydrological landscape ingredients thus becomes the foundation of hydrology and soil erosion research. According to the geomorphologic principle, Hutchinson put forward a hydrologically correct DEM algorithm (Hutchinson 1988; Hutchinson and Gallant 1999, 2000) and used corresponding software ANUDEM for characterizing the hydrological landscape after exploring various interpolation methods and the data structure. Although many techniques are available to produce DEMs (Moore et al. 1991; Hutchinson 2011), ANUDEM has become one of the most well-known, reliable and computationally efficient tools for generating hydrologically sound DEMs from various data sources, including contour lines, spot heights and stream lines. (Yang et al. 2007) performed a study in the Coarse Sandy Hilly Catchments of the Loess Plateau, China, which demonstrated the effectiveness of a hydrologically correct digital elevation model algorithm for improving the DEM quality by identifying and correcting source topographic data errors and optimizing ANUDEM algorithm parameters. For this reason, ANUDEM has already been integrated into ArcGIS as one functional interpolation module invoked by Topogrid or Topo to Raster. A DEM surface is normally generated by applying interpolation models (bilinear, bicubic, B-spline, Kringing, etc.) based on a given set of node elevations. Several effective interpolation approaches have been reported (Pumar 1996; Caselles et al. 1998; Almansa et al. 2002; Shi and Tian 2006; Maune 2007). Of these, the bilinear interpolation method (Maune 2007) is the simplest, while other higher order interpolations (for example the Hermite method) may require more terrain information than just elevations, such as slope. Comprehensive studies have been conducted on evaluating DEM interpolation methods. Rees (2000) showed that simple bi-linear or bi-cubic convection was an adequate approach for DEM interpolation according to some case studies. However, based on different case studies, Kidner (2003) argued that high-order interpolation techniques were more accurate than the bilinear algorithm. Florinsky (2002) presented a new interpretation of the following three classes of digital terrain model (DTM) errors: (a) errors in interpolation of DEMs caused by the Gibbs phenomenon, (b) errors in DTM derivation from DEMs with enhanced resolution due to noise increase following DEM differentiation and (c) errors in DTM derivation caused by displacement of a DEM grid. Several preventative measures for handling these DEM errors have been recommended. Sinha and Schunck (1992) proposed a two-stage algorithm for discontinuity-preserving surface reconstruction using

123

Math Geosci (2014) 46:445–481

447

a weighted bicubic spline as a surface descriptor. By using zonal analysis, Shan et al. (2003) showed that the standard deviation of the 1-degree DEM can be largely approximated by a linear function of the slope change. Based on Grid-based Contours, Gousie and Franklin (2005) improved thin plate DEM approximations by the addition of tension into the biharmonic interpolation operator to counter the effects of Gibbs phenomena. The contour-based approach has also been employed for clods identification and characterization on a soil surface (Taconet et al. 2010). Because terrain factors are closely related to the applied terrain analyses and the accuracy of these factors is closely related to the accuracy of the original DEM, research on the relationships between DEM accuracy and terrain factors is very necessary. Extensive research on this subject has been reported in the literature (Holmes et al. 2000; Zhou and Liu 2004; Raaflaub and Collins 2006). In addition, the accuracy of a derived terrain slope is dependent on the accuracy of the original DEM, with larger errors associated with the areas of steepest slope (Toutin 2002; Wu et al. 2012). The tendency for slope accuracy to vary with grid size has also been reported (Chang and Tsai 1991; Gao 1997; Tang et al. 2001). Shi and Tian (2006) combined the bilinear and bi-cubic methods and proposed a hybrid interpolation method for DEM generation. His experimental results demonstrate that the hybrid method is effective for interpolating DEMs for various types of terrain. In those studies, it was assumed that the original data sources from which a DEM was generated were error-free and that the main error source affecting the accuracy of the interpolated DEM was the interpolation model error itself. Li and Zhu (2003) derived the mean error for a regular grid DEM interpolated by bilinear interpolation methods. Shi and Tian (2005) provided a model for estimating DEM propagation errors resulting from high-order interpolation algorithms, with the simplification that the node derivatives are replaced approximately by first order differentials. Zhu et al. (2005) derived the average DEM accuracy of the triangulated irregular network (TIN) model interpolated by linear methods with the assumption that the errors associated with the given nodes are independent of each other. (Kyriakidis and Goodchild 2006) further derived average line accuracy, the TIN model and a rectangle model, using bilinear interpolation methods where the errors at given nodes are supposed to be independent. In a practical study, Yamazaki et al. (2012) developed a new algorithm for adjusting a space-borne DEM in floodplain hydrodynamic modeling using drainage network information. A robust multiquadric method (Chen and Li 2013) was put forward to reduce the impact of outliers on the accuracy of DEM construction. In addition, the support vector machine (SVM) classifier is adopted for the adaptive selection of appropriate terrain modeling methods based on the terrain complexity (Jia et al. 2013). Generally, methods for characterizing those factors (for example, resolution and terrain slope) which affect DEM accuracy can be classified into two categories. The first is based on simulation while the second is analytical. Many studies have been performed using the simulation approach. Focusing on the TK-350 stereo-scenes of the Zonguldak test field in northwestern Turkey, Büyüksalih et al. (2004) showed that DEM accuracy depends mainly on the surface structure and slope of the local terrain. By analyzing different levels of terrain detail, Tang et al. (2001) built a linear model for DEM error estimation, DEM resolution and DEM mean profile curvature. By utilizing a quadratic polynomial model, Wang et al. (2004) provided a modified version of

123

448

Math Geosci (2014) 46:445–481

Tang’s model. Vaze et al. (2010) found that both DEM accuracy and resolution affect the quality of DEM-derived hydrological features as do the significant differences between the elevation and slope values derived from high resolution LiDAR DEM and their counterparts from coarse resolution contour derived DEM. However, these researchers conducted experimental studies without constructing analytical modeling. The principle of the simulation approach is straightforward and easy to understand, but the approach cannot provide predictive analytical expressions for identifying the relationships between DEM accuracy and the corresponding related factors. Furthermore, simulation results cannot be guaranteed as valid for all possible cases that might arise. The second approach for studying relationships between DEM accuracy and the above related factors is analytical in nature and based on theoretical analysis and mathematical inference. Theoretically, this approach is more rigorous, and the mathematical relationships of the analytical models are more generic in their applicability. For example, elevation data for floodplain mapping, published by the National Research (2007), investigated the relationship between DEM accuracy and terrain slope. However, the relationships between DEM accuracy and other related factors based on theoretical analyses have not received much attention. In addition, these analytic findings first establish the rules for directly correlating the horizontal resolution of DEM data with vertical accuracy, which has always been regarded as an unresolved issue in the DEM manual (Maune 2007). This research has focused on modeling DEM accuracy as a function of the DEM resolution and the average terrain slope. Unlike other studies, the research was mainly based on an analytical approach with the aim of deriving quantitative model relationships, expressed in the form of inequalities, and derived using mathematical theories. Inequality relationships are then further rectified to equality relationships on the basis of simulation experiments. The common bilinear interpolation method is adopted as a first step in studying the relationship between DEM accuracy and the resolution and terrain complexity. Accuracy analysis for nonlinear interpolation methods will be addressed in future studies. To facilitate a better understanding of the relationships to be studied, a one-dimensional case is first introduced and then our analysis is extended into two dimensions, which are of more concern. The rest of this paper is organized as follows. In Sect. 2, the principle of the bilinear model is briefly introduced. In Sect. 3, theoretical findings of DEM accuracy in relation to DEM resolution and terrain slope factors are concluded for both one-dimensional and two-dimensional cases for raster and scattered DEM data. Sections 4.1–4.3 are devoted to conducting experimental studies based on both real DEM data and synthetic data aimed at verifying the validity of the analytical findings obtained in Sect. 3. Meanwhile, Sect. 4.4 further derives a more accurate quantitative model of the relationship between the sampling ratio and average terrain slope by fitting of the DEM RMSE error and its corresponding upper bound based on the experimental analysis on Gaussian surfaces. Conclusions drawn from this research are presented in the last section. 2 Bilinear Interpolation Model The method applied for DEM interpolation is a significant factor affecting the accuracy of the DEM. In general, DEM interpolation methods can be classified as linear and non-

123

Math Geosci (2014) 46:445–481

449

linear. Each of these interpolation methods has its own characteristics. For instance, linear interpolation is applicable for areas with lower levels of terrain change and is easily implemented, while nonlinear interpolation seems to be more appropriate for higher levels of terrain change. However, this latter method is normally timeconsuming in implementation (Shi and Tian 2006). In this research, a bilinear linear interpolation model was used as an example to study the impact of DEM resolution and terrain slope on DEM accuracy. In the following equations, the bilinear interpolation method is briefly introduced. The bilinear polynomial interpolation model for a regular grid takes the following form z = a1 + a2 x + a3 y + a4 x y,

(1)

where ai (i = 1, 2, 3, 4) are coefficients that need to be determined. For the onedimensional case, the above equation is simplified to z = a1 + a2 x,

(2)

where ai (i = 1, 2) are coefficients that need to be estimated in advance. There are two ways of implementing the bilinear interpolation model. One is to provide four points in advance and substitute these points into the bilinear formula to obtain a group of four equations. Their solutions give the values of the four parameters ai (i = 1, 2, 3, 4). These linear equations can be solved using Cramer’s rule. The second method is to divide the bilinear interpolation into three one-dimensional linear interpolations. Stated more precisely, the first two linear interpolations are carried out along the vertical direction of the original cell. Based on these interpolations, the third linear interpolation is then implemented along the horizontal direction of the original cell passing through the point in question. This type of bilinear interpolation is frequently adopted due to its simplicity and practicality. 3 Analytical Modeling of DEM Accuracy In this section, we analyze the impact of DEM resolution and terrain slope on the interpolated DEM accuracy. The result of this analysis is an analytical model depicting the relationships between the variables involved. To facilitate an understanding of the derivation process, our analysis started with the one-dimensional case and proceeded to the two-dimensional case. The analysis procedure for both cases is as follows. First, the original DEM with higher resolution is provided. Then, a re-sampled DEM with lower resolution is generated on the basis of the original DEM by a down-sampling method. Next, the interpolated DEM is obtained by applying the linear interpolation model (for the one-dimensional case) or the bilinear interpolation model (for the two-dimensional case) on the re-sampled DEM. Finally, the relationships between the accuracy of the interpolated DEM and (a) DEM resolution and (b) terrain slope are modeled. 3.1 Accuracy Analysis for One-Dimensional, Uniformly Distributed Points For a general re-sampling ratio r in the one-dimensional case, the original DEM and the corresponding re-sampled DEM are denoted by f (ri + j) and f¯(ri + j), respectively.

123

450

Math Geosci (2014) 46:445–481

As illustrated in Fig. 1, points represented as , such as x0 , x4 , x8 , x12 , denote the horizontal coordinates of the sampling data, while the circled stars , such as x1 , x2 , x3 , mark the coordinates of the true elevation values. The notation is explained as follows: n is the total number of terrain elevation items of the data, i denotes the region number, such as [x0 , x4 ], to be interpolated, and j is the point index located within each region. Therefore, xk = xri+ j ; 0 ≤ i < n/r, i ∈ N ; 0 ≤ j ≤ r, j ∈ N . The linear interpolation method for the one-dimensional case leads to the following equations   j j ¯ f (ri) + f (ri + r ) . (3) f (ri + j) = 1 − r r Hence, the difference between the true elevation value and the interpolated one is    j j ¯ f (ri + j) − f (ri + j) = f (ri + j) − 1 − f (ri) + f (ri + r ) r r   j = 1− ( f (ri + j) − f (ri)) r j + ( f (ri + j) − f (ri + r )) . (4) r Therefore, the root mean square error (RMSE) between the original true elevation and the interpolated DEM values can be calculated as (RMSE)2 ≤

  2 r4 − 1 · s 2 E{(Slope)2 }. 15r 4

(5)

See Appendix A (part 1) for the derivation of Eq. (5) in detail. In terms of this finding, if the terrain slope can be estimated beforehand, then the upper bound of the square of RMSE by bilinear interpolation is directly proportional to the square of the re-sampled DEM resolution (re-sampling DEM post spacing) s 2 and inversely proportional to the average value of the square of the DEM slope E{(Slope)2 } with a scale multiple of f

x0 x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 i-1

i

xn

x

i+1

Fig. 1 Schematic diagram of the linear interpolation of the DEM with re-sampling ratio r = 4 for onedimensional case

123

Math Geosci (2014) 46:445–481

451

2(r 4 −1) . This conclusion is consistent with our intuitive observation of DEM accuracy 15r 4

in relation to the slope and resolution factors. Therefore, a mathematical model (the inequality of Eq. 5) has been provided to quantify the relationship between DEM accuracy and (a) DEM resolution and (b) terrain slope. 3.2 Accuracy Analysis for Two-Dimensional Grid-Based DEM For the general re-sampling ratio r in the two-dimensional case, the original DEM and the corresponding re-sampled DEM are denoted by f (ri + k, r j + l) and f¯(ri + k, r j + l), respectively. As illustrated in Fig. 2, the points represented as , such as (x0 , y0 ), (x0 , y4 ), (x4 , y0 ) and (x4 , y4 ), denote the projection coordinates of the sampling data, while the points represented by , such as (x1 , y1 ), are those coordinate points of the true elevation values. The notation is as follows: m is the total number of data points in the x direction while n is the total number of data points in the y direction. Therefore, mn is the total number of points in this area. Meanwhile, i denotes the region number of points, such as [x0 , x4 ], to be interpolated along the x direction, while j denotes the region number of points, such as [y0 , y4 ], to be interpolated along the y direction, k is the point index located within each region i and l is the point index located within each region j. Additionally, the original DEM cell is assumed to be square, that is, the cell size is h along both the x and y directions. Therefore, 0 ≤ i < m/r, i ∈ N ; 0 ≤ k ≤ r, k ∈ N , 0 ≤ j < n/r, j ∈ N ; 0 ≤ l ≤ r, l ∈ N . It is easy to obtain the following equations for the one dimensional linear interpolation case

y3

y4

y

f

y0

y1

y2

j

x0

x1

x2

x3

x4

x

i Fig. 2 Schematic diagram of the linear interpolation DEM with re-sampling ratio r = 4 for twodimensional case

123

452

Math Geosci (2014) 46:445–481

     l k k 1− f (ri, r j) + 1− f¯ (ri + k, r j + l) = 1 − r r r       k l k l + 1− f (ri, r j + r )+ r r r r

l r

 f (ri + r, r j)

f (ri + r, r j + r ) . (6)

Similarly, but by a more complex formula deduction process, the mean square error (RMSE)2 can be calculated as follows (a more rigorous derivation can be found in Appendix A (part 2) of this paper)    2 8r 2 − 9r + 4 r 4 − 1 h 2 E{(Slope)2 } (RMSE) ≤ 4 45r    2 8r 2 − 9r + 4 r 4 − 1 s 2 = E{(Slope)2 }. 45r 6 2

That is to say    2 8r 2 − 9r + 4 r 4 − 1 · s 2 E{(Slope)2 }. (RMSE) ≤ 45r 6 2

(7)

A more accurate model of the relationship between sampling ratio and average terrain slope is obtained by adjustment of the analytical model presented as Eq. (13) based on the experimental analysis of Gaussian surfaces given below in Sect. 4.3. Note that for the two-dimensional case, the slope at point (i, j) is defined as ( f x2 + f y2 )1/2 , which is the Euclidean norm for the gradient of the terrain function f at location (i, j). Therefore, the average slope squared for the entire terrain surface should be defined as E{(Slope)2 } =

1  2 2 f (i, j) + f (i, j) . x y n2

(8)

i, j

Based on Eq. (7), it can be concluded that the smaller the slope value or the higher resolution of the original DEM (i.e., the smaller of s value), the corresponding RMSE of the bilinear interpolated DEM becomes smaller. 3.3 Accuracy Analysis for Two-Dimensional Randomly Scattered DEM In this section, consideration has shifted to the two-dimensional, randomly scattered DEM data. Connection of the scattered data points into triangles forms a triangulated irregular network (TIN) rather than a regular square grid. Thus, a straightforward numerical difference method cannot be applied for estimating the information of the topographic slope. Therefore, one solution for this dilemma is the conversion of the randomly scattered DEM elevation data into the regular square grid data. The similar result presented in Eq. (7) certainly remains valid after such conversion because the global DEM interpolation error indicated by RMSE is inherently determined by the

123

Math Geosci (2014) 46:445–481

453

Fig. 3 Convex combination of the triangular derivatives around node P in the triangulation

π2

π1

P π5

π3 π4

specific topographic character, interpolation method and resolution size. However, such a new directly transformed DEM may inevitably contain additional errors, for instance, arranging extra temporary sampling locations when converting TIN to new mesh grids, as well as further numerical differences for estimating the partial derivatives along the x and y directions. To this end, we first adopt one method of local derivative estimation by using a convex combination of all derivatives on related triangular planes Goodman et al. (1995). Pick up a node P from the triangulation of the scattered DEM data and let πi , i = 1, . . . , k be the triangles taking P as its one vertex as illustrated in Fig. 3. We denote gi as the gradient of the linear interpolation to the DEM data at vertices of πi . One method for estimating the gradient at P is given by a convex combination of the gradients of its located triangles, g1 , . . . , gk , that is DP =

k  i=1

λi gi /

k 

λi ,

(9)

i=1

where the weights λi are selected to be the reciprocal of the base triangle’s area i of the corresponding triangle πi , that is, λi = 11/i . In addition, the gradient gi and the area i can be obtained according to the following computational procedure. Suppose that three vertices (x j , y j ) with corresponding DEM elevation values z j , j = 1, 2, 3, produce a triangle plane πi given by ax + by + cz + d = 0, where (a, b, c)T constitutes the plane’s normal vector a = y1 (z 3 − z 2 ) + y2 (z 1 − z 3 ) + y3 (z 2 − z 1 ) b = x1 (z 2 − z 3 ) + x2 (z 3 − z 1 ) + x3 (z 1 − z 2 ) c = x1 (y3 − y2 ) + x2 (y1 − y3 ) + x3 (y2 − y1 ) .

123

454

Math Geosci (2014) 46:445–481

Therefore, the gradient for this triangle plane πi is  gi =

∂z ∂z , ∂x ∂y

T

  b T a = − ,− , c c

and the corresponding base area is i = 21 |c|. As mentioned above, Eq. (7) remains valid for the newly generated regular square grid from the randomly scatted DEM elevation data. Meanwhile, the global average of the terrain slope square (indicator for estimating the degree of DEM rugged bumpiness) can be approximated by averaging all local derivative estimations by means of the above generalized numerical finite difference calculation method E{(Slope)2 } =

n 1 Di 22 , n

(10)

i=1

where n is the number of the total nodes constituting TIN and Di represents the convex combination of the gradients at node i according  to Eq. (9). In addition, the square of the resolution size is adjusted to be s 2 = n1 i∈T i , where T is the triangulation coverage of the concerned district, based on the fact that the area cannot be changed by such conversion from TIN into Rectangle. Therefore, the formula for characterizing the relationship between the RMSE and the sampling ratio, the terrain slope for TIN, becomes (RMSE)2 ≤

   n 2 8r 2 − 9r + 4 r 4 − 1 1  1 Di 22 . ·  · i 45r 6 n n i∈T

(11)

i=1

Suppose that there are N original √randomly scattered DEM data and n sampling points making the sampling ratio r = n/N . It should be noted that here the sampling ratio is supposed to represent the reciprocal of the sampling proportion in terms of the number of DEM points. However, introducing the Sqrt operation is nothing less than maintaining the consistency with the previous sampling ratio as employed in Sects. 3.1 and 3.2 referring to the sampling proportion in terms of the one-dimensional interval length. Indeed, it may still introduce some tiny DEM errors after such conversion from scattered data into the regular values. However, at the outset, our manuscript (Sect. 3.2) only provides the initial inequality relationships between the RMSE indicator (DEM interpolation error estimation) and the resolution and terrain slope. Apparently, these intrinsically qualitative relationships between the DEM error and its influencing factors (resolution and topographic relief represented by the terrain slope) would not be affected after such conversion from TIN to Rectangle. Meanwhile, more precise quantitative relationships can be rectified in Sect. 4.4 by experimental study of a specified geomorphic district because these estimated relationships between RMSE, resolution and terrain slope most likely vary with different resolutions and rugged bumpiness of the various landforms.

123

Math Geosci (2014) 46:445–481

455

4 Experimental Analysis 4.1 Experimental Analysis for Realistic Grid-Based DEM Data Sets To verify the validity of the findings and further adjust the results presented by Eq. (7), which depict the analytical relationships between DEM accuracy and (a) the DEM resolution and (b) terrain slope, an experimental study was carried out to validate the theoretical findings using practical DEM data sets with different levels of complexity. The three test areas are located within Shannxi province (left part of Fig. 4), and the three test areas/sample regions represent the three typical types of terrain: plain, hill and mountain. Sample region A lies within 33.577◦ –33.844◦ N and 107.465◦ –107.785◦ E, situated in the Qianlong Mountains, in the southern part of Shanxi province. Sample region B is located within 34.231◦ –34.497◦ N and 108.135◦ –108.455◦ E and is situated on the central Shaanxi plain. Sample region C is located within 36.065◦ –36.322◦ N and 109.487◦ –109.804◦ E. It is a loess hill and gully area in northern Shannxi, a hilly topography. All three sample regions are illustrated in Fig. 4 and their characteristics are described in Table 1.

Fig. 4 Three different test areas within Shannxi province representing mountain, hill and plain terrains

123

456 Table 1 Characteristics of the topographic surfaces of the studied area

All surfaces are grid DEMs composed of 106 points and with 30×30 m spacing

Math Geosci (2014) 46:445–481 Terrain descriptive statistics

High mountain

Hill

Plain

Min (m)

1,067

1,005

224

Max (m)

3,069

1,464

401

Average elevation (m)

2,002

1,266

338

Hmax − Hmin (m)

2,002

459

177

Std dev. of elevation (m)

386.70

78.28

11.26

Elevation coefficient of variation (%) Min slope (◦ )

19.32

15.98

2.95

0

0

0

Max slope (◦ )

69.75

61.46

53.25

Average slope (◦ )

26.38

17.78

7.10

Std dev. of slope (◦ )

10.89

8.46

4.56

Slope coefficient of variation (%)

41.28

27.58

14.23

The newly released DEM data, advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model Version 2 (GDEM2), produced by The Ministry of Economy, Trade, and Industry of Japan and the United States National Aeronautics and Space Administration (NASA), were used for the study. Because of the 260,000 additional scenes, GDEM2 is reported to have a much higher spatial resolution than its previous version, GDEM1. Moreover, the negative 5m overall bias detected in GDEM 1 has been removed. GDEM2 is in GeoTIFF format with geographic lat/long coordinates and a 1 arc-second (30 m) grid of elevation postings. GDEM2 is referenced to the WGS84/EGM96 geoid. Pre-production estimated accuracies for this global product were 20 m at 95% confidence levels for vertical data and 30 m at 95% confidence levels for horizontal data. When calculating the terrain gradient value, the projection was converted to UTM/WGS84. Three topographic surfaces measuring 1,000 by 1,000 pixels were extracted from the original dataset to be the test data. Based on our theoretical findings above, for the experiment a quarter of the original terrain data were first evenly selected (that is to say, a sampling ratio of r = 2 in this experiment), and then interpolated by bilinear interpolation. Therefore, the RMSE can be calculated by comparison of the elevation differences between these two DEM surfaces. Meanwhile, slopes for every point were extracted based on the octagonal linear differences. Figure 5 displays the original DEM and the bilinear interpolated DEM, as well as the corresponding DEM residuals of a 60 × 60-point area for the top left corner of the three plain, hill and mountain districts. It is difficult to obtain a perfect picture if all 1000 × 1000 points are presented together in a single image for each terrain case. The bilinear interpolation method was applied on the re-sampled DEM and the interpolated DEM accuracy was assessed by comparing the difference between the original DEM data and the interpolated data, that is, | f (x, y) − f (x, y)|, for all the points on the DEM. A striking visual contrast was obtained by the above experiment as the terrain complexity increased significantly in each of the three cases. Table 2 below presents the corresponding results for this experimental study.

123

6

3.718

3.7175

3.991

3.7905

3.9915

3.791

3.9905

3.79

3.645

7.29

2.37

7.295

3.65

2.375

7.3

2.38

5

5

x 10

x 10

3.655

5

x 10

2.385 6

x 10

6

3.7185

1600

1800

2000

6

x 10

3.992

1400 1350 1300 1250

x 10

600 550

3.718

3.9915

3.791

3.7175

3.991

3.7905

3.9905

3.79

Fig. 5 Comparative visualization of the original and interpolated terrain surface for mountain, hill and plain areas within Shannxi province

x 10

6

3.7185

1600

1800

2000

x 10

6

3.992

1400 1350 1300 1250

x 10

600 550

7.29

3.645

2.37

7.295

3.65

2.375

7.3

5

5

x 10

5

x 10

2.385

x 10

3.655

2.38

Math Geosci (2014) 46:445–481 457

123

458

Math Geosci (2014) 46:445–481

Table 2 The experimental results for three different terrain types E{(Slope)2 }

s0 = 2h 0

Terrain type

(RMSE)2

3/16s 2 E{(Slope)2 }

Plain

57.4824

0.0227

9.0416

14.0885

Hill

56.2585

0.1364

15.2410

80.9341

Mountain

57.7782

0.3372

25.9608

211.0785

Table 3 The experimental results for rectangle with different levels of terrain complexity (r = 2) Terrain parameters A

B

s0 = 2h 0

E{(Slope)2 }

(RMSE)2

3/16s 2 E{(Slope)2 }

C

9

30

1

300

0.0005

0.8208

8.6449

50

150

9

300

0.0512

82.0424

863.9056

200

500

25

300

0.1538

500

1,500

50

300

1.3087

249.4699 2109.714

2594.7536 22084.4387

3 2 For the plain region, the values of E{(Slope)2 }, {RMSE}2 and 16 s E{(Slope)2 } are 3 2 0.0227, 9.0416 and 14.0885, respectively. The value 16 s E{(Slope)2 } was calculated for the re-sampling DEM with a resolution twice that of the original cell size. In this case, 9.0416 < 14.0885; that is to say inequality of Eq. (7) and (RMSE)2 < 3 2 2 16 s E{(Slope) } is valid according to the experiment. In fact, similar experimental findings were obtained for the hill and mountain cases. As observed from Table 3, the levels of terrain complexity steadily increase from the top to the bottom in Fig. 5, as do 3 2 s E{(Slope)2 }. their corresponding RMSEs and their upper bounds in the column 16 This result nicely demonstrates the validity of the findings represented by Eq. (7). It 3 2 s E{(Slope)2 } values are often larger than the RMSEs should be noted that these 16 for the introduction of inequality techniques and the approximation of terrain slope with its mean value for estimating the DEM accuracies.

4.2 Experimental Analysis for Synthetic Grid-Based DEM Data Sets This section is devoted to verifying the validity of the proposed model with respect to different levels of simulated terrain complexity using mathematical surfaces. Actually, the terrain complexity is difficult to fully quantify due to the many factors involved, such as the frequency of changes in terrain surface orientation, the average slope and the amplitude of the surface. Below only the surface amplitude considered for characterizing the terrain complexity, as also adopted by (Zhou and Liu 2002, 2003) in generating terrain surfaces of different levels of complexity. The Gaussian surface function is expressed by  x 3 y 5  y 2 x 2 x x 2 −( x )2 −( y +1)2 n e−( m ) −( n ) e m −B − z = A 1− − m 5m m n 2 y 2 x − Ce−( m +1) −( n ) ,

123

(12)

Math Geosci (2014) 46:445–481

459

Fig. 6 Comparative visualization of the original and interpolated Gaussian surfaces with different complexity levels (rectangle)

where A, B and C are parameters determining the amplitude and m and n are the parameters which control the spatial extent of the surface. The true values of the elevation and slope can be computed using the above expression. For the sake of convenience, m = n = 500 was set, meaning the terrain was located within x, y ∈ [−1,000, 1,000]. In this study, four different levels of terrain complexity were generated by assigning the Gaussian surface function with four groups of parameter values. In addition, the sampling ratio r is still set at 2 for this experimental study. The visualization of the interpolation results is shown in Fig. 6, including a comparison between the original and the interpolated surfaces together with the residuals. As observed in the above table, the levels of terrain complexity increase in turn from the top to the bottom in Fig. 6, as do their corresponding RMSEs and 3 2 s E{(Slope)2 }. This result fully demonstrates the their upper bounds in the column 16 3 2 s E{(Slope)2 } values are often much validity of Eq. (7). It should be noted that the 16 larger than the RMSEs for the introduction of inequality techniques and the approximation of terrain slope with its mean value for estimating the DEM accuracies. More experiments with real data sets are presented in the next section expressing DEM RMSE values in terms of both the average terrain slope and the re-sampling ratio.

123

460

Math Geosci (2014) 46:445–481

4.3 Experimental Analysis for Randomly Scattered DEM Based on Mathematical Surfaces This section is devoted to verifying the theoretical findings from Sect. 3.2. The Gaussian surface together with its configuration is adopted again for carrying out experiments on randomly scattered DEM data. Still, the terrain district was restricted within x, y ∈ [1,000, 1,000] and the same four groups of parameter values generate four different Gaussian surface topographies. Figure 7 illustrates the comparison of the original image and the bilinear interpolation image, together with their interpolated residuals placed in the middle column. For each surface, 400 locations are randomly selected out from the above experimental district to create a triangulation mesh, among which only 100 locations and their elevations are randomly chosen as the seed points for spanning the interpolation surface, thereby approximating the elevations of the other 300 locations. It should be noted that only these 100 DEM seed points as well as their elevations can be adopted to calculate the upper bound of the DEM error, the right term in Eq. (11), while in practical applications usually no surplus check points are available for exam-

Fig. 7 Comparative visualization of the original and interpolated Gaussian surfaces with different complexity levels (TIN)

123

Math Geosci (2014) 46:445–481

461

Table 4 The experimental results for TIN with different levels of terrain complexity (r = 2, 4) r=

Terrain parameters A 9 50 200 500

B

√

n/N

(RMSE)2

C 30 150 500

1,500

1 9 25 50

2(8r 2 −9r +4)(r 4 −1) 6  45r n Di 22 · n1 i∈T i · n1 i=1

Multiple

2

0.8740

4.8818

5.586

4

0.8378

6.6595

7.949

2

22.299

128.17

5.748

4

21.617

171.64

7.940 6.042

2

257.10

1553.4

4

253.38

2020.2

7.973

2

2222.1

12697

5.714

4

2148.7

17043

7.932

ining the real DEM interpolation error. Although the theoretical findings in Sect. 3.2 cannot characterize the realistic DEM error in terms of RMSE, they still provide the upper bound estimation as a function of the sampling ratio, resolution size and terrain slope. Table 4 lists the experimental results corresponding to four different levels of terrain complexity with the sampling ratio at 2 and 4. From this table, it can be observed that the upper error bounds fully confine the square of the RMSE. Moreover, an apparent finding emerges that a steady constant multiple exists between the square of the RMSE of the DEM error and its upper bound, which is related to the sampling ratio, but is independent of terrain complexity. The following subsection further extends the discussion of these findings for the scenario of regular square grids. 4.4 Extension of the Study and Analysis of the General Re-Sampling Ratio Case Based on Gaussian Surfaces Functionalities of the experimental studies in Sects. 4.1 and 4.2 are mainly to validate the theoretical findings using a re-sampling ratio of 2 for regular mesh grids. In this section, attention is turned to the general re-sampling ratio r . Furthermore, as mentioned above, a large gap exists between the DEM interpolation error in terms of the RMSE and its upper bound as formulated by Eqs. (5) and (7). To this end, the following subsection is devoted to narrowing this gap by the data fitting method. 17 groups of terrain parameters for the Gaussian surface were selected for similar experiments (see Tables 5, 8, 9 in Appendix B) to obtain a more precise relationship description between DEM accuracy (RMSE) and its upper bound with re-sampling ratios from 2 to 4 and even larger. Because the inequality of Eq. (7) leads to a large gap between the RMSE and its upper bound, the RMSE upper bound does qualitatively account for the DEM accuracy. However, by conducting a large number of numerical simulation experiments, we find that it is reasonable to assume a linear relationship between the DEM accuracy expressed by the (RMSE)2 indicator and its upper bound 2(8r 2 −9r +4)(r 4 −1) 2 s E{(Slope)2 } consisting of the re-sampling ratio and the average 45r 6

123

462

Math Geosci (2014) 46:445–481

Table 5 The experimental results for different levels of terrain complexity with r = 2, k = 10.4125 Terrain parameters A

B

3

s0 = 2h 0

E{(Slope)2 }

300

0.0001

(RMSE)2

3/16s 2 E{(Slope)2 }

C 10

1/3

0.0912

0.9605

9

30

1

300

0.0005

0.8208

8.6449

30

100

3

300

0.0057

9.1158

95.9895

300

750

30

300

0.0228

36.4633

383.9581

50

150

9

300

0.0512

82.0424

863.9056

60

200

6

300

0.0910

145.8531

1535.8322

90

300

9

300

0.1427

228.5045

2408.8325

120

400

12

300

0.3454

560.5679

5829.2307

400

900

45

300

0.5194

843.6800

8765.2211

150

500

24

300

0.0132

21.1793

222.0079

200

500

25

300

0.1538

249.4699

2594.7536

450

1,200

50

300

0.8658

1402.4979

14611.0981 22084.4387

500

1,500

50

300

1.3087

2109.7141

600

1,600

60

300

1.5382

2491.9384

25957.5923

700

1,800

80

300

1.9719

3197.3737

33275.0838

800

2,000

90

300

2.4583

3988.8640

41483.6484

1,000

2,000

100

300

2.7194

4407.5113

45889.5424

terrain slope. In other words, the analytical inequality (Eq. 7) does not quantify DEM vertical errors (by RMSE2 ) with an explicit function of the terrain slope and DEM horizontal resolution. However, it has actually clarified this issue with only one multiple to be determined. According to Eq. (7)    2 8r 2 − 9r + 4 r 4 − 1 · s 2 E{(Slope)2 } = k (RMSE)2 . 45r 6

(13)

For each fixed re-sampling ratio r , a linear coefficient k can be calculated based on the experimental results for the 17 terrain types by means of the least squares method. For example, the calculated result for a re-sampling ratio of 2:1 is presented below in Table 5. Figure 8 is the corresponding fitted relationship line represented by Eq. (13). Two other comparison figures are provided in Appendix B as Figs. 10 and 11, corresponding to Tables 8 and 9 for each re-sampling ratio from Eqs. 3 and 4, respectively. The experiments show that the parameter k becomes ever larger with the increase of the re-sampling ratio r . Further experiments were then conducted to more explicitly establish their relationship. For this reason, all re-sampling ratios below 30 were tested, leading to a series of linear coefficients k for different re-sampling ratios as shown in Table 6.

123

Math Geosci (2014) 46:445–481

463

Fig. 8 The fitted relationship between {RMSE}2 and its upper bound with r = 2, k = 10.4125 Table 6 The experimental coefficient k corresponding to each re-sampling ratio r r

2

3

4

5

6

7

8

9

10

11

k

10.41

12.40

13.55

14.30

14.82

15.20

15.50

15.74

15.93

16.09

r

12

13

14

15

16

17

18

19

20

21

k

16.23

16.34

16.44

16.53

16.61

16.68

16.74

16.79

16.84

16.88

r

22

23

24

25

26

27

28

29

30

k

16.92

16.96

17.00

17.03

17.06

17.08

17.11

17.13

17.16

The following expression can be the optimized curve providing the relationship between the adjusted coefficients k and each corresponding re-sampling ratio r kˆ = 6.48 ln(r ) − 2.77 × 10−4 r 3 + 2.30 × 10−2 r 2 − 0.86r + 7.59.

(14)

The curve of Fig. 9 demonstrates visually the closeness of fit. Thus far, by combining Eqs. (13) and (14), the exact relationship between DEM accuracy and both re-sampling ratio (resolution) and average terrain slope for a DEM that is generated by bilinear interpolation can be given as follows    1 2 8r 2 − 9r + 4 r 4 − 1 (RMSE) = · · s 2 E{(Slope)2 }. 45r 6 kˆ 2

(15)

By deformation of Eq. (15), this model can also be adopted to determine the necessary DEM resolution for practical applications when a prescribed RMSEp value is given supposing the fluctuation of the concerned terrain indicated by the average squared

123

464

Math Geosci (2014) 46:445–481

Fig. 9 The fitted relationship between the linear coefficients k and re-sampling ratio r

slope has been approximated. To be exact, the precise re-sampling ratio r can be obtained by solving its implicit function (Eq. 15). In other words, a suitable re-sampling ratio r should be determined by means of Eq. (15); it must be an integer and also a larger ratio r corresponding to a lower resolution because s = r h leads to higher DEM errors in terms of (RMSE)2 . What concerns us, therefore, is the maximum re-sampling ratio r , which can be expressed as an integer optimization problem as follows max r RMSE  RMSEp , s.t. r∈N

(16)

in which RMSE is calculated according to Eq. (15). The optimized re-sampling ratio can be found by computing several re-sampling ratios below a large prescribed integer value, such as 20, as an approach to solving the above optimization problem. This finding should have great potential in DEM data generation and applications. For a specified topographic region, its average value of the global terrain slope would be certain and can be calculated by numerical finite differences of all DEM data or partial samplings. In this case, the suitable re-sampling ratio r can be calculated in terms of Eq. (16) in order to meet the prescribed accuracy requirement in terms of RMSEp . Actually, this formula has often been urgently required in many practical engineering DEM applications. The result given by Eq. (15) first establishes the rules for directly correlating the horizontal resolution of DEM data with vertical accuracy, which has always been regarded as an unresolved issue in the DEM manual (Maune 2007). As a response to this issue, Eq. (15) directly correlates the horizontal resolution sˆ of the digital elevation data with the vertical accuracy RMSEp . Thus, it could be tabulated for comparison and reference with regard to the common re-sampling ratios for the requirement of practical implementation shown below in Table 7. The second column indicates the

123

Math Geosci (2014) 46:445–481

465

Table 7 The practical DEM accuracy indicator and resolution with sampling ratio 2 ≤ r ≤ 5 Sampling ratio r

DEM accuracy (RMSE)2

Re-sampling ratio rˆ

2

(RMSE)2 = 0.0180s 2 E{(Slope)2 }

3

(RMSE)2 = 0.0193s 2 E{(Slope)2 }

max r RMSE ≤ RMSEp s.t. r∈N

4

(RMSE)2 = 0.0196s 2 E{(Slope)2 }

5

(RMSE)2 = 0.0197s 2 E{(Slope)2 }

DEM accuracy (RMSE)2 interpolated by the bilinear method when the resolution size and the average terrain slope are given. The last column lists the required re-sampling ratio rˆ when the prescribed RMSEp and the average terrain slope are given.

5 Conclusions In this research, the relationship between the accuracy of an interpolated DEM and (a) the DEM resolution and (b) slope of terrain has been modeled analytically from the theoretical analysis. For both one-dimensional and two-dimensional DEMs, relationship models were built and presented first in the form of two inequalities. These relationship models reveal the fact that the accuracy of an interpolated DEM depends on (a) its resolution and (b) the slope of the terrain, an indication of terrain complexity level. When the slope value is smaller and the DEM resolution is greater, the generated RMSE for the bilinear interpolated DEM will become smaller. Two types of experimental study were conducted based on both real DEM data sets and mathematically generated DEM surfaces. The experimental results demonstrate that the inequality expressions developed for the relationships between DEM accuracy, DEM resolution and terrain slope are valid. Experimental studies verified the validity of the analytical relationship models. Furthermore, a more precise expression of DEM accuracy in terms of RMSE with respect to both DEM resolution and average terrain slope was obtained. In addition, a practical optimized model was derived for calculating the required re-sampling ratio r when a prescribed DEM accuracy upper bound in RMSE is given. The proposed theoretical models providing the relationships between the interpolated DEM accuracy and its resolution and terrain slope can be used for determining DEM resolution in practical DEM applications. This new knowledge might be included in future DEM manuals for determining DEM resolution or to estimate DEM accuracy in engineering applications involving DEM. For a given terrain with a fixed slope and DEM accuracy level, the DEM resolution needed can be deduced using the relationship model developed in this study. On the other hand, if both DEM resolution and terrain slope are known, an accuracy estimation for the interpolated DEM can be provided in terms of the maximum possible RMSE for the interpolated DEM using the developed relationship models, that is, the inequalities of Eq. (7), together with the more precise estimated expression using Eq. (15).

123

466

Math Geosci (2014) 46:445–481

As mentioned in the DEM manual (Maune 2007), no established criterion exists to depict the vertical accuracy of digital elevation data by its horizontal resolution. The work presented here can be viewed as an attempt to provide a solution to this issue. Indeed, our findings do directly correlate the horizontal resolution s = r h of the digital elevation data with the vertical accuracy RMSE as expressed by Eq. (15). Future development of this study will investigate the accuracy of interpolated DEMs for a wider range of terrain types using many other interpolation models, for instance, the bi-cubic model, as a step towards the bilinear interpolation-based model. A possible approach is to first transform the nonlinear model into a linear one and then to apply to that linear model the techniques developed in this study related to linear models. Acknowledgments The research presented in this paper is support by Ministry of Science and Technology, China (Project No. 2012AA12A305).

Appendix A: Derivations of Eqs. (5) and (7) for the General Re-Sampling Ratio A.1 Formula Derivations for One-Dimensional Uniformly Distributed Points According to the elevation difference Eq. (4), together with 1st order Taylor expansion for the elevation function f at the sampling nodes 2 1  f (ri + j) − f¯ (ri + j) n i, j  2   j 1 j 1− = ( f (ri + j) − f (ri)) + ( f (ri + j) − f (ri + r )) n r r i, j

    2 j 2 2 j 1− ≤ ( f (ri + j) − f (ri))2 + ( f (ri + j) − f (ri + r ))2 n r r i, j

    j 2 f (ri + j) − f (ri) 2 2h 2  2 = j 1− n r jh i, j  2    j f (ri + j) − f (ri + r ) 2 2 + ( j − r) r ( j − r) h     2  2 2 2 j   j   2h 2  2 2 f (ri) + ( j − r ) f (ri + r ) ≈ j 1− n r r i, j     2   2

j 2   2h 2  2 f (ri) + f (ri + r ) = j 1− n r i, j

   n/r −1 r −1 2 2  2  j 2    2h 2 f (ri) + f  (ri + r ) j 1− = n r

(RMSE)2 =

j=1

123

i=0

Math Geosci (2014) 46:445–481

467

 2  2   n/r −1    4 f (ri) + f  (ri + r ) r − 1 h 2 i=0 = 15r 2 n/r  2    4 2 r 4 − 1 s2 2 r −1 h 2 E{(Slope) } = E{(Slope)2 }. = 15r 2 15r 4 Note, for the above derivation, when j = 0, the results should be zero for the elevation values of those interpolation nodes that are exactly the same as the original true values. Because forward difference quotients have been introduced, the sampling size s should be restricted within a very small range; otherwise, the average slope cannot be approximated in this fashion and could not represent the terrain map in actual applications.

A.2 Formula Derivations for Two-Dimensional Grid-Based DEM According to the linear interpolation method for the two-dimensional case, the following equation can be easily obtained as follows       l k l k ¯ 1− f (ri, r j) + 1− f (ri + r, r j) f (ri + k, r j + l) = 1 − r r r r       k l k l + 1− f (ri, r j + r )+ f (ri + r, r j + r ) . r r r r Therefore, f (ri + k, r j + l) − f¯ (ri + k, r j + l)    l k 1− = 1− ( f (ri + k, r j + l) − f (ri, r j)) r r    l k 1− + ( f (ri + k, r j + l) − f (ri + r, r j)) r r    l k + 1− ( f (ri + k, r j + l) − f (ri, r j + r )) r r    l k + ( f (ri + k, r j + l) − f (ri + r, r j + r )) . r r According to the Cauchy–Schwarz inequality x, y ≤ x y, if y = ( 1 1 · · · 1 )T , then the following formula will hold  n 2 n   xi ≤n xi2 . i=1

i=1

123

468

Math Geosci (2014) 46:445–481

The final square of the RMSE can be divided into two types of series summations as follows 1  ( f (ri + k, r j + l) − f¯(ri + k, r j + l))2 mn i, j,k,l ⎛ 1  ⎝ = ( f (ri + k, r j + l) − f¯(ri + k, r j + l))2 mn i, j kl=0   + ( f (ri + k, r j + l) − f¯(ri + k, r j + l))2 .

(RMSE)2 =

kl=0

and  kl=0

( f (ri +k, r j +l)− f¯(ri +k, r j +l))2 ⎛

⎜ ⎜ ⎜ ⎜ = ⎜ kl=0 ⎜ ⎝

  ⎞2 1− rk 1− rl ( f (ri +k, r j +l)− f (ri, r j)) ⎟    + rk 1− rl ( f (ri +k, r j +l)− f (ri +r, r j )) ⎟ ⎟ ⎟ l   k + 1− r r ( f (ri +k, r j +l)− f (ri, r j +r )) ⎟ ⎟ ⎠    + rk rl ( f (ri +k, r j +l)− f (ri +r, r j +r ))

⎛

⎞  2 (1− rk )2 1− rl ( f (ri +k, r j +l)− f (ri, r j))2 ⎜ ⎟  2  2 ⎜ ⎟ + rk 1− rl ( f (ri +k, r j +l)− f (ri +r, r j))2 ⎟ ⎜ ⎜ ⎟ ≤4 ⎜ ⎟   ⎜ +(1− k )2 l 2 ( f (ri +k, r j +l)− f (ri, r j +r ))2 ⎟ r r kl=0 ⎜ ⎟ ⎝ ⎠  k 2  l 2 2 + r ( f (ri +k, r j +l)− f (ri +r, r j +r )) r ⎛ ⎜ ⎜ ⎜ ⎜ 2 ≈ 4h ⎜ ⎜ kl=0 ⎜ ⎝ ⎛

⎞ 2  2 1− rk 1− rl ((k ∂∂x +l ∂∂y ) f (ri, r j))2 ⎟  2  2 ⎟ + rk 1− rl (((k −r ) ∂∂x +l ∂∂y ) f (ri +r, r j ))2 ⎟ ⎟ ⎟   k 2 l 2 ∂ ∂ 2 ⎟ +(1− r ) r ((k ∂ x +(l −r ) ∂ y ) f (ri, r j +r )) ⎟ ⎠  k 2  l 2 ∂ ∂ 2 + r (((k −r ) ∂ x +(l −r ) ∂ y ) f (ri +r, r j +r )) r

⎞  2 (1− rk )2 1− rl (k 2 ( f x (ri, r j))2 +l 2 ( f y (ri, r j))2 ) ⎜ ⎟ 2  2  ⎜ ⎟ 1− rl ((k −r )2 ( f x (ri +r, r j))2 +l 2 ( f y (ri +r, r j))2 ) ⎟ ⎜ + rk ⎜ ⎟ ⎜ ⎟   ⎜ +(1− k )2 l 2 (k 2 ( f x (ri, r j +r ))2 +(l −r )2 ( f y (ri, r j +r ))2 ) ⎟ ≤ 8h 2 r r ⎜ ⎟ ⎟ kl=0 ⎜ ⎜ +  k 2  l 2 ((k −r )2 ( f (ri +r, r j +r ))2 ⎟ ⎜ ⎟ x r r ⎝ ⎠ +(l −r )2 ( f y (ri +r, r j +r ))2 )

123

Math Geosci (2014) 46:445–481

⎛

469

(k 2 (r −k)2 (r −l)2 ( f x (ri, r j))2 +l 2 (r −k)2 (r −l)2 ( f y (ri, r j))2 )

⎞

⎟ ⎜ ⎜ +(k 2 (r −k)2 (r −l)2 ( f x (ri +r, r j))2 +k 2 l 2 (r −l)2 ( f y (ri +r, r j ))2 ) ⎟ ⎟ 8h 2  ⎜ ⎟ ⎜ = 4 ⎜ +(k 2 l 2 (r −k)2 ( f x (ri, r j +r ))2 +l 2 (r −k)2 (r −l)2 ( f y (ri, r j +r ))2 ) ⎟, r ⎟ ⎜ kl=0 ⎝ ⎠ +(k 2 l 2 (r −k)2 ( f x (ri +r, r j +r ))2 +k 2 l 2 (r −l)2 ( f y (ri +r, r j +r ))2 ) ⎛ 8h 2

=

r4

k 2 l 2 (r −k)2 (( f x (ri, r j +r ))2 +( f x (ri +r, r j +r ))2 )

⎜ ⎟ 2 2 2 2 2 ⎟ ⎜  ⎜ +k l (r −l) (( f y (ri +r, r j)) +( f y (ri +r, r j +r )) ) ⎟ ⎜ ⎟ ⎜ +k 2 (r −k)2 (r −l)2 (( f x (ri, r j))2 +( f x (ri +r, r j ))2 ) ⎟ ⎟ kl=0 ⎜ ⎝ ⎠ +l 2 (r −k)2 (r −l)2 (( f y (ri, r j))2 +( f y (ri, r j +r ))2 ) ⎛

kl=0 k

8h 2 = 4 r

⎞

2 l 2 (r −k)2 (( f

2 2 x (ri, r j +r )) +( f x (ri +r, r j +r )) )

⎞

⎜ ⎟ 2 2 2 2 2 ⎟ ⎜ + kl=0 k l (r −l) (( f y (ri +r, r j)) +( f y (ri +r, r j +r )) ) ⎟ ⎜ ⎜ ⎟ 2 2 2 2 2 ⎟. ⎜ + kl=0 k (r −k) (r −l) (( f x (ri, r j)) +( f x (ri +r, r j)) ) ⎟ ⎜ ⎝ ⎠  + kl=0 l 2 (r −k)2 (r −l)2 (( f y (ri, r j))2 +( f y (ri, r j +r ))2 )

Because these following equations hold up 

k 2 l 2 (r − k)2 =

kl=0

k=1



r −1 

k 2 l 2 (r − l)2 =

kl=0



r −1 

k (r − k) (r − l) = 2

2

2

kl=0

k=1



r −1 

l 2 (r − k)2 (r − l)2 =

kl=0

(r − 1)(2r − 1)r 2 (r 4 − 1) , 180

l 2 (r − l)2 =

(r − 1)(2r − 1)r 2 (r 4 − 1) , 180

l=1

k2

k=1 r −1 

r −1 

l2 =

k 2 (r − k)2 r −1  l=1

k (r − k) 2

2

r −1 

(r − l)2 =

(r − 1)(2r − 1)r 2 (r 4 − 1) , 180

l 2 (r − l)2 =

(r − 1)(2r − 1)r 2 (r 4 − 1) . 180

l=1

(r − k)2

k=1

r −1  l=1

therefore, we have 

( f (ri + k, r j + l) − f¯(ri + k, r j + l))2

kl=0

=

⎛ 2(r −1)(2r −1)(r 4 −1)h 2 45r 2

( f x (ri, r j))2 + ( f y (ri, r j))2

⎞

⎜ ⎟ ⎜ +( f x (ri, r j + r ))2 + ( f y (ri, r j + r ))2 ⎟ ⎜ ⎟ ⎜ ⎟. ⎜ +( f x (ri + r, r j))2 + ( f y (ri + r, r j))2 ⎟ ⎜ ⎟ ⎝ ⎠ +( f x (ri + r, r j + r ))2 + ( f y (ri + r, r j + r ))2

123

470

Math Geosci (2014) 46:445–481

Meanwhile  kl=0

=

2

f (ri + k, r j + l) − f¯(ri + k, r j + l) 

( f (ri + k, r j + l) − f¯(ri + k, r j + l))2

k(l=0)



+

( f (ri + k, r j + l) − f¯(ri + k, r j + l))2

l(k=0)

=



2   2 f (ri + k, r j) − f¯(ri + k, r j) + f (ri, r j + l) − f¯(ri, r j + l)

k



l

    2 k k f (ri + k, r j) − ( 1 − f (ri, r j) + f (ri + r, r j)) = r r k      2  l l + f (ri, r j + l) − 1− f (ri, r j) + f (ri, r j + r ) r r l    2   k k = 1− ( f (ri + k, r j)− f (ri, r j))+ ( f (ri +k, r j)− f (ri +r, r j)) r r k     2  l l 1− ( f (ri, r j +l)− f (ri, r j))+ ( f (ri, r j + l)− f (ri, r j +r )) + r r l    2   k k k 1− f x (ri, r j) + (k − r ) f x (ri + r, r j) ≈ h2 r r k      2  l l 2 l 1− f y (ri, r j) + (l − r ) f y (ri, r j + r ) +h r r l      2  k 2 k 2 2 2 2 2 = 2h ( f x (ri, r j)) + (k − r ) ( f x (ri + r, r j)) k 1− r r k      2  l 2 l 2 2 2 2 2 ( f y (ri, r j)) + (l − r ) ( f y (ri, r j + r )) +2h l 1− r r l     r −1

 k 2 2 2 = 2h ( f x (ri, r j))2 + ( f x (ri + r, r j))2 k 1− r k=1     r −1

 l 2 2 2 ( f y (ri, r j))2 + ( f y (ri, r j + r ))2 +2h l 1− r l=1

=

 − 1)h 2  ( f x (ri, r j))2 + ( f x (ri + r, r j))2   15r + ( f y (ri, r j))2 + ( f y (ri, r j + r ))2 (r 4

123

Math Geosci (2014) 46:445–481

471

Note that in this case, 0 ≤ i ≤ m/r, i ∈ N ; 0 ≤ j ≤ n/r, j ∈ N . Therefore, (RMSE)2 =

1  ( f (ri + k, r j + l) − f¯(ri + k, r j + l))2 mn i, j,k,l

1  = mn



kl=0 ( f (ri

+

i, j

⎛



+ k, r j + l) − f¯(ri + k, r j + l))2

kl=0 ( f (ri



+ k, r j + l) − f¯(ri + k, r j + l))2 ⎛

( f x (ri, r j))2 + ( f y (ri, r j))2

⎞⎞

⎜ ⎜ ⎟⎟ 2 2 ⎜ ⎜ ⎟⎟ ⎜ 2(r −1)(2r −1)(r 4 −1)h 2 ⎜ +( f x (ri, r j + r )) + ( f y (ri, r j + r )) ⎟⎟ ⎜ ⎜ ⎟⎟ ⎜ ⎜ +( f x (ri + r, r j))2 + ( f y (ri + r, r j))2 ⎟⎟ 45r 2 ⎜ ⎟⎟ 1 ⎜ ⎜ ⎝ ⎠⎟ ≤ ⎜ 2 + ( f (ri + r, r j + r ))2 ⎟ +( f (ri + r, r j + r )) ⎜ ⎟ x y mn i, j ⎜ ⎟   ⎜ ⎟ 4 2 2 + ( f (ri + r, r j))2 ⎜ + (r −1)h ⎟ ( f (ri, r j)) x x 15r ⎜ ⎟ ⎝ ⎠   + ( f y (ri, r j))2 + ( f y (ri, r j + r ))2 ⎛

⎛

⎞⎞

( f x (ri, r j))2 + ( f y (ri, r j))2

⎜ ⎜ ⎟⎟ ⎜ ⎜ +( f x (ri, r j + r ))2 + ( f y (ri, r j + r ))2 ⎟⎟ ⎜ ⎜ ⎟⎟ ⎜ ⎟⎟ ⎜2(r −1)(2r −1)r 2 (r 4 −1)h 2 m/r −1 n/r −1 ⎜ ⎟ 2 2 ⎜ +( f x (ri + r, r j)) + ( f y (ri + r, r j)) ⎟ ⎜ ⎜ ⎟⎟ 4 i=0 j=0 ⎜ ⎟ 45r ⎜ ⎟ ⎜ +( f x (ri + r, r j + r ))2 ⎟⎟ 1 ⎜ ⎜ ⎜ ⎟⎟ = ⎜ ⎝ ⎠⎟ ⎟ mn ⎜ ⎜ ⎟ +( f y (ri + r, r j + r ))2 ⎜ ⎟  ⎜ ⎟    4 2 m/r −1 n/r ⎜ + (r −1)h ⎟ 2 2 ⎜ ⎟ i=0 j=0 ( f x (ri, r j)) + ( f x (ri + r, r j)) 15r ⎜ ⎟   m/r n/r −1  ⎝ ⎠ + i=0 j=0 ( f y (ri, r j))2 + ( f y (ri, r j + r ))2 ≈

1 mn 

= =



! 2(r −1)(2r −1)(r 4 −1)h 2 4mn (r 4 − 1)h 2 2mn 2 2 ( ( E{(Slope) }) + E{(Slope) }) 45r 2 r2 15r r2

 2r (r 4 − 1)h 2 8(r − 1)(2r − 1)(r 4 − 1)h 2 2 2 E{(Slope) } + E{(Slope) } 45r 4 15r 4

2(8r 2 − 9r + 4)(r 4 − 1)h 2 E{(Slope)2 }. 45r 4

Appendix B: Numerical Results Involved in Sect. 4.4 See Tables 8 and 9, and Figs. 10 and 11.

123

472

Math Geosci (2014) 46:445–481

Table 8 The experimental results regarding different levels of terrain complexity with r = 3, k = 12.4001 s0 = 3h 0

Terrain parameters A

B

E{(Slope)2 }

(RMSE)2

3/16r 2 E{(Slope)2 }

C

3

10

1/3

300

0.0001

0.1007

1.2639

9

30

1

300

0.0005

0.9066

11.3755

30

100

3

300

0.0059

10.0688

126.3093

300

750

30

300

0.0136

23.3973

292.1342

50

150

9

300

0.0235

40.2752

505.2372

60

200

6

300

0.0529

90.6192

1136.7836 2020.9487

90

300

9

300

0.0940

161.1007

120

400

12

300

0.1474

252.4158

3169.7007

400

900

45

300

0.1587

275.6253

3414.4029

150

500

24

300

0.3566

619.3115

7670.6090

200

500

25

300

0.5362

932.3065

11534.1245

450

1,200

50

300

0.8939

1549.3610

19226.4759

500

1,500

50

300

1.3511

2330.3482

29060.3081

600

1,600

60

300

1.5880

2752.8284

34157.1181

700

1,800

80

300

2.0357

3532.3778

43786.1728

800

2,000

90

300

2.5379

4406.9707

54587.7962

1,000

2,000

100

300

2.8075

4871.9846

60386.4213

Table 9 The experimental results regarding different levels of terrain complexity with r = 4, k = 13.5462 s0 = 4h 0

Terrain parameters A

B

E{(Slope)2 }

(RMSE)2

3/16r 2 E{(Slope)2 }

C

3

10

1/3

300

0.0001

0.1042

1.4276

9

30

1

300

0.0005

0.9374

12.8486

30

100

3

300

0.0060

10.4101

142.6657

300

750

30

300

0.0138

24.1914

329.9654

50

150

9

300

0.0239

41.6405

570.6629

60

200

6

300

0.0537

93.6912

1283.9914

90

300

9

300

0.0955

166.5621

2282.6514

120

400

12

300

0.1498

260.9816

3580.1618

400

900

45

300

0.1613

284.9758

3856.5941

150

500

24

300

0.3624

640.3104

8664.0109

200

500

25

300

0.5450

963.9659

13027.9458

450

1,200

50

300

0.9084

1601.8853

21716.3895

500

1,500

50

300

1.3730

2409.3180

32823.5924

600

1,600

60

300

1.6138

2846.1302

38580.6148 49456.7532

700

1,800

80

300

2.0688

3652.1598

800

2,000

90

300

2.5791

4556.4362

61657.3305

1,000

2,000

100

300

2.8531

5037.7686

68207.6695

123

Math Geosci (2014) 46:445–481

473

Fig. 10 The fitted relationship between {RMSE}2 and its upper bound with r = 3, k = 12.4001

Fig. 11 The fitted relationship between {RMSE}2 and its upper bound with r = 4, k = 13.5462

Appendix C: Supplementary Explanation for Sect. 3 This self-contained supplementary material is provided for the special case with resampling ratio r = 2 for a better understanding of the detailed derivation process of Sect. 3.

123

474

Math Geosci (2014) 46:445–481

C.1 Accuracy Analysis for One-Dimensional, Uniformly Distributed Points with Re-Sampling Ratio r = 2 For a one-dimensional case, the original DEM, the corresponding re-sampled DEM and the linear interpolated DEM are denoted by f (i), f D (i) and f¯(i), respectively. It is assumed that the original DEM is composed of a point set {xi }, (i = 0, 1, 2, . . . , n), the resolution of the original DEM is h, the re-sampling ratio is r , and the resolution of the re-sampled DEM is s. A schematic diagram of the original DEM (Fig. 12a), the re-sampled DEM (Fig. 12b), and the interpolated DEM (Fig. 12c) are given in Fig. 12 where the re-sampling ratio is set as 2:1 (i.e., s = r h, r = 2) for the sake of simplicity to start the construction of our theoretical models. Using the linear interpolation method for the one-dimensional case, it is easy to obtain the following equations f (i) − f¯(i) = 0 (i = 0, 2, . . . , 2k, k ∈ N , 2k ≤ n), 1 f (i) − f¯(i) = f (i) − [ f (i + 1) + f (i − 1)] 2 1 = f (i) − f (i − 1) − [ f (i + 1) − f (i − 1)] 2 (i = 1, 3, . . . , 2k + 1, k ∈ N , 2k + 1 ≤ n).

(17)

(18)

For the last term of Eq. (18), we notice that f (i + 1) − f (i − 1) = f (i + 1) − f (i) + f (i) − f (i − 1).

(19)

Therefore,   f (i) − f (i − 1) 1 f (i + 1) − f (i) f (i) − f (i − 1) f (i) − f¯(i) = − + h h 2 h h 1 = f  (i − 1) − ( f  (i) + f  (i − 1)) 2 1  1 = f (i − 1) − f  (i). (20) 2 2 where f  (i − 1) and f  (i) denote the forward difference quotients which approximate the derivatives of function f (x) at the locations i − 1 and i, respectively. Because f (x) is the terrain function, f  (i) can be viewed as a depiction of the terrain slope at point i. Remembering that for the one-dimensional case, the error indicator, the root mean square error (RMSE), has the form " RMSE =

123

 i

( f (i) − f¯(i))2 . n

(21)

Math Geosci (2014) 46:445–481

475

f

x0 x1 x2 x3 x4 x5 x6 x7 x8 x9 x10

xn

x

xn

x

xn

x

(a) The original DEM f(i)

f

x0 x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 (b) The re-sampled DEM fD(i)

f

x0 x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 –

(c) The linear interpolated DEM f (i) Fig. 12 Original one-dimensional DEM (a), the re-sampled DEM (b), and the linear interpolated DEM (c) on the re-sampled DEM (b)

123

476

Math Geosci (2014) 46:445–481

Additionally, considering Eqs. (17) and (18), we have   ¯ 2 ( f (i)− f¯(i))2 ( f (i)− f¯(i))2 i=even ( f (i)− f (i)) = + i=odd (RMSE) = n n n  2 ¯ ( f (i) − f (i)) = i=odd . (22) n 

2

i

Combining Eqs. (21) and (22) yields (RMSE)2 =

h2  h2  ( f  (i − 1) − f  (i))2 ≤ {[ f  (i − 1)]2 + [ f  (i)]2 }. 4n 2n i=odd

i=odd

(23) Notice that the slope of point i − 1 is defined as | f  (i − 1)|, and the average value of the square of the DEM slope can therefore be defined as E{(Slope)2 } =

n 1  2 1  [ f (i)] = {[ f  (i − 1)]2 + [ f  (i)]2 }. n n i=1

(24)

i=odd

If the resolution of the re-sampled DEM is denoted by s, then obviously s = 2h for this case. Therefore, Eq. (23) is simplified as (RMSE)2 ≤

1 h 2 E{(Slope)2 } = s 2 E{(Slope)2 } 2 8

(25)

In terms of Eq. (25) for the one-dimensional case, if the terrain slope can be estimated beforehand, then the upper bound of the RMSE of the bilinear interpolation result could be obtained. From Eq. (25), the upper bound for the square of the RMSE is directly proportional to the square of the re-sampled DEM resolution s 2 and inversely proportional to the average value of the square of the DEM slope E{(Slope)2 } with a scale factor of 1/8. This conclusion is consistent with our intuitive observation of DEM accuracy in relation to the slope and resolution factors. Thus, we have here provided a mathematical model (the inequality of Eq. 25) to quantify the relationship between DEM accuracy and (a) DEM resolution and (b) terrain slope. C.2 Accuracy Analysis for the Two-Dimensional Grid-Based DEM with Re-Sampling Ratio r = 2 For the two-dimensional case (i.e., two-dimensional DEM), it is assumed that the original DEM, the re-sampled DEM, and the interpolated DEM are denoted by f (i, j), f D (i, j) and f¯(i, j) (i, j = 1, 2, . . . , n), respectively, as shown in Fig. 13. The resolution of the original DEM is also assumed to be h, the re-sampling ratio is denoted by r , and the resolution of the re-sampled DEM is denoted by s. For the sake of simplicity, the re-sampling ratio along each of the xand y directions is set to be 2:1,

123

Math Geosci (2014) 46:445–481

(a) The original DEM

477

(b) The re-sampled DEM

(c) The bilinear interpolated DEM

Fig. 13 An illustration of two-dimensional data sets: the original two-dimensional DEM, the re-sampled DEM, and the bilinear interpolated DEM

i.e., s = r h, (r = 2) to begin the two-dimensional case discussion of the construction of the theoretical models. For the bilinear interpolation method in the two-dimensional case, ⎧ f (i, j − 1) + 21 [ f (i, j + 1) − f (i, j − 1)], i = even, j = odd ⎪ ⎪ ⎪ ⎪ 1 ⎪ ⎪ ⎨ f (i − 1, j) + 2 [ f (i + 1, j) − f (i − 1, j)], i = odd, j = even (26) f¯(i, j) = 41 [ f (i − 1, j − 1) + f (i + 1, j − 1) + f (i − 1, j + 1) ⎪ ⎪ ⎪ ⎪ + f (i + 1, j + 1)], i = odd, j = odd ⎪ ⎪ ⎩ f (i, j), i = even, j = even Now we calculate the value | f (i, j) − f¯(i, j)|2 in the following four cases: Case (1) i = even, j = odd 1 f (i, j) − f¯(i, j) = { f (i, j) − f (i, j − 1) − [ f (i, j + 1) − f (i, j − 1)]}/ h h 2 1 = { f (i, j) − f (i, j − 1) − [ f (i, j + 1) − f (i, j) 2 + f (i, j) − f (i, j − 1)]}/ h  1 = f y (i, j − 1) − f y (i, j) + f y (i, j − 1) 2 1 = [ f y (i, j − 1) − f y (i, j)]. 2 By the inequality (a − b)2 ≤ 2(a 2 + b2 ) (a, b ∈ R), 1 | f (i, j) − f¯(i, j)|2 = h 2 [ f y (i, j − 1) − f y (i, j)]2 4 1 ≤ h 2 [| f y (i, j − 1)|2 + | f y (i, j)|2 ], 2

(27)

where f y (·, ·) is the forward difference quotient with respect to the first variable y, which can be regarded as an approximation for the partial derivative of the function f , that is, the terrain slope along the direction of coordinate y.

123

478

Math Geosci (2014) 46:445–481

i-1, j+1

i, j+1

i+1, j+1

i-1, j

i, j

i+1, j

i-1, j-1

i, j-1

i+1, j-1

Fig. 14 Explanation for the difference scheme

Case (2) i = odd, j = even Following a similar approach to that of case (1), 1 | f (i, j) − f¯(i, j)|2 ≤ h 2 [| f x (i − 1, j)|2 + | f x (i, j)|2 ], 2

(28)

where f x (·, ·) is the forward difference quotient with respect to the first variable x, which could be viewed as an approximation for the partial derivative of the function f , that is, the terrain slope along the direction of coordinate x. Case (3) i = odd, j = odd As denoted in Fig. 14 above, 1 f (i, j) − f¯(i, j) = f (i, j) − [ f (i − 1, j + 1) + f (i − 1, j − 1) 4 + f (i + 1, j − 1) + f (i + 1, j + 1)] 1 = f (i, j) − {[ f (i − 1, j + 1) − f (i − 1, j)] + [ f (i − 1, j − 1)− f (i, j −1)] 4 + [ f (i + 1, j − 1) − f (i + 1, j)] + [ f (i + 1, j + 1) − f (i, j + 1)] + f (i − 1, j) + f (i, j − 1) + f (i + 1, j) + f (i, j + 1)} 1 = − {[ f (i − 1, j + 1) − f (i − 1, j)] + [ f (i − 1, j − 1) − f (i, j − 1)] 4 + [ f (i + 1, j − 1) − f (i + 1, j)] + [ f (i + 1, j + 1) − f (i, j + 1)] + [ f (i, j) − f (i − 1, j)] + [ f (i, j) − f (i, j − 1)] + [ f (i, j) − f (i + 1, j)] + [ f (i, j) − f (i, j + 1)]}

123

Math Geosci (2014) 46:445–481

479

Therefore, 1 f (i, j)− f¯(i, j) =− f y (i −1, j)− f x (i −1, j −1)− f y (i +1, j −1)+ f x (i, j +1) h 4  + f x (i − 1, j) + f y (i, j − 1) − f x (i, j) − f y (i, j) 1 = [ f x (i, j) − f x (i, j + 1) − f x (i − 1, j) + f x (i − 1, j − 1) 4  + f y (i, j) − f y (i, j − 1) − f y (i − 1, j) + f y (i + 1, j − 1) Thus, | f¯(i, j) − f (i, j)|2 1  ≤ h 2 [ f x (i, j) − f x (i, j + 1) − f x (i − 1, j) + f x (i − 1, j − 1)]2 8  2  + f y (i, j) − f y (i, j − 1) − f y (i − 1, j) + f y (i + 1, j − 1) 1  ≤ h 2 [ f x (i, j)]2 + [ f x (i, j + 1)]2 + [ f x (i − 1, j)]2 + [ f x (i − 1, j − 1)]2 2  2  2  2  2  · + f y (i, j) + f y (i, j − 1) + f y (i − 1, j) + f y (i + 1, j − 1) (29) where the meanings of f x (·, ·) and f y (·, ·) are the same as those of the previous two cases as explained above. Case (4) i = even, j = even By Eq. (26), obviously, | f (i, j) − f¯(i, j)|2 = 0.

(30)

For the two-dimensional case, the root mean square error (RMSE) is given as " RMSE =

i, j

( f (i, j) − f¯(i, j))2 n×n

.

Therefore, 1  | f (i, j) − f¯(i, j)|2 n2 i, j ⎛  1 = 2⎝ | f (i, j) − f¯(i, j)|2 + n

(RMSE)2 =

i=even, j=odd

+

 i=odd, j=odd

| f (i, j) − f¯(i, j)|2 +



| f (i, j)− f¯(i, j)|2

i=odd, j=even



⎞

| f (i, j) − f¯(i, j)|2⎠ .

i=even, j=even

123

480

Math Geosci (2014) 46:445–481

By using Eqs. (27)–(30) with s = 2h (s is resolution of the re-sampled DEM), (RMSE)2 ≤

3 2 3 2 h E{(Slope)2 } = s E{(Slope)2 } 4 16

(31)

References Almansa A, Cao F, Gousseau Y, Rougé B (2002) Interpolation of digital elevation models using AMLE and related methods. IEEE Trans Geosci Remote 40:314–325 Armstrong RN, Martz LW (2003) Topographic parameterization in continental hydrology: a study in scale. Hydrol Process 5:3763–3781 Büyüksalih G, Koçak G, Oruç M, Akcin H, Jacobsen K (2004) Accuracy analysis, DEM generation and validation using Russian TK-350 stereo-images. Photogramm Rec 19:200–218 Caselles V, Morel JM, Sbert C (1998) An axiomatic approach to image interpolation. IEEE Trans Image Process 7:376–386 Chang K, Tsai B (1991) The effect of DEM resolution on slope and aspect mapping. Cartogr Geogr Inf Sci 18:69–77 Chen C, Li Y (2013) A Robust multiquadric method for digital elevation model construction. Math Geosci 1–23 Florinsky IV (2002) Errors of signal processing in digital terrain modeling. Int J Geogr Inf Sci 16:475–501 Gao J (1997) Resolution and accuracy of terrain representation by grid DEMs at a micro-scale. Int J Geogr Inf Sci 11:199–212 Goodman T, Said H, Chang L (1995) Local derivative estimation for scattered data interpolation. Appl Math Comput 68:41–50 Gousie MB, Franklin WR (2005) Augmenting grid-based contours to improve thin plate dem generation. Photogramm Eng Rem S 71:69–79 Holmes K, Chadwick O, Kyriakidis PC (2000) Error in a USGS 30-meter digital elevation model and its impact on terrain modeling. J Hydrol 233:154–173 Hutchinson MF (1988) Calculation of hydrologically sound digital elevation models. In: Proceedings of the third international symposium on spatial data handling, Sydney Hutchinson MF, Gallant JC (1999) Representation of terrain. In: Longley PA, Maguire DJ, Goodchild MF (eds) Geographical information systems—principles and technical issues, 2nd edn. Wiley, New York, pp 105–124 Hutchinson MF, Gallant JC (2000) Digital elevation models and representation of terrain shape. Terrain analysis: principles and applications, pp 29–49 Hutchinson MF (2011) ANUDEM Version 5.3. Australia’s Knowledge Gateway Web. http://fennerschool. anu.edu.au/files/usedem53_pdf_16552.pdf. Accessed 20 August 2013 Jia G, Wang X, Wei H (2013) An effective approach for selection of terrain modeling methods. IEEE Geosci Remote S 10:875–879 Kidner DB (2003) Higher-order interpolation of regular grid digital elevation models. Int J Remote Sens 24:2981–2987 Kyriakidis PC, Goodchild MF (2006) On the prediction error variance of three common spatial interpolation schemes. Int J Geogr Inf Sci 20:823–855 Li ZL, Zhu Q (2003) Digital elevation models, 2nd edn. Wuhan, China Maune D (2007) Digital elevation model technologies and applications: the DEM user’s manual, 2nd edn. Bethesda, USA Moore ID, Grayson RB, Ladson AR (1991) Digital terrain modeling: a review of hydrological, geomorphological and biological applications. Hydrol Process 5:3–30 National Research Council Committee on Floodplain Mapping Technologies (2007) Elevation data for floodplain mapping. National Academy Press, Washington Pumar M (1996) Zooming of terrain imagery using fractal-based interpolation. Comput Graph 20:171–176 Raaflaub LD, Collins MJ (2006) The effect of error in gridded digital elevation models on the estimation of topographic parameters. Environ Modell Softw 21:710–732 Rees W (2000) The accuracy of digital elevation models interpolated to higher resolutions. Int J Remote Sens 21:7–20

123

Math Geosci (2014) 46:445–481

481

Shan J, Zaheer M, Hussain E (2003) Study on accuracy of 1-degree DEM versus topographic complexity using GIS zonal analysis. J Surv Eng-Asce 29:85–89 Shi WZ, Li Q, Zhu CQ (2005) Estimating the propagation error of DEM from higher-order interpolation algorithms. Int J Remote Sens 26:3069–3084 Shi WZ, Tian Y (2006) A hybrid interpolation method for the refinement of a regular grid digital elevation model. Int J Geogr Inf Sci 20:53–67 Sinha SS, Schunck BG (1992) A two-stage algorithm for discontinuity-preserving surface reconstruction. IEEE Trans Pattern Anal Mach Intell 14:36–55 Taconet O, Vannier E, Le Hégarat-Mascle S (2010) A contour-based approach for clods identification and characterization on a soil surface. Soil Till Res 109:123–132 Tang GA, Gong JY, Chen ZJ, Cheng YH, Wang ZH (2001) A simulation on the accuracy of DEM terrain representation. Acta Geod 30:361–365 Toutin T (2002) Impact of terrain slope and aspect on radar grammetric DEM accuracy. ISPRS J Photogramm 57:228–240 Vaze J, Teng J, Spencer G (2010) Impact of DEM accuracy and resolution on topographic indices. Environ Modell Softw 25:1086–1098 Wang GX, Zhu CQ, Shi WZ, Zhang GQ (2004) The further study on the error accuracy of DEM terrain representation. Acta Geod 33:168–173 Wu H, Sun Y, Shi WZ, Chen X, Fu D (2013) Examining the satellite-detected urban land use spatial patterns using multidimensional fractal dimension indices. Remote Sens 5:5152–5172 Wu H, Zhou L, Chi X, Li Y, Sun Y (2012) Quantifying and analyzing neighborhood configuration characteristics to cellular automata for land use simulation considering data source error. Earth Sci Inform 5:77–86 Yamazaki D, Baugh CA, Bates PD, Kanae S, Alsdorf DE, Oki T (2012) Adjustment of a spaceborne DEM for use in floodplain hydrodynamic modeling. J Hydrol 436:81–91 Yang Q, McVicar TR, Van Niel TG, Hutchinson MF, Li L, Zhang X (2007) Improving a digital elevation model by reducing source data errors and optimising interpolation algorithm parameters: An example in the Loess Plateau, China. Int J Appl Earth Obs 9:235–246 Zhou QM, Liu XJ (2002) Error assessment of grid-based flow routing algorithms used in hydrological models. Int J Geogr Inf Sci 16:819–842 Zhou QM, Liu XJ (2003) The accuracy assessment on algorithms that derive slope and aspect from DEM In: Shi WZ, Goodchild MF, Fisher PF (Eds) Proceedings of the second international symposium on spatial data quality, 19–20 March, Hong Kong, pp 275–285 Zhou QM, Liu XJ (2004) Analysis of errors of derived slope and aspect related to DEM data properties. Comput Geosci 30:369–378 Zhu C, Shi W, Li Q et al (2005) Estimation of average DEM accuracy under linear interpolation considering random error at the nodes of TIN model. Int J Remote Sens 26:5509-5523

123