Estimating slope from raster data--a test of eight algorithms at different resolutions in flat and steep terrain.

Tang, Jing ; Pilesjo, Petter ; Persson, Andreas 等

Introduction

A Digital Elevation Model (DEM) is a digital representation of a portion of Earth's surface, and is widely utilized to extract topographic variables, to be used in hydrological, geomorphological and biological applications (Wilson, Gallant 2000). The accuracies of the extracted variables, such as slope, aspect, specific catchment area, and flow length, are directly related to DEM quality as well as algorithm selections, and those factors have been studied by many researchers (Woods et al. 1997; Guentner et al. 2004; Kopecky, Cizkova 2010). Slope, as one of the most commonly used variables in many models, has been estimated from DEMs by applying a large number of different algorithms (Zhou, Liu 2004; Jones 1998; Florinsky 1997; Warren et al. 2004). Different Geographical Information System (GIS) software (e.g. ArcInfo, ERDAS Imagine, GRASS, and PCRaster) also use different slope algorithms. The choice of algorithm can result in different slope estimations, even when the estimation is based on the same DEM data (Warren et al. 2004; Skidmore 1989), making further modelling and interpretation using slope values uncertain and possibly faulty. Renard et al. (1991) reported that a deviation of 10% in slope estimation can result in a 20% error in soil erosion estimation. An increase in (estimated) slope and a decrease in specific drainage area estimation can result in a decrease in estimated topographic wetness index (TWI), and thus an underestimation of soil wetness (Hasan et al. 2012). A study of the correlation between TWI and vascular plant species richness, soil pH, ground level, and soil moisture was conducted in northern Sweden (Sorensen et al. 2006), concluding that diverse slope estimations, resulting in varying estimated TWI values, do influence the correlation with the aforementioned variables. The choice of suitable slope algorithms should thus be done with caution, especially when using high accuracy data (Irfan Ashraf et al. 2012).

The most frequently used numerical slope estimation algorithms are calculating the elevation differences of the DEM within a moving 3*3 windows (Jones 1998). In this study, eight frequently used algorithms are selected for comparison (see slope algorithms below). The selected algorithms correspond to those implemented in commonly available GIS software (see above). Algorithms not based on 3*3 moving windows, like the downslope index (Hjerdt et al. 2004), depending on downslope elevation differences and flow path, are not included in this paper. It is widely accepted that DEM quality is an important source of errors in slope estimation, but, additionally, the errors produced by different algorithms cannot be overlooked. Some authors have evaluated the accuracy of various slope algorithms, using different assessment methodologies (Zhou, Liu 2004). In those studies, the methodology of ranking the 'best' algorithms is normally based on artificial surfaces, like the Morrison surface (Jones 1998), the Gauss synthetic surface, and the ellipsoid (Zhou, Liu 2004). However, it might be more effective and meaningful to rank algorithms based on their performance with respect to the real world surfaces that resemble those to which they will be operationally applied. It is also desirable to discuss intrinsic features of the individual algorithms. Additionally, the estimations of slope are always dependent on the resolution and quality of the DEM (Thompson et al. 2001; Chang, Tsai 1991; Kienzle 2004), and this must be taken into account.

There is a lack of studies that aim to investigate possible differences between different slope algorithms, as well as algorithmic sensitivity to different types of terrain in connection with different degrees of generalisation (resolution) differences in DEMs. The previous findings of 'better' or 'worse' slope algorithms do not provide sufficient information when it comes to selection of algorithm for a particular application, and it is not unusual that researchers are unaware of, or neglect, the importance of selecting an appropriate slope algorithm. This might significantly influence the results.

Given the importance of slope estimations in many applications and the wide availability of Airborne Light Detection And Ranging (LiDAR) data, e.g. the dataset used in this study (see Fig.1), there is a strong demand to study possible differences between different slope algorithms and investigate how they behave in different sorts of terrain, represented with different DEM resolutions. In this study, eight frequently used slope algorithms are evaluated against each other. They are applied in two distinct terrain forms, one flat area and one steep area in the Stordalen catchment, located in northern Sweden. Four levels of resolution, 0.5 m, 1 m, 2 m and 5 m of the DEM are created in order to investigate the effects of DEM resolution on slope estimation. The aim of this study is to examine the internal differences and characteristics of the selected algorithms, and their sensitivities to the different DEM resolutions and terrain forms. The results will highlight differences and thus help researchers to choose appropriate slope algorithms when modelling terrain.

Specifically, the following research questions are addressed:

1. How does the terrain form (flat and steep areas) influence the relative differences of estimated slope between different algorithms?

2. With changing DEM resolution, are these algorithms sensitive to the variation in information content of the DEM and if so, how do they change? How do the relative differences between the algorithms change at different resolutions?

3. Since no grading of 'best' or 'worst' algorithm is performed in this study, the inherent characteristics of eight algorithms will be studied. How do these algorithms estimate the slope values, referring to overestimation or underestimation in relation to others?

1. Study area and data

In this study LiDAR data are used to generate DEMs of different resolution. LiDAR is a laser-based active remote sensing technique, and considered to be one of the most effective and reliable methods of terrain data collection (Liu 2008), resulting in a dense sample of highly accurate elevation points. Generally LiDAR derived DEMs are more accurate than products based on traditional methods like topographic maps and photographic interpretations, and thus are widely used in environmental modelling. In general, a high quality DEM results in better estimates of slope values. In this study, different slope algorithms are applied on a LiDAR generated DEM from northern Sweden. This DEM has been used in many studies, including a test of incorporation of topographic indices into the dynamic ecosystem model, LPJ-GUESS (Tang et al. manuscript).

Our study area is located in the subarctic Stordalen peatland, northern Sweden. The LiDAR data points cover the whole area and the total number of the measured elevation points is 76 940 341, with an average point density of approx. 13 points/[m.sup.2]. The retrieved data were post processed using 54 known points from the national geodetic network, resulting in an average vertical magnitude of errors of 0.022 m and an RMSE of 0.031 m (Hasan et al. 2012). The Stordalen peat-land has been included in many research projects and has been thoroughly studied since the early 20th century (Ryden et al. 1980, Olefeldt, Roulet 2012, Christensen et al. 2004). The acquisition of high resolution LIDAR data has mainly been used for analyses of hydrological processes, as well as examining the relationships between extracted topographic variables and other environmental processes. For our study, two different areas covering 100*100 m (see Fig. 1), are selected from the Stordalen peatland. The flat area (see Fig. 1, left) is characterized by smaller changes in elevation (max 0.7 m), while the terrain in the steeper area (see Fig. 1, right) is characterized by larger differences in elevation (max 42 m).

[FIGURE 1 OMITTED]

2. Methodology

2.1. Slope algorithms

At every point in a DEM the slope (S) can be defined as a function of gradients in the X and Y directions:

S = arctan [square root of ([f.sub.x.sup.2] + [f.sub.y.sup.2])]. (1)

The key in slope estimation is the computation of the perpendicular gradients [f.sub.x] and [f.sub.y]. Different algorithms use different techniques to estimate [f.sub.x] and f , resulting in a diversity in estimated slope values. As mentioned above, from a gridded DEM the common approach when estimating [f.sub.x] and f is to apply a moving 3*3 window (see Fig. 2) to derive the finite differential or local polynomial surface fit for the estimations (Florinsky 1997; Zhou, Liu 2004).

Below, eight frequently used slope algorithms tested in this study are briefly presented. Algorithms 1-5 and 8 are convolutional methods, based on approximation of differential operators by finite differences. Algorithm 6 compares the central elevation with its eight neighbours, adopting the largest difference in elevation. Algorithm 7 uses a quadratic regression surface constrained to go through the central elevation point of the local 3*3 window sampling kernel (Jones 1998). Before the descriptions of the different methods the numbering of the cells in the 3*3 window is defined (see Fig. 2), and the cell size (spatial resolution) is denoted as g.

In the following mathematical equations of slope, [z.sub.i] (i = 1, 2, ... 9) is the elevation value in cell i (defined in Fig. 2 above). The selected slope algorithms are the:

1. Second-order Finite Difference (2FD), (Fleming, Hoffer 1979)

[f.sub.x] = [([z.sub.6] - [z.sub.4])/2g], (2)

[f.sub.y] = [([z.sub.8] - [z.sub.2])/2g]. (3)

2. Third-order Finite Difference Weighted by Reciprocal of Distance (3FDWRD), (Unwin 1981)

[f.sub.x] = [([z.sub.3] - [z.sub.1] + [square root of (2)][[z.sub.6] - [z.sub.4]) + [z.sub.9] - [z.sub.7])/(4 + 2[square root of (2)])g], (4)

[f.sub.y] = [([z.sub.7] - [z.sub.1] + [square root of (2)]([z.sub.8] - [z.sub.2]) + [z.dub.9] - [z.sub.3])/(4 + 2[square root of (2)])g]. (5)

3. Third-order Finite Difference, linear regression plane (3FD), (Sharpnack, Akin 1969)

[f.sub.x] = ([z.sub.3] - [z.sub.1] + [z.sub.6] - [z.sub.4] + [z.sub.9] - [z.sub.7])/6g, (6)

[f.sub.y] = ([z.sub.7] - [z.sub.1] + [z.sub.8] - [z.sub.2] + [z.sub.9] - [z.sub.3])/6g. (7)

4. Third-order Finite Difference Weighted by Reciprocal of Squared Distance (3FDWRSD), (Horn 1981)

[f.sub.x] = [([z.sub.3] - [z.sub.1] + 2([z.sub.6] - [z.sub.4]) + [z.sub.9] - [z.sub.7])/8g], (8)

[f.sub.y] = [([z.sub.7] - [z.sub.1] + 2([z.sub.8] - [z.sub.2]) + [z.sub.9] - [z.sub.3])/8g]. (9)

[FIGURE 2 OMITTED]

5. Frame Finite Difference (FFD), (Chu, Tsai 1995)

[f.sub.x] = [([z.sub.3] - [z.sub.1] + [z.sub.9] - [z.sub.7])/4g], (10)

[f.sub.y] = [([z.sub.7] - [z.sub.1] + [z.sub.9] - [z.sub.3])/4g]. (11)

6. Maximum (Max), (O'Callaghan, Mark 1984; Travis et al. 1975)

max(abs(([z.sub.5] - [z.sub.1])/[square root of (2g)]), abs(([z.sub.5] - [z.sub.2])/g),

abs(([z.sub.5] - [z.sub.3])/[square root of (2g)]), abs(([z.sub.5] - [z.sub.9])/[square root of (2g)]),

abs(([z.sub.5] - [z.sub.7])/[square root of (2g)]), abs(([z.sub.5] - [z.sub.6])/g),

abs(([z.sub.5] - [z.sub.8])/g),abs(([z.sub.5] - [z.sub.4])/g)). (12)

7. Constrained Quadratic Surface (Quadsurface), (Wood 1996)

F(x,y) = [ax.sup.2] + [by.sup.2] + cxy + dx + ey + j, (13)

F(x,y) = Z = AX, (14)

where A has been defined according to Eq. 15 below,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (15)

X stands for the unknown vector of parameters, X = [(a b c d e f).sup.T] and Z is the elevation vector, Z = [([z.sub.9] [z.sub.8] [z.sub.7] [z.sub.6] [z.sub.5] [z.sub.4] [z.sub.3] [z.sub.2] [z.sub.1]).sup.T], so X can be calculated as:

X = [([A.sup.T]A).sup.-1] [A.sup.T]Z. (16)

Since the number of equations is more than the number of unknown parameters, the least-squares method is thus used to minimize the squares of the errors from each single equation and to determine the indices of the constrained quadratic surface. And it is relatively easy to estimate the [f.sub.x] and [f.sub.y] values at the centre of the 3*3 window (see Eq.17 and 18).

[f.sub.x] [|.sub.x=0,y=0] = d (17)

[f.sub.y] [|.sub.x=0,y=0] = e. (18)

8. Simple Difference (SimpleD), (Jones 1998)

[f.sub.x] = ([z.sub.5] - [z.sub.4])/g, (19)

[f.sub.y] = ([z.sub.5] - [z.sub.2])/g. (20)

2.2. Comparisons between estimated slope values

Since the aim of this study is to analyse the relative differences between the algorithms, the following preconditions are identified in order to strengthen the outcome of the study:

1. The elevation errors of the different DEMs are equivalent, since they are based on the same data set and interpolated in the same way;

2. All slope values are estimated based on a 3*3 moving window;

3. The comparisons between the estimated values are processed in the same way, independent of algorithm, terrain form and resolution;

4. The comparisons are aimed at highlighting the algorithms' behaviour applied on various terrains and resolutions. No reference values are used to rank the accuracy of the different algorithms.

Based on the original LiDAR data points the linear inverse distance interpolation method was applied on the two selected areas to generate DEMs with grid cell sizes of 0.5, 1, 2 and 5 m (Hasan et al. 2012). To distinguish the relative differences between the algorithms relating to terrain form, the terrain roughness of different DEMs can be quantified by implementing a scaling diversity index (SDI). SDI is a statistical index originally designed to estimate richness and evenness of ecological diversity, and it has been widely used in ecosystem studies (Yue et al. 2007). In this study it is used to estimate elevation diversity differences between the two terrain types at various resolutions.

SDI([epsilon],r) = -ln ([[summation].sup.i=m.sub.i=1] [[P.sub.i.sup.1/2]).sup.2]/ln[epsilon], (21)

[epsilon] = [(e + A).sup.-1], (22)

where [p.sub.i] is the proportion of each elevation value in relation to the whole investigation area; m is the total number of different elevation values; A is the area in hectares, e is the constant value of 2.71828, and r is the DEM resolution (Yue et al. 2007).

After quantifying the terrain roughness of each DEM at the different resolutions, slope values were estimated using the eight selected algorithms. The statistical analyses were conducted on these estimated slope values. The initial stage of the analysis compares the mean and standard deviation (STD) of the slope values for the eight different cases (four resolutions, two terrain types). This describes the general characteristics and relative differences between the different algorithms. After this, a one-way ANOVA analysis is carried out. The eight algorithms at one level of resolution are treated as a group when running the ANOVA analysis. The repetitive ANOVA testing goes through all the resolutions of DEMs, and the ANOVA F statistic is presented in Table 2. Furthermore, an assessment of the different slope algorithms by pairwise comparison is carried out. This is, from a statistical point of view, a comparison of the intrinsic differences of slope algorithms' results, instead of the more common approaches that use 'reference' or average values as 'ground-truth' data. A pairwise comparison reveals which pairs of algorithms are different (and which ones are not) at different resolutions and in different terrain types. A Matlab program was written to perform the multiple pairwise comparisons. A 95% confidence interval around the means is used as the standard for significant separation of the algorithms. The figures showing mean values and an interval around the mean values are provided for flat and steep terrains respectively. Multiple comparisons between different resolutions are carried out to depict the different algorithms' sensitivity to the changing information content of the DEMs. By combining the results of the analyses presented above, both the terrain form and resolution influences on the estimates produced by the different algorithms are documented.

3. Results

3.1. Scaling Diversity Index (SDI) for different DEM resolutions

Statistical analysis of elevation diversity (estimated by applying the SDI) at different spatial resolutions strongly indicates a reduced terrain roughness when reducing DEM resolution. The rates of reduction for flat and steep areas are not equal when changing resolution from 0.5 m to 5 m. The SDI decreases by 28.26% for the steep area, but only by 4.4% for flat area (see Fig. 3). Knowing there is a dramatic difference in the change of elevation diversity between flat and steep areas when changing the resolution, one can also expect differences in estimated slope. The changes of terrain roughness at the four resolutions are combined with the estimated slope values from different algorithms to reveal the effects of DEM generalisation on slope estimation.

3.2. Descriptive statistics

The statistical values of mean and standard deviation of estimated slope for flat and steep terrain at four levels of resolution are presented in Table 1. The mean values indicate possible overestimation or underestimation of each algorithm relative to the others over the same area. Generally a reduced DEM resolution results in lower estimated slope values for all algorithms, both in flat and steep terrain. All algorithms show a remarkable reduction in estimated mean slope in the flat area, where especially values estimated by the use of e.g. Max and SimpleD algorithms decrease by 84% (6.7 to 1.039) and 83.6% (5.056 to 0.826) respectively. Comparatively, the reduced percentage for Max and SimpleD are 19.74% (30.327 to 24.341) and 8.47% (26.134 to 23.92) in the steep area. The Max algorithm always produces the highest mean slope values, the SimpleD algorithm estimates the 2nd highest slope values, and the Quadsurface algorithm produces the lowest values, independent of terrain form and resolution. The differences between the Max and Quadsurface algorithms' mean values range from 0.669 (5 m resolution) to 3.931 (0.5 m resolution) in the flat area, and from 2.360 (5 m resolution) to 7.991 (0.5 m resolution) in the steep area.

[FIGURE 3 OMITTED]

Regarding comparisons of standard deviations (STD), reduced DEM resolution also decreases the variation in estimated slope values (see Table 1). In general, the SimpleD algorithm produces the highest values compared to the other algorithms at the same resolution, indicating a wide spread of estimated slope values. The Max algorithm produces the second highest variation in estimated slope. The Quadsurface method produces the lowest mean values and also produces the lowest variation in estimated slope. Besides, the SimpleD algorithm undergoes the greatest reductions of estimated variations from finer to coarser resolution. Those trends are the same both for the flat and the steep areas.

The comparisons of mean and standard variation values extract broad, non-specific information about the estimated slope values, indicating large differences between the tested algorithms. In the following ANOVA and pairwise comparisons, a more detailed analysis of the differences between the eight algorithms for different combinations of terrain form and DEM resolution is carried out.

3.3. Analysis of variance

The outputs of the ANOVA in Table 2 demonstrate that there is at least one of the slope algorithm estimations that is significantly different from the others (see P-value columns in Table 2) at the 0.01 level of significance, except for the case of the steep area at 5 m resolution. The mean square values of Between Groups ([MS.sub.b]) and Within Groups ([MS.sub.w]) describe the variability among algorithms and the variability within the estimated slope datasets, respectively. Comparing the same resolution outputs in flat and steep areas, the greater variation within each algorithm estimation (see column [MS.sub.w]) causes greater variances among algorithms (see column [MS.sub.b]). The F ([MS.sub.b]/[MS.sub.w]) value in the ANOVA test depicts the strength of those two sources of variances and more distinctive variations between groups comparing to the variations within each algorithm estimations are found in the flat area for the same resolution. Furthermore, for the steep area, the decreasing rate of the variability among algorithms (99.93% reduction of [MS.sub.b] value from 0.5 m to 5 m resolution) is more pronounced than the decreasing rate of the variability in the estimated slope (51.9% reduction of [MS.sub.w] value), which shows that no statistical significance is found among algorithms at 99% probability (P-value = 0.039). When compared to the steep area, the lesser topographical variation present in the flat area makes the smoothing effects from the resolution generalization less effective, so the differences between the algorithms become more prominent.

But regardless of whether we consider the flat or steep area, the increase of DEM grid cell size does influence the differences among algorithms. Even so, we still need more information about how particular algorithms differentiate themselves from others. The next pairwise comparison will analyse the detailed algorithmic differences.

3.4. Pairwise comparisons

The pairwise comparison depicts significant differences between the eight algorithms at different resolutions. Algorithms falling in the same group show no significant difference from each other, while algorithms falling in separate groups do differ.

In Fig. 4 and 5 the outputs of the multiple pair-wise comparisons are presented. The slope values are plotted along the x-axis, and the different algorithms are separated along the y-axis. Each algorithm is represented by a circle, and a horizontal line. The middle of the circle corresponds to the mean value and the line shows the interval (min-max value) around the mean. In order to highlight differences between different resolutions (0.5, 1, 2 and 5 m) the scale on the x-axis is maintained. Partly due to this choice of presentation, the intervals around the circle are relatively small for the flat area. The grey colour between the two dashed vertical lines represents the same group as the algorithm in blue, and the rest of algorithms in red represent one group each.

[FIGURE 4 OMITTED]

For the flat area (see Fig. 4), the change in resolution from finer to coarser does not change the groupings between the algorithms, except for the FFD method at the 2 m resolution, where it forms a group by itself. To summarize, the results show that the 2FD, Max, and SimpleD algorithms are significantly different from the other, producing higher slope values. The Max algorithm produces the highest slope values, independent of resolution. It can also be noted that the lowest values produced from Quadsurface algorithm are not statistically different from the rest of its group (3FDWRD, 3FD, 3FDWRSD, and FFD). It is also expected that the algorithms in this group, all estimating slope using horizontal and vertical elevation differences involving 8 or 4 neighbouring cells in the 3*3 window, are not significantly different. The Max algorithm, only considering the steepest slope from the centre cell to one of the neighbours, and the SimpleD including the two cells to the East and South in the slope estimation, both form individual groups. Also, the 2FD algorithm which estimates slope in a simplistic way by involving 4 neighbouring cells on the perpendicular direction forms a group by itself for all resolutions on the flat area.

Fig. 4 (a)-(d) shows that the differences between the algorithms (the span between the minimum value estimated by the algorithm producing the lowest values and the maximum value estimated by the algorithm producing the highest values) decrease when the resolution is changed from finer to coarser. This confirms the trend in the [MS.sub.b] values presented in Table 2.

Also for the steep area (see Fig. 5, (a)-(d)) it is notable that estimated slope values from the eight algorithms are more concentrated and lower when the resolution becomes coarser. At the 5 m resolution, all eight algorithms belong to the same group, confirming the results of the ANOVA. This is logical, since steeper terrain in general diminishes differences between algorithms (a pronounced trend in slope will evenly be captured by e.g. the Max and SimpleD algorithms), and a generalized surface (when decreasing resolution) results in a more pronounced/generalized slope. This effect of generalization is also confirmed by the SDI values presented in Fig. 3. Summarizing the differences between the algorithms for the steep area, it is concluded that the Max algorithm produces the highest slope values for all resolutions, even if it is not significant for the 5 meter DEM. The Quadsurface algorithm produces significantly lower slope values than all other algorithms for all resolutions but 5 meter. It should also be noted that the Max and SimpleD algorithm sensitivity to resolution is much higher than that of the other algorithms, changing more than 2 units when changing resolution from 0.5 to 5 metres. This also confirms the generalisation of the DEM when increasing cell size. It is worth noting that the range of estimated slope values (the length of the horizontal lines in Fig. 5) increases when DEM resolution is decreased. This is also explained by the generalization, where some areas become 'flatter', and other 'steeper.

[FIGURE 5 OMITTED]

4. Discussion

4.1. General characteristics of algorithms

In comparison with previous studies that rank slope algorithms based on 'ground truth' data, our study provides more detailed information regarding the characteristics and behaviours of the eight tested slope algorithms. The analysis of mean slope values (see Table 1), together with the pairwise comparisons presented above, strongly indicates that the Quadsurface algorithm produces the lowest slope values, while the Max algorithm produces the highest values, independent of terrain and resolution. For the Quadsurface method, the surface is fit to all eight surrounding cells, using a least squares of residuals approach, and this results in a reduction of the surface roughness when estimating the slope, consequently producing lower values than other algorithms. The Max algorithm estimates slope based on maximum difference in elevation between the centre cell and the eight surrounding neighbour cells. This exaggerates slope, resulting in higher values than other algorithms. Only involving one neighbouring cell in the estimation also makes the Max algorithm more sensitive to generalisation (i.e. change in resolution). Referring to the comparisons between the standard deviation values of the different algorithms, the SimpleD method produced the highest values, explained by the bias introduced only including the neighbour cells east and south of the centre cell in the estimations. The other six algorithms are based on the two perpendicular partial gradients, and are less diverse in estimating slope. However, the results indicate that a larger number of neighbouring cells included in the estimation (8 instead of 4 or 2) results in relatively low estimated slope, as well as lower sensitivity to generalisation/changes in DEM resolution. This can be explained by the fact that the algorithm itself generalises the surface, by including as many as 2-8 neighbour cells.

4.2. Terrain influence

Two areas with different topographical characteristics were used in this study; one smoother surface (referred to as the flat area in the text) with small absolute differences in elevation and one steeper area with a higher degree of surface roughness (see Fig. 1). The flat surface with smaller variations in slope generally makes the absolute values of differences between the tested slope algorithms smaller than the steep area (see [MS.sub.b] in Table 2), which is also confirmed by previous findings (e.g. Carter 1992). However, the smaller absolute differences between the slope estimates in the flat area are statistically more significant than the ones in the steep area (see F in Table 2). The two most simplified algorithms, Max and SimpleD, differ significantly from the other algorithms both in flat terrain and in high resolution steep terrain. This is because of the 'built-in' generalisations (see above) present in the other algorithms, but not in these two. Independent of terrain form, the Max and SimpleD algorithms seem to overestimate slope compared to the other algorithms. Especially for the steep area with high resolution, the Max algorithm produces much higher estimates than other algorithms. In the steep area, the Quadsurface algorithm seems to underestimate slope. This is due to a more pronounced 'built-in' generalizing effect in the Quadsurface algorithm compared to the 3FDWRD, 3FD, 3FDWRSD, and FFD algorithms, but this is still not big enough to be significant when elevation differences are relatively small (i.e. flat areas).

4.3. Generalisation (DEM resolution)

Regarding the influence of generalization (DEM resolution) on estimated slope values, the results strongly indicate that an increase in cell size results in a decrease in estimated slope values. This also confirms the results by Hasan et al. (2012) and Chang and Tsai (1991), and is valid for all algorithms. Also, the relative differences between different slope algorithms decrease with decreased DEM resolution, (see [MS.sub.b] in Table 2), supporting findings by Thompson et al. (2001). Thompson concluded that certain landscape features are less discernible at the coarser resolution. A coarser resolution generally suppresses surface roughness, and thus represents the terrain in a less complex way, with more gentle slopes. This results in smaller differences between algorithms, where e.g. the number of neighbouring cells included in the estimations is less important. The effect is obvious when applying the different algorithms on a steep surface with a coarse resolution (5 meter, see Fig. 5), where no significant differences between the algorithms are found. Even if the same trend is observed for all algorithms when the resolution is changed, it can be noted that the Max and SimpleD algorithms seem to have a more distinct decrease in estimated slope than the other algorithms. This is due to their relatively higher sensitivity to surface roughness, where a large deviation in elevation of one or two neighbour cell significantly influences the estimated slope value. The Quadsurface algorithm, on the other hand, seems to be less sensitive to cell size in steeper areas. This is explained by the 'built-in' generalisation in the algorithm itself. Referring to sensitivity in both terrain form and generalisation, one should be very cautious if applying either the Max or SimpleD algorithms when estimating slope.

It is also worth noting, referring to the pairwise comparison (see Figs 4 and 5), that the differences between algorithms in the flat area do not change significantly when applied to DEMs with different resolutions, while in the steep area the significant differences decease with an increase in cell size. At the 5 meter resolution no significant differences are found. The more pronounced generalisations in steeper terrain DEMs, with sharp reductions of spatial variation, result in declined differences between different algorithms. During the generalization processes of the DEMs (moving from 0.5 to 5 meter resolution), the flat surface loses less information than the steep one (see Fig. 3).

4.4. Concluding remarks

Based on the above discussions, it is not recommended to use the Max and SimpleD algorithms unless the study or the dataset suffers from 'special circumstances'. For instance, the Max method can be useful when estimating the slope of a channel, because the finite slopes adjacent to the channel will be ignored using this method (Wilson, Gallant 2000). The SimpleD method can, as is the case of the Max algorithm, lead to a relative overestimation of the surface roughness and further increase e.g. the uncertainty of soil erosion modelling (Renard et al. 1991). However, for some applications, e.g. focusing on risks and maximum drainage, Max and SimpleD can be justified. Regarding the Quadsurface algorithm, caution should be paid when implementing it in heterogeneous terrain, mainly due to the smoothing effect producing different results (compared with other algorithms) on steep surfaces.

Apart from the challenges discussed in this study, there are a number of other variables that most probably influence estimated slope. For example, the effects of different interpolation algorithms to create DEMs from scattered LiDAR data, resulting in varying accuracy (Liu 2008; Hasan et al. 2012), have to be taken into account. In this study, the same linear inverse distance interpolation algorithm is implemented throughout, in order to create the different DEMs. The reason for this has been to limit the number of variables, keeping the sources of errors in the DEMs to a consistent minimum, and thus to increase the robustness of the analyses of algorithm differences. Since slope estimation by nature is heavily cell size and terrain dependent, it is desirable to compare the differences between algorithms without a large number of influencing factors.

Other things that can be discussed are the effects of the source of primary elevation data and choice of algorithms. This study is based on high accuracy LiDAR data, making the results applicable also where other less accurate datasets are being used. Regarding the slope algorithms, the selected eight algorithms are commonly used and all based on a 3*3 moving window technique. Since this approach is dominant operationally, the results are relevant, meaningful, and widely applicable to other studies. Even though the analyses and results presented here are based on a case study, they highlight the importance of awareness when choosing slope estimation algorithm. In nature there is no 'true' value of slope, making 'ground truth' comparisons inappropriate. Slope is a matter of scale. Instead of trying to evaluate accuracy, relative performance and sensitivity of different algorithms should be the focal points.

Conclusions

In this study, eight frequently used slope algorithms in flat and steep terrain areas at four levels of resolution are compared. Significant differences are found among the eight algorithms in both flat and steep areas. With reduced resolution of the DEM, the differences among algorithms are decreased, until no statistical differences at the significance level of 0.01 were found at 5 m resolution in the steep area. The Max and SimpleD algorithms always produce higher slope values than the rest, and are not recommended for application, except in special cases. The Quadsurface algorithm, with its strong smoothing effects, always shows relatively lower values, and could easily remove existing roughness and terrain details. One needs to be cautious when applying this algorithm in a steeper area with finer resolution. In general, the choice of slope algorithm at finer resolution becomes more influential for further modeling than at coarse resolutions, due to the greater diversity among algorithms. Overall, the results from this study, illustrating the potential effects from different slope algorithms, could be important for modeling analysis. The results also suggest a need for awareness of the appropriateness of various algorithms' applications at different resolution and terrain forms.

Caption: Fig. 1. LiDAR DEMs with a spatial resolution of 0.5 m covering 100*100 meters areas. On the left is the flat area with less topographic deviation and on the right is the steep area

Caption: Fig. 2. A 3*3 window with numbered cells. The number in the cell is to identify each cell in presented equations, and the [f.sub.x] and [f.sub.y] represent the x and y perpendicular gradients

Caption: Fig. 3. Scaling diversity index for different resolutions of DEM

Caption: Fig. 4. Pair-wise comparisons of eight algorithms at four levels of resolution in the flat area. Flat-R05 represents 0.5 meter resolution in the flat area.

Caption: Fig. 5. Pair-wise comparisons of eight algorithms at four levels of resolution in the steep area. Steep-R05 represents 0.5 meter resolution in the steep area

doi: 10.3846/20296991.2013.806702:

Acknowledgements

Jing Tang was financed by the European Union program Erasmus Mundus 'External Cooperation' (EMECW lot 14). The authors would like to thank Abdulghani Hasan for providing the LiDAR data.

References

Carter, J. R. 1992. The effect of data precision on the calculation of slope and aspect using gridded DEMs, Cartographica: the International Journal for Geographic Information and Geovisualization 29(1): 22-34.

Chang, K.-T.; Tsai, B.-W. 1991. The effect of DEM resolution on slope and aspect mapping, Cartography and Geographic Information Science 18(1): 69-77. http://dx.doi.org/10.1559/152304091783805626

Christensen, T. R.; Johansson, T.; Akerman, H. J.; Mastepanov, M.; Malmer, N.; Friborg, T.; Crill, P.; Svensson, B. H. 2004. Thawing sub-arctic permafrost; effects on vegetation and methane emissions, Geophysical Research Letters 31(4): L04501. http://dx.doi.org/10.1029/2003GL018680

Chu, T. H.; Tsai, T. H. 1995. Comparison of accuracy and algorithms of slope and aspecct measures from DEM, in Proceedings of the GIS AM/FM ASIX95, to 21-24 August, 1995, Bangkok.

Fleming, M. D.; Hoffer, R. M. 1979. Machine processing of landsat mss data and dma topographic data for forest cover type mapping, LARS Technical Report 062879. USA, Purdue University, West Lafayette: Laboratory for Applications of Remote Sensing.

Florinsky, I. V. 1997. Accuracy of local topographic variables derived from digital elevation models, International Journal of Geographical Information Science 12: 47-62. http://dx.doi.org/10.1080/136588198242003

Guentner, A.; Seibert, J.; Uhlenbrook, S. 2004. Modeling spatial patterns of saturated areas; an evaluation of different terrain indices, Water Resources Research 40.

Hasan, A.; Pilesjo, P.; Persson, A. 2012. On generating digital elevation models from liDAR data--resolution versus accuracy and topographic wetness index indices in northern peatlands, Geodesy and Cartography 38(2): 57-69. http://dx.doi.org/10.3846/20296991.2012.702983

Horn, B. K. P. 1981. Hill shading and the reflectance map, in Proceedings of the IEEE 69(1): 14-47. http://dx.doi.org/10.1109/PROC.1981.11918

Irfan Ashraf, M.; Zhao, Z.; Bourque, C. P. A.; Meng, F.-R. 2012. GIS-evaluation of two slope-calculation methods regarding their suitability in slope analysis using high-precision LiDAR digital elevation models, Hydrological Processes 26(8): 1119-1133. http://dx.doi.org/10.1002/hyp.8195

Jones, K. H. 1998. A comparison of algorithms used to compute hill slope as a property of the DEM, Computers and Geosciences 24(4): 315-323. http://dx.doi.org/10.1016/S0098-3004(98)00032-6

Kienzle, S. 2004. The effect of DEM raster resolution on first order, second order and compound terrain derivatives, Transactions in GIS 8(1): 83-111. http://dx.doi.org/10.1111/j.1467-9671.2004.00169.x

Kopecky, M.; Cizkova, S. 2010. Using topographic wetness index in vegetation ecology: does the algorithm matter?, Applied Vegetation Science 13(4): 450-459. http://dx.doi.org/10.1111/j.1654-109X.2010.01083.x

Liu, X. 2008. Airborne lidar for DEM generation; some critical issues, Progress in Physical Geography 32(1): 31-49. http://dx.doi.org/10.1177/0309133308089496

O'Callaghan, J. F.; Mark, D. M. 1984. The extraction of drainage networks from digital elevation data, Computer Vision, Graphics and Image Processing 28: 323-344. http://dx.doi.org/10.1016/S0734-189X(84)80011-0

Olefeldt, D.; Roulet, N. T. 2012. Effects of permafrost and hydrology on the composition and transport of dissolved organic carbon in a subarctic peatland complex, Journal of Geophysical Research G: Biogeosciences 117.

Renard, K. G.; Foster, G. R.; Weesies, G. A.; Porter, J. P. 1991. RUSLE: revised universal soil loss equation, Journal of Soil and Water Conservation 46(1): 30-33.

Ryden, B. E.; Fors, L.; Kostov, L. 1980. Physical properties of the tundra soil-water system at Stordalen, Abisko, Ecological Bulletins 30: 27-54.

Sharpnack, D. A.; Akin, G. 1969. An algorithm for computing slope and aspect from elevations, Photogrammetric Survey 35: 247-248.

Skidmore, A. K. 1989. A comparison of techniques for calculating gradient and aspect from a gridded digital elevation model, International Journal of Geographical Information Systems 3(4): 323-334. http://dx.doi.org/10.1080/02693798908941519

Sorensen, R.; Zinko, U.; Seibert, J. 2006. On the calculation of the topographic wetness index; evaluation of different methods based on field observations, Hydrology and Earth System Sciences (HESS) 10: 101-112. http://dx.doi.org/10.5194/ hess-10-101-2006

Thompson, J. A.; Bell, J. C.; Butler, C. A. 2001. Digital elevation model resolution: effects on terrain attribute calculation and quantitative soil-landscape modeling, Geoderma (1-2) 100: 67-89. http://dx.doi.org/10.1016/S0016-7061(00)00081-1

Travis, M. R.; Elsner, G. H.; Iverson, W. D.; Johnson, C. G. 1975. VIEWIT: computation of seen areas, slope, and aspect for land-use planning, USDA Forest Services General Technical Report PSW-11/1975. Berkeley, California: Pacific Southwest Forest and Range Experimental Station.

Unwin, D. J. 1981. Introductory Spatial Analysis. London and New York: Methuen.

Warren, S. D.; Hohmann, M. G.; Auerswald, K.; Mitasova, H. 2004. An evaluation of methods to determine slope using digital elevation data, Catena 58(3): 215-233. http://dx.doi.org/10.1016/j.catena.2004.05.001

Wilson, J. P.; Gallant, J. C. 2000. Terrain Analysis: Principles and Applications. New York: Wiley.

Wood, J. D. 1996. The geomorphological characterisation of digital elevation models. PhD. University of Leicester.

Woods, R. A.; Sivapalan, M.; Robinson, J. S. 1997. Modeling the spatial variability of subsurface runoff using a topographic index, Water Resources Research 33(5): 1061-1073. http://dx.doi.org/10.1029/97WR00232

Yue, T. X.; Ma, S. N.; Wu, S. X.; Zhan, J. Y. 2007. Comparative analyses of the scaling diversity index and its applicability, International Journal of Remote Sensing 28(7): 1611-1623. http://dx.doi.org/10.1080/01431160600887714

Zhou, Q.; Liu, X. 2004. Analysis of errors of derived slope and aspect related to DEM data properties, Computers and Geo sciences 30(4): 369-378. http://dx.doi.org/10.1016/j.cageo.2003.07.005

Jing TANG is a PhD student in the Department of Physical Geography and Ecosystem Science, Lund University. Her main area of interest is studying distributed hydrological modeling in northern peatland catchments. She dedicates her PhD to develop the hydrological processes in the dynamic ecosystem model, LPJ-GUESS and try to investigate the influences of hydrological processes on the carbon cycling.

Petter PILESJO. Professor Dr. Director of Lund University GIS Centre. His main research interests are hydrological modelling and interdisciplinary applications of geographical information science. Prof. Pilesjo is heavily involved in research and teaching in Africa and Middle East. Migration, land degradation, and SDI are topics on the agenda for Iraq, Uganda, and Tanzania.

Andreas PERSSON. Ph.D. Assistant Professor in Physical Geography and Ecosystem Sciences at Lund University, Sweden. MSc and Ph.D. from Lund University. His research is focused on distributed hydrological modelling with development of new techniques in GIS. Fieldwork in climates ranging from arid to subarctic to apply hydrological models and new techniques in the context of climate change is a major part of his research. Teaching includes hydrology, distributed modelling, GIS and remote sensing.

Jing Tang (1), Petter Pilesjo (2), Andreas Persson (3)

(1, 2, 3) Department of Physical Geography and Ecosystem Science, Lund University, Lund, Sweden

(1) Department of Resource and Environmental Science, Wuhan University, Wuhan, China

E-mail: (2) [email protected] (corresponding author)

Received 11 April 2013; accepted 16 May 2013

Table 1. The mean and standard deviation of estimated slope from
different algorithms in flat and steep area at four resolutions

Selected   Algorithms                  Mean Slope
area
                         R05m       R1m        R2m        R5m

Flat           2FD       3.533      2.104      1.118      0.493
             3FDWRD      2.815      1.502      0.755      0.346
               3FD       2.777      1.475      0.746      0.340
             3FDWRSD     2.874      1.549      0.778      0.357
               FFD       2.810      1.559      0.849      0.368
               MAX       6.700#     4.312#     2.501#     1.039#
           Quadsurface   2.769@     1.474@     0.746@     0.340@
             SimpleD     5.056#     3.479#     2.015#     0.826#

Steep          2FD       24.232     23.783     23.583     23.378
             3FDWRD      23.786     23.596     23.467     23.289
               3FD       23.783     23.594     23.463     23.284
             3FDWRSD     23.805     23.603     23.475     23.296
               FFD       23.914     23.646     23.473     23.281
               MAX       30.327#    27.124#    25.465#    24.341#
           Quadsurface   22.336@    22.188@    22.098@    21.981@
             SimpleD     26.134#    24.446#    23.971#    23.920#

Selected   Algorithms          Standard Deviation (STD)
area
                         R05m       R1m        R2m        R5m

Flat           2FD       2.651      1.623      0.779      0.320
             3FDWRD      2.197      1.202      0.514      0.218
               3FD       2.166      1.177      0.512      0.215
             3FDWRSD     2.240      1.240      0.528      0.225
               FFD       2.138      1.177      0.598      0.232
               MAX       3.943#     2.611#     1.555#     0.518#
           Quadsurface   2.144@     1.173@     0.512@     0.215@
             SimpleD     4.206#     2.680#     1.637#     0.552#

Steep          2FD       12.471     10.871     9.765      8.590
             3FDWRD      11.453     10.229     9.231      8.202
               3FD       11.424     10.196     9.198      8.170
             3FDWRSD     11.516     10.280     9.278      8.241
               FFD       11.602     10.217     9.177      8.102
               MAX       12.318#    10.891#    9.847#     8.178#
           Quadsurface   10.001@    8.841@     7.950@     7.090@
             SimpleD     14.890#    12.368#    10.876#    9.514#

Values in bold represent higher slope estimations, and values
in Italic represent the lowest slope estimations

Note: Values in bold represent higher slope estimations is
indicated with #, and values in Italic represent the lowest
slope estimations is indicated with @.

Table 2. ANOVA testing outputs of eight slope
algorithms in flat and steep areas at four
levels of resolution

Selected    Resolution   Between Groups   Within Groups
area                      ([MS.sub.b])    ([MS.sub.w])

Flat           R05         83100.510          7.995
                R1         11632.140          2.970
                R2          1077.872          0.891
                R5           23.530           0.115
Steep          R05         238357.251        144.735
                R1         19073.693         110.826
                R2          1953.901         89.241
                R5          147.174          69.611

Selected    Resolution   F ([MS.sub.b]/    P-value    df1    df2
area                       [MS.sub.w])

Flat           R05          10393.412         0        7    313624
                R1          3916.947          0        7    76824
                R2          1210.111          0        7    18424
                R5           204.679       4.1E-242    7     2584
Steep          R05          1646.849          0        7    313624
                R1           172.105       6.9E-254    7    76824
                R2           21.895        1.13E-29    7    18424
                R5            2.114         0.039      7     2584

* df1 and df2 are the degrees of freedom of
groups and data points, respectively. Value
in bold represents no significance at the level
of 0.01.