The effect of error in gridded digital elevation models on the estimation of topographic parameters

Environmental Modelling & Software 21 (2006) 710e732 www.elsevier.com/locate/envsoft

The effect of error in gridded digital elevation models on the estimation of topographic parameters Lynn D. Raaflaub, Michael J. Collins * Department of Geomatics Engineering, University of Calgary, Calgary, AB, Canada T2N 1N4 Received 21 May 2004; received in revised form 12 December 2004; accepted 16 February 2005 Available online 27 April 2005

Abstract Digital elevation models (DEMs) provide the basic information required to characterise the topographic attributes of terrain. The primary derived topographic parameters associated with DEMs are slope and aspect. Slope and aspect maps are used in a wide variety of applications. Slope and aspect can be used to calculate other significant topographic parameters such as upslope area and topographic index. The topographic index, in turn, can be used by distributed hydrological models to characterise the spatial distribution of terrain moisture. Many algorithms have been developed to calculate slope, aspect and upslope area from DEMs e specifically from gridded DEMs e but little work has gone into determining the uncertainty in these parameters, or the affect of this uncertainty in further applications. The accuracy of these parameters is dependent both on the algorithm and on the errors associated with the DEM itself. Since it is almost impossible to model all the errors associated with a given slope/aspect algorithm and since a DEM is normally only provided with a single rms error, simple error propagation is not adequate to determine the error associated with the derived topographic parameters. A more rigorous method of determining the affect of DEM errors on derived topographic parameters is with statistical analysis using Monte Carlo simulation and error realisations of the DEMs. In this research we demonstrate that the error sensitivity of slope decreases as the number of neighbours used in the algorithm increases, hence steepest neighbour algorithms, which are common in hydrology are more sensitive to DEM error than algorithms that use four or more neighbours. In contrast, the average error sensitivity of aspect to DEM error is not dependent on the algorithm used. However, while the mean variability of this sensitivity was lower for the steepest neighbour algorithms, their errors were spread over a greater variety of slopes while the eight neighbour algorithms had errors confined to flat regions. The error sensitivity of upslope area and topographic index is related to the use of steepest neighbour flow routing algorithm. Ó 2005 Elsevier Ltd. All rights reserved. Keywords: Digital elevation model; Topography; Error analysis

1. Introduction Digital elevation models (DEMs) are representations of the terrain elevation as a function of geographic location. As such, they provide the basic information required to characterise the topographic attributes of

* Corresponding author. Tel.: C1 403 220 4952; fax: C1 403 284 3697. E-mail address: [email protected] (M.J. Collins). 1364-8152/$ - see front matter Ó 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.envsoft.2005.02.003

terrain. DEMs are typically represented in two formats: contour maps, where the surface is represented by lines of constant elevation at even intervals; or point heights, where the surface elevation is sampled on either regular or irregular bases. The choice of DEM format is dependent on the application and, of course, on the availability of data. The main topographic parameters derived from DEMs are slope and aspect. These parameters are local descriptions of the downhill slope, in the sense that they are the magnitude and direction of the vector tangent to

L.D. Raaflaub, M.J. Collins / Environmental Modelling & Software 21 (2006) 710e732

the downhill. While the definitions of slope and aspect are fairly straightforward, the actual calculation of these variables is not. Difficulties arise in trying to model a specific element of a terrain surface from a model of a terrain surface. Because of the problems associated with deriving slope and aspect from DEMs, many algorithms have been developed to accomplish this goal. Topographic information such as slope and aspect have wide applicability in the environmental sciences. Since the spatial routing of water over terrain is strongly dependent on the shape of the downhill surface, one such application is distributed hydrological modelling. Hydrologically, slope is an indication of the amount of gravitational energy available to drive water flow, hence it influences the rate of water flow. Aspect, which defines the slope direction, can be used to determine the direction of water flow, which is the information required to determine other hydrologically significant variables such as upslope area (Zevenbergen and Thorne, 1987). Upslope area is the total area upslope of a given location, which is an indication of how much water can flow through the location. From a hydrological standpoint, the effectiveness of a slope and especially an aspect calculation is determined by how accurately upslope area can be resolved. Like slope and aspect, there are different algorithms available for the calculation of upslope area, which are essentially dependent on the choice of algorithm used to calculate slope and aspect. The slope, along with the upslope area (A), can be used to calculate a parameter known as the topographic index (ln(A/slope)). This parameter reflects the spatial distribution of soil moisture, surface saturation and runoff generation processes, hence it is one of the parameters used in distributed hydrological modelling. TOPMODEL is one such distributed hydrological model that uses surface topography, in the form of the topographic index, along with soil properties to predict the distribution of soil moisture and runoff (Beven and Kirkby, 1979). It uses the assumption that the spatial distribution of the topographic index approximates the spatial distribution of the depth to the water table in a watershed. Hence any variations in the calculation of the topographic index will affect the outcome of the model’s predictions. If all sources of error associated with a DEM and with the algorithms used to derive the topographic parameters were known, the propagation of error through the topographic models and subsequent models utilising the topographic parameters could be done using analytical error propagation. The difficulty arises in the errors associated with the assumptions inherent to the algorithms of the topographic parameters. It is almost impossible to model all the errors associated with a given algorithm, therefore a simple error propagation is not adequate to represent the error of values derived

711

from a DEM with associated error. The affect of error in a DEM on its derived parameters is best determined with statistical analysis using Monte Carlo stimulation and error realisations of the DEM (Veregin, 1997). 1.1. Research objectives Our objective was to examine the affect of model and data uncertainty on the estimation of topographic parameters. The affect of both correlated and uncorrelated gridded DEM data error will be investigated using Monte Carlo methods. Uncorrelated error is seen as the worst case scenario, where completely random error is applied to each elevation in the DEM. The more realistic correlated error, on the other hand, is smoother due to the high correlation between neighbouring error values. We hypothesize that correlated error should produce less variation in the topographic parameters compared with uncorrelated error. As mentioned above, many slope and aspect algorithms have been developed to calculate topographic parameters from DEMs, specifically from gridded DEMs. These algorithms are generally based on neighbourhood operations. The primary difference between the algorithms is the number of neighbours used in the calculation. For a three by three cell neighbourhood, between two and nine grid cells may be used. In order to investigate how model uncertainty affects topographic analysis and hydrological models, nine slope and aspect algorithms will be examined. This is in essence performing an error sensitivity analysis on the various algorithms presented. The more neighbours used the more constrained the algorithms become. Therefore, it is hypothesised that the more neighbours utilised in the algorithm the less affect error will have on the derived parameters and hence on applications that use the parameters. To accomplish the objectives presented, this paper has been organised in the following way: Section 2 reviews current algorithms for estimating slope, aspect, upslope area and topographic index, Section 3 describes the study area, the data set used and presents the experimental methodology and implementation, Section 4 presents the experimental results and finally Section 5 provides conclusions and recommendations.

2. Theory 2.1. Surface representation by gridded elevations A DEM may be defined as any numeric or digital representation of terrain elevation given as a function of geographic location. As mentioned earlier, the type of DEM that is used in a given application depends on the requirements of the application and on the availability of the DEM. With the extensive use of computers in

712

L.D. Raaflaub, M.J. Collins / Environmental Modelling & Software 21 (2006) 710e732

modelling, DEMs that can be easily stored and manipulated in digital form are preferable. Discrete DEMs sampled on a regular interval, known as gridded DEMs, have these advantages, and hence are utilised in the majority of environmental models. Although there have been a few algorithms published for analysing other forms of DEM (Goodrich et al., 1991; Moore et al., 1988), their use is relatively rare. In gridded DEMs, elevations, z, relative to some data are posted at regular planimetric intervals. These regular intervals may be orientated in any direction, but usually take the form of an eastewest, x, and northesouth, y, representation. The distance between two samples is known as the grid spacing. Usually, this distance is the same for both grid directions. The elevations in a DEM normally represent the actual height of the ground at the location where they are given. A gridded representation allows for easy access and manipulation since the elevations can be stored as a simple matrix. Another advantage to using gridded DEMs as opposed to other terrain models is that they are normally easy to obtain in digital form for a given area without the need of further manipulation. Because of the fixed sample spacing, gridded DEMs suffer from sampling problems. Terrain is usually undersampled in rough (rapidly fluctuating elevation) areas and oversampled in flat areas, and a balance must be found between the sampling interval and other requirements. The choice of grid spacing is not normally available to the user of DEMs, therefore, when looking at DEM uncertainty the issue of possible undersampling should be considered. Though gridded DEMs have their advantages, there are a few things that need to be taken into account when using them. Gridded DEM data may be acquired in a variety of ways. These include measurement of the terrain elevation directly using land surveying or indirectly using photogrammetry (or perhaps interferometric radar or inSAR), or interpolation from an existing gridded DEM or derivation from either a digital or paper contour map. It is important to recognise that if a DEM of one type is generated from a DEM of another type, instead of from measurements, it will suffer from compounded modelling error. Regardless of how a DEM is generated, it is important to remember that what is being produced is a discrete representation of a continuous surface. In other words, even the highest quality DEM is an approximation of the surface and not an exact representation. One of the most important characteristics of a gridded DEM is the grid spacing. Another term that is occasionally and incorrectly interchanged with grid spacing is resolution. This terminology stems from such fields as image processing, where a grid value, or pixel, represents an area average, derived from the impulse response of the imaging system. Similarly, a gridded DEM will only have a resolution if the elevation values

are grid cell averages. For example, when a DEM has been derived from satellite imagery or inSAR data, each elevation will represent an area average. However, in many if not most, gridded DEMs the elevation values represent point elevations, and consequently the term, resolution, is not part of the metadata. In any event, the use of resolution in place of grid spacing is always incorrect. It should be further noted that some topographic algorithms treat grid values as area averages. 2.2. Pit removal The generation of a DEM invariably introduces a number of ‘‘artificial’’ topographic features that should be detected and corrected. Hydrologically, the most serious of these features are pits (sink features) (O’Callaghan and Mark, 1984; Band, 1986; Jenson and Dominque, 1988), and to a lesser extent, dam features (Quinn et al., 1991). Dam features occur when the elevation values are artificially raised above the actual elevation value, while pits occur at points that do not have any neighbours with a lower elevation. While it is possible for pits to be found in topographic surfaces, especially in recently glaciated areas or in limestone areas, most pits can be considered errors since fluvial erosion processes will not normally produce such features at the scale resolved by DEMs (Band, 1986). Pits generally appear in flatter areas where even a 1-m error in elevation can be enough to produce a closed depression, while on steeper slopes a higher variation would be required. Artificial pits in a DEM will cause serious problems for any subsequent algorithm that depends on mapping hydrologically connected regions, such as in the calculation of upslope area. Because of the errors that pits can introduce, their removal is an important first step when using DEMs for hydrological applications. The method used for pit removal can influence the effectiveness of other algorithms applied to the DEM. There are a variety of techniques to perform pit removal e examples of their use can be found in O’Callaghan and Mark (1984), Band (1986), Martz and de Jong (1988), Jenson and Dominque (1988), Hutchinson (1989), Fairfield and Leymarie (1991), Freeman (1991), Tribe (1992), Tarboton (1997). In the simplest method, the pits are removed manually by flagging them in a preprocessing stage. Unfortunately, this process is tedious and labour intensive, and consequently a more automated approach is preferred. The most straightforward method of automated pit removal is to numerically smooth the terrain. This method reduces the number of pits, especially shallow ones, but it also tends to oversmooth the topography, which may lead to a significant loss of information. An alternative approach to pit-removal attempts is to preserve flow paths and drainage directions. These other methods primarily involve finding

L.D. Raaflaub, M.J. Collins / Environmental Modelling & Software 21 (2006) 710e732

a pit’s outflow point and directing flow out through this point. This can be accomplished by raising the elevation of the pit to that of the outflow point, or by redirecting flow lines to neighbouring cells that have defined flow directions that do not point back into the pit itself. Several of the pit removal techniques outlined above require knowledge of the surface’s aspect, and, by extension, its slope. Additionally, the calculation of upslope area also requires aspect. Like pit removal, there are several techniques to calculate these quantities. In the following section, several of these techniques are outlined. 2.3. Slope and aspect algorithms Traditionally, slope and aspect are calculated manually from contour data. In this method a line tangent to the contours is drawn along with its corresponding perpendicular bisector. The slope is found by dividing the difference in height along the perpendicular bisector by the length of the perpendicular bisector. Even though this method has been used to represent the ‘true’ slope and aspect values of a terrain (Skidmore, 1989), it is almost impossible to obtain the same values for slope and aspect twice for any given location due to the manual calculation. It is not surprising then that other methods have been developed that depend solely on digital data and yielding repeatable results. There are many different ways in which slope and aspect can be calculated from gridded DEMs. Generally, their determination is based on neighbourhood operations where calculations are made for a cell based on the values of the cells that are spatially adjacent in the grid. Some methods are better than others in representing the actual slope and aspect of the terrain, while others are more suited to hydrological analysis. As mentioned previously, slope and aspect are just the magnitude and direction of the vector tangent to the topographic surface. The goal then is to create a localised model of the topographic surface from which slope and aspect values can be derived. Differences in slope and aspect algorithms have to do with how these localised surface models are defined. These definitions are based on the number of neighbouring points used, whether or not the centre point contributes, and whether or not the aspect calculation will be restricted to a finite number of directions (Guth, 1995). In all slope/aspect algorithms, except those utilising the steepest neighbours, some form of numerical differentiation is employed. In terms of a topographic surface, the numerical derivative is calculated from two quantities: the elevation difference and the ground distance (Eyton, 1991). This is an example of differentiation using finite difference calculus, where the first derivative of elevation describes the rate of change of elevation, which is the slope (Eyton, 1991; Dozier and

713

Strahler, 1983; Horn, 1981). Together, the slope in the x direction and the slope in the y direction (the partial derivatives of z with respect to the x and y directions), define the gradient vector of the surface. The maximum slope can be determined by taking the norm of this vector, sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2  2ffi vz vz slopeZ : C vx vy

ð1Þ

The aspect, which is the direction of the maximum slope, is the angle between the slope defined in x and the slope defined in y, which is given by the relationship aspectZtan
1



 vz=vx : vz=vy

ð2Þ

The difference in the slope and aspect algorithms which use Eqs. (1) and (2) are related to how they determine the local slope in x (vz/vx) and y (vz/vy). The determination of the local slope in x and y is dependent on one to eight of the neighbouring grid points. As a consequence, consideration must be made for those grid points that lie on the edge of the DEM, which do not have eight neighbours. In order to make the calculations consistent across the entire DEM and to avoid edge effects, the slope and aspect values are not normally calculated at the edges. In essence, this creates slope and aspect maps that are slightly smaller than the original DEM. However, this is preferable to the alternative of obtaining erroneous edge values. Therefore, all algorithms presented are for internal grid points. Since the algorithms are all detailed in the literature we will not reproduce them here. In order to emphasise the differences between algorithms, we present a graphical example of a single calculation in Fig. 1. In the local neighbourhoods illustrated, the cells used in the calculation of slope and aspect are highlighted in grey, while the dark blue arrow depicts the resulting aspect (flow direction) from the centre cell. For comparison, the numerical values for the slope and aspect are also provided. The algorithms used are as follows. 2.3.1. Three point plane (TPP) O’Neill and Mark (1987) defined a localised topographic surface using a three point plane. This method defined a triangle using three grid points: the point for which the slope and aspect value is to be calculated, and its neighbouring points to the east and north. The choice of using the eastern and northern neighbouring grid points is arbitrary, since the DEM could be orientated in any direction.

714

L.D. Raaflaub, M.J. Collins / Environmental Modelling & Software 21 (2006) 710e732

TPP

FCN

965 957 955

965 957 955

953 950 947

953 950 947

945 942 937

945 942 937

slope = 0.305 aspect = 157o

slope = 0.323 aspect = 158o

FDN

ENU

965 957 955

965 957 955

953 950 947

953 950 947

945 942 937

945 942 937

slope = 0.420 aspect = 155o

OODS

slope = 0.357 aspect = 171o

MDG

965 957 955

965 957 955

953 950 947

953 950 947

945 942 937

945 942 937

slope = 0.347 aspect = 168o

slope = 0.368 aspect = 135o

MDN

MAG

965 957 955

965 957 955

953 950 947

953 950 947

945 942 937

945 942 937

slope = 0.237 aspect = 158o

slope = 0.424 aspect = 135o

Fig. 1. Examples of slope and aspect algorithms investigated in this paper.

2.3.2. Four closest neighbours (FCN) Uses the four cardinal neighbours, those to the north, south, east and west, representing a second order finite difference relationship (Hoffer et al., 1979; Unwin, 1981; Sharpnack and Akin, 1969; Ritter, 1987; Eyton, 1991; Dozier and Strahler, 1983; Skidmore, 1989; Guth, 1995; Hodgson, 1995; Jones, 1998). The elevations at these four closest neighbours can be used to define two orthogonal components of slope, the slope in x and y, which define the steepness and downhill direction, but does not define a surface at the point of interest (Guth, 1995). 2.3.3. Four diagonal neighbours (FDN) A variation to the four closest neighbours algorithm is the four diagonal neighbours algorithm (Jones, 1998). The difference in this algorithm is that the perpendicular gradients are taken at an angle of 45  to the principal axis of the DEM grid. By using the diagonal neighbours, those points to the northwest, northeast, southwest and southeast of the central point, along with a correction for the increased distance to the diagonal neighbours, allows for investigation into the possible effects of grid bias. 2.3.4. Eight neighbours unweighted (ENU) When more than four grid points are used, it is possible to define a topographic surface from which slope and aspect values can be determined. One such algorithm uses the eight neighbouring grid points to calculate the slope and aspect, representing a third order finite difference algorithm (Sharpnack and Akin, 1969; Papo and Gelbman, 1984). The slope and aspect values obtained using this method are the same as those for a least squares linear trend surface fitted to these eight neighbouring grid points. 2.3.5. One over distance squared (OODS) Horn (1981) implemented a slope and aspect algorithm in which grid point neighbours were weighted proportional to the reciprocal of the square of the distance from the central point. 2.3.6. Maximum downhill gradient (MDG) A model introduced by O’Callaghan and Mark (1984), finds the steepest downhill neighbour in order to determine topographic parameters such as slope, aspect, upslope area and drainage networks (Band, 1986; Jenson and Dominque, 1988; Fairfield and Leymarie, 1991). This method is based on an eight neighbour connectivity where the drainage direction (aspect) of the central point is taken as the direction to the neighbour with the maximum drop, i.e. the neighbour which is lower than the central point and is the minimum value. Problems arise when more than one neighbour fits the criteria

L.D. Raaflaub, M.J. Collins / Environmental Modelling & Software 21 (2006) 710e732

simultaneously. There are, of course, different ways to deal with this problem, some of the more popular include arbitrarily choosing the first neighbour encountered when going through the neighbours clockwise from north, or flagging the point for later inspection. Problems are encountered when the central point is a pit and no such minimum neighbour can be found. Specific implementations of this algorithm vary in terms of how they handle this special case. In this study, the first instance of the lowest neighbour was selected when going through the neighbours clockwise from north. 2.3.7. Multiple downhill neighbours (MDN) One of the main disadvantages of the steepest downhill neighbour algorithm is that it is greatly affected by grid orientation bias. Truly accurate values for slope and aspect are not obtained because of the need to follow the major points on the compass. Multiple downhill neighbours’ methods (Quinn et al., 1991; Freeman, 1991; Wolock and McCable, 1995) are a variation of the steepest downhill neighbour algorithm which attempt to solve the limitations present when dealing with a one dimensional representation of flow. This is accomplished by distributing flow from a pixel amongst all of its lower elevation neighbour pixels. Slope and aspect are calculated by taking the average of all the downslope neighbours, i.e. the average slope and aspect values. 2.3.8. Maximum adjacent gradient (MAG) Another variation on the maximum downhill gradient algorithm is the maximum adjacent gradient algorithm (Skidmore, 1989; Guth, 1995). In this algorithm, the steepest slope is found among the eight neighbours, compensating for the distance to the neighbour cell, and assigned to the slope value. If the slope is uphill, the downhill aspect is considered to be in the directly opposite direction. This algorithm uses the same equations applied in the MDG algorithm, and hence suffers from the same limitations. 2.3.9. Comparison of slope and aspect algorithms Different slope and aspect algorithms will, of course, result in different slope and aspect values. Guth (1995), found that the choice of algorithm can vary the average slope by as much as 25%. Such a large variation emphasises the care that should be taken in selecting an algorithm that will produce an accurate representation of the terrain. The criteria for correctness need to be based on the desired application. Solving for basic drainage networks could effectively be done using the steepest neighbour algorithms, but if something closer to the actual aspect value is required then the eight neighbours unweighted algorithm might be more applicable, since it provides results which more closely represent the terrain.

715

Comparing the ‘‘correctness’’ of slope and aspect algorithms can be done in a variety of ways. One method is to compare results for a surface with known slope and aspect values (Hodgson, 1995; Jones, 1998). The only way that known slope and aspect maps can be generated is if the terrain surface used for comparison is analytic, e.g. it has been generated by a large term polynomial surface. Given the equation for the surface, actual values for the slope and aspect can be determined. Another method of comparison is to choose a method of determining slope and aspect that is considered to produce the ‘‘true’’ value, such as hand derived values from contour maps (Skidmore, 1989). If no truth values are available, the last alternative is to do statistical comparisons of the results obtained by the differing algorithms (Guth, 1995). However, a statistical comparison between methods really only investigates the variation between methods and cannot be used to compare estimates with ground measurements. Based on the various comparisons done on slope and aspect algorithms (Guth, 1995; Skidmore, 1989; Hodgson, 1995; Jones, 1998), a few generalisations on the algorithms can be made. Surprisingly, algorithms utilising only four neighbouring cell values tend to be better at estimating slope and aspect values than algorithms using eight neighbouring cell values. Algorithms that use only one to two of the neighbouring cell values along with the centre value do not, in general, perform well. These two to three point algorithms normally perform significantly worse than the four closest neighbour algorithm, and should not be used when accurate slope and aspect values are desired. 2.4. Upslope area Upslope area is an important parameter in determining hydrologically significant terrain characteristics such as drainage networks. It represents not only the flow direction of water, but also the accumulated area draining through a point. As mentioned previously, the calculation of upslope area is dependent on the determination of aspect. However, slope and aspect values can be determined using localised neighbourhood operations, upslope area calculations are fundamentally dependent on the entire DEM. Fortunately, they can be determined through the use of modified neighbourhood operations where values are passed down cell by cell. For a gridded DEM, the upslope area can be generalised as the number of cells that drain into a specified cell multiplied by the area of a grid cell. For a given grid cell, the upslope area can be expressed as, upslope areaZNDd 2 ;

ð3Þ

where N is the number of upslope cells and Dd is the grid spacing. Based on this relationship, a ridge point, i.e.

716

L.D. Raaflaub, M.J. Collins / Environmental Modelling & Software 21 (2006) 710e732

a point that has no drainage inputs, would have zero upslope area, while a pit could have a very high upslope area since it can be drained into from all sides. The estimation of upslope area for each cell in a DEM can be difficult since it involves tracing flow paths through the entire DEM. While there are different ways to attack this problem, one of the simplest computationally is a recursive algorithm (Tarboton, 1997; Freeman, 1991). The first step in this process is, of course, to determine the two dimensional flow direction of each cell in the DEM. Besides the aspect grid, two other grids are required, an inflow count grid and an upslope count grid. The cells in the inflow count grid contain a count of the number of cells that flow into them, while the cells in the upslope count grid contain a count of the number of upslope cells. To determine the inflow count grid, each cell is initialised to zero. Then, using the aspect grid, the values of the cells in the inflow count grid are incremented for every cell that drains into them, in other words, every time the flow is directed towards the cell. In this way, the values of the inflow count grid can range from zero to eight (none to all of the neighbours). For the upslope count grid, the cells are initialised to one. The inflow count grid is then gone through until a zero value is found. This indicates that a cell that has no upslope neighbours. The count of the upslope count grid cell (or cells) that the zero valued cell flows into is then increased by adding the zero cell’s upslope area (or portion of). This flow direction is determined using the grid of aspect values. The inflow grid cells involved, including the zero valued cell, all have their values decreased by one. This process is repeated, adding cell outflows iteratively to lower neighbours, until all the cells in the inflow count grid are negative. The values in the upslope count grid indicate how many cells are upslope of a given cell. The upslope area for each cell is found by multiplying the number of upslope cells by the area of each grid cell (Eq. 3). While the procedure outlined above is the general method used to determine upslope area, the calculation of upslope area, like the calculation of slope and aspect, can be determined in different ways. The main difference between these methods is how they deal with flow divisions, i.e. flow can either be restricted to exactly one downslope neighbour, or divided according to either the number of downslope neighbours or the aspect direction. The choice of how to divide flow depends on how the flow direction (aspect) is estimated. Any of the methods can be used to determine the flow direction of a cell, though some algorithms are more suited to specific upslope area algorithms than others. Algorithms used to determine upslope area are presented in the following sections along with the aspect algorithms they best compliment. Also included are some of the algorithms’ strengths and weakness.

2.4.1. Single flow e steepest downhill neighbour The steepest downhill neighbour upslope area algorithm (O’Callaghan and Mark, 1984) goes with the MDG algorithm for determining slope and aspect, although it can also be used with the MAG algorithm (Band, 1986; Jenson and Dominque, 1988; Fairfield and Leymarie, 1991). In this method flow can only be routed to one downslope neighbour. The flow direction, as with the MDG and MAG algorithms, is restricted to one of the cardinal or diagonal directions. This allows flow to be traced with very little dispersion, hence making the determination of drainage networks much easier. Because of the simplicity of this method and the lack of dispersion, this is generally the choice of upslope algorithms from a hydrological standpoint. Despite the popularity of the steepest downhill neighbour upslope area algorithm, many authors have drawn attention to its drawbacks (O’Callaghan and Mark, 1984; Band, 1986; Fairfield and Leymarie, 1991; Costa-Cabral and Burges, 1994; Skidmore, 1989; Guth, 1995; Tarboton, 1997; Wolock and McCable, 1995). The main problem with this method is that it suffers from grid orientation bias. This can dramatically throw a drainage path off course by as much as 44  . In the extraction of drainage networks, this discretisation can prevent stream formation due to the generation of parallel flow lines especially in relatively flat areas. The rougher and steeper the terrain, however, the less significant this problem is. The limitations of this method have prompted the development of other algorithms.

2.4.2. Stochastic method Fairfield and Leymarie (1991) applied a stochastic version of steepest downhill neighbour upslope area algorithm in which the drainage direction was chosen probabilistically. Using their algorithm, flow direction is randomly assigned to one of the downslope neighbours, with the probability proportional to slope. In this method the greatest downhill slope is not always taken in order to more faithfully follow the true aspect of the land. This algorithm attempts to alleviate some of the problems of grid bias, where the flow path must follow one of the cardinal or diagonal direction resulting arbitrarily from grid orientation. While this method does provide more appropriate flow path directions, it introduces problems of its own (CostaCabral and Burges, 1994). Randomness does not ensure reproducible results, making it difficult to get comparable values. There are also problems in flat areas where flow should be parallel. In the stochastic method adjacent flow paths are not parallel but ‘‘wiggle’’ randomly often converging laterally with one another. Once converged, flow paths are nearly impossible to separate, hence flow becomes increasingly concentrated downslope. Because of the difficulties with achieving

L.D. Raaflaub, M.J. Collins / Environmental Modelling & Software 21 (2006) 710e732

reproducible results, the stochastic method was not implemented in this study. 2.4.3. Multiple flow e steepest downhill neighbours The multiple downhill neighbours (MDN) algorithm for slope and aspect is utilised in multiple flow upslope area algorithms. In this method, flow directions are still restricted to the cardinal or diagonal grid cells, but unlike the single flow method, flow is partitioned to all the downslope neighbours. The amount of flow partitioned to each downslope neighbour is based on the size of the slope. The larger the slope to a particular cell the more of the flow that is partitioned to it. For a grid cell with an upslope area of A and n downslope neighbours, the upslope area partitioned (Ap) to each downslope cell i is given by, slopei Api ZA P : n slopei

ð4Þ

1

From a drainage network perspective, the multiple flow upslope area algorithm can significantly disperse the water flow depending on the type of terrain, which can be problematic (Skidmore, 1989). The multiple direction method tends to perform poorer in flat areas than the single direction approach, but can produce more realistic flow paths in steeper terrain (Quinn et al., 1991). While it is possible to more accurately resolve the flow path in a DEM using a multiple flow direction algorithm, the significant dispersion introduced does not make this method a popular choice in some terrains. 2.4.4. Power method In an attempt to compensate for the problems associated with the single and multiple flow algorithms, Holmgren (1994) applied a power weighting scheme to the multiple flow algorithm (Eq. 4) such that x

ðslopei Þ Api ZA P : n ðslopei Þx

ð5Þ

1

The exponent x is used to narrow or widen the flow path of water to obtain the best possible drainage network for a watershed. When the exponent is equal to one, the algorithm degrades into the multiple flow algorithm, and as the exponent gets larger it is essentially the same as the single flow algorithm. Therefore, the single and multiple flow algorithms can be seen as the two extremes of the power method. While it has been shown that the power method, when properly calibrated for a specific watershed, is better at accurately delineated the stream network (Holmgren, 1994; Quinn et al., 1995), the objective of this

717

study is to investigate the effect of DEM error on different algorithms, not to find the algorithm best suited for stream network delineation. Since the value of x will change depending on the watershed, and since this study was only interested in generalised methods, the power method of determining upslope area was not implemented. 2.4.5. Plane fitting Plane fitting uses the aspect associated with each pixel to specify flow direction (Tarboton, 1997), and is conducive to any of the slope and aspect algorithms in which the aspect does not necessarily point directly to one of the eight neighbours (including TPP, FCN, FDN, ENU and OODS). Flow is routed as though it were a ball rolling on a plane released from the centre of each grid cell. Flow in this algorithm is not restricted to the cardinal or diagonal directions. Therefore, this method has the advantage of being able to, in a sense, specify flow continuously. If the aspect angle does not lie directly on one of the cardinal or diagonal directions, flow is partitioned, based on the angle, between two cells. Problems encountered with this method are normally related to discontinuous representations of the surface at pixel edges that lead to inconsistent or counterintuitive flow. 2.5. The topographic index The topographic index is used to represent the soil saturation (Beven and Kirkby, 1979; Wolock and McCable, 1995; Quinn et al., 1995). It is an estimate of the accumulated water flow at any point in a watershed. A high topographic index indicates that the ground in a region has a higher potential to be saturated, while a low topographic index indicates the ground is less likely to be saturated. Areas of saturation are the source areas for overland flow. Generally, high values of the topographic index occur in flat areas with little slope and/or in regions with a large upslope area. The topographic index is generally defined as the natural log of upslope area divided by the slope (ln(A/ slope)). This is, however, somewhat of an oversimplification. More precisely, it is the natural log of the specific catchment area (a) divided by the slope, 

 a topographic indexZln : tanb

ð6Þ

By convention, in the equation for the topographic index the slope is represented by tan b, where b is the slope angle. The slope is used to approximate, under steady state conditions, the local hydraulic gradient, i.e. the slope of the water table.

718

L.D. Raaflaub, M.J. Collins / Environmental Modelling & Software 21 (2006) 710e732

The specific catchment area is defined as the upslope area per contour length, A aZ ; l

ð7Þ

where l is the contour length, and, as before, A is the upslope area. Another way of describing the specific catchment area is to say that it is the area above a certain contour length that contributes flow across the contour (Gallant and Wilson, 2000). Fig. 2 illustrates such an area. For gridded data, the contour length can be given a value equivalent to the grid spacing. The topographic index is calculated for each cell in a DEM. To do this, the slope, aspect and upslope area must first be calculated using, for example, one of the algorithms discussed in Section 2.3. The tan b term is equivalent to the slope estimated by these algorithms. 3. Methods 3.1. Study area and data set The DEM used in this study was taken from a series of gridded DEMs provided by the province of Alberta, Canada. These DEMs have extents that coincide with provincial 1:20000 scale map sheets. Depending on the roughness of the terrain, the DEMs were originally derived from aerial photographs at spacings of 25 m, 50 m, or 100 m. They were then interpolated to 25 m grids by linear least squares prediction using the SCOP program (Ackermann, 1992). For the northing and easting coordinates, the resulting DEMs have a specified rms accuracy of G10 m 90% of the time. The elevations have a specified rms accuracy of G3 m 90% of the time. The remaining 10% not within this bound has a specified accuracy of G5 m 90% of the time. For this study, the error in planimetry was ignored, while the rms error in the elevation was taken to be 3 m. From the available DEMs a watershed was extracted that served as the study area. The watershed selected, Marmot Creek main stem near Seebe, was drawn from the DEM that coincided with the provincial map sheet

Fig. 2. Specific area (a): the upslope area per contour length (a Z A/l ).

C82J14NE. The region covered by this map sheet is located in the Kananaskis region of Alberta, in the front ranges of the southern Canadian Rockies. The watershed, depicted in Fig. 3, covers an area of approximately 10 km2 (15 484 grid cells). Elevations in the watershed DEM ranged between 1584.904 m and 2837.548 m. The stream outlet gauge, shown as a blue dot in Fig. 3, is located at UTM coordinates (5645700 N, 629825 E) and is at an elevation of 1585 m. Extraction of the watershed from the DEM was performed using watershed selection tools available in the GIS software Arc/Info. 3.2. Experimental methodology and implementation 3.2.1. Propagation of error DEMs have errors in both planimetry and elevation. In the algorithms for slope and aspect e the values upon which other topographic parameters are based e only the grid spacing and elevation values are required. If the grid spacing is constant then errors in planimetry may be wrapped into elevation errors. Systematic errors in planimetry may still exist; however, they would probably be in the form of a translation error, where the entire grid has been shifted from its true location, or a scale error, where the grid interval is slightly smaller or larger than the specified interval. We do not treat systematic errors in this work. There are a variety of methods that can be used to propagate error. These methods include an analytical approach, Taylor series expansions, and Monte Carlo simulation (Heuvelink, 1998). In the analytical approach, the mean and standard deviation of a function is explicitly derived. This is dependent on knowing the exact form of the function involved. In methods based on Taylor series expansions the function is approximated by a truncated Taylor series from which the mean and standard deviation are determined. Neither of these methods of error propagation are suitable to the error analysis of this study. While it is possible to use these methods to determine the error in some of the topographic parameters, they cannot be used for all those investigated. Problems arise with the slope and aspect algorithms that are based on the steepest neighbour principle (MDG, MAG and MDN). These algorithms are decision-based, involving the selection of neighbours based on specific properties. Such algorithms do not have the functional form required by the analytical and Taylor series approaches. Another limitation of these approaches is that they are restricted to neighbourhood operations. They are not suitable for global operations, such as the calculation of upslope area, and, by extension, the calculation of the topographic index (Heuvelink, 1998). Also, because of complex model structure and dynamics, these error propagation methods are difficult to implement for a hydrological model. Therefore, the affect of error in

L.D. Raaflaub, M.J. Collins / Environmental Modelling & Software 21 (2006) 710e732

719

Fig. 3. DEM of the Marmot Creek Main Stem watershed in Kananaskis, Alberta. The blue dot is the location of the stream outlet gauge. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

a DEM on its derived parameters and on hydrological modelling is best determined with statistical analysis using Monte Carlo simulation and error realisations of the DEM (Veregin, 1997). 3.2.1.1. Monte carlo method. In the Monte Carlo method of error analysis, a given algorithm is repeatedly applied to sets of random input data. Each set of random input data is termed a realisation, and all the realisations are generated from the same error distribution. The realisations represent possible error-free data sets, and not just a data set with an added error component. Using the results from all the realisations, it is possible to calculate a mean and standard deviation that correspond to the algorithm under analysis. The mean is considered the ‘true’ value that the algorithm would produce, while the standard deviation is considered the error. Because the analysis of results from a Monte Carlo simulation is done using statistical, and not analytical, quantities, the method is said to use a stochastic approach to error propagation (Heuvelink, 1998). A simple Monte Carlo simulation of DEM error simulates values for each grid cell independent of the rest of the surface. The error simulation process begins with the generation of random error grids that have the same dimensions as the original DEM. The values in the error grids are generated from a normal probability distribution that has a standard deviation (stdv) equal to the rms error of original DEM and a mean of zero. The error grids are then added to the original DEM, creating a realisation of the original DEM. The Monte Carlo method of error propagation is advantageous for this study because it is not affected by the exact formulation of the algorithms. Instead, it simply treats the algorithms as ‘black boxes’ whose response to the error in input is studied from the

resulting outputs. Therefore, with Monte Carlo simulation it is possible to determine the actual effect of error, regardless of global operations or a selection process in the algorithms. The primary disadvantage of the Monte Carlo method is the computational requirements. In order to produce statistically reliable results, a large number of realisations should be run. Obviously, a difficulty with the approach is determining exactly how many times an algorithm should be implemented in order to minimise the computational load, while still insuring reliable results. A few studies have been done where only a few runs have been used (20e50) (Fisher, 1991; Holmes et al., 2000), but other research has shown that this is not enough to produce reliable results (Veregin, 1997). A potential problem with a Monte Carlo error simulation that we have described is that it does not preserve the spatial structure of the landscape elevation. Neighbouring elevations in a landscape are correlated. The degree of correlation is driven by the ruggedness of the terrain which is in turn driven by the geological forces that have shaped the surface of the earth. Different landscapes each with a different geological history will have different degrees of correlation between neighbouring elevations. One might describe this correlation with geostatistical tools, such as the variogram, and attempt to preserve the spatial structure of the elevation by preserving the variogram. However, the variogram estimated from a DEM does not necessarily reflect that of the landscape. There have been no studies that we know of that compare the geostatistics of a carefully measured set of elevations, say by conventional surveying or GPS, with that of a DEM product. So, while we could estimate the variogram from our DEM and preserve it in each realisation, this would, in our opinion, be an arbitrary exercise. We recognise the

720

L.D. Raaflaub, M.J. Collins / Environmental Modelling & Software 21 (2006) 710e732

need for measuring and modelling the spatial structure of elevations, but we did not wish to implement this in an arbitrary fashion. 3.2.1.2. Generation of error realisations. In the literature, choosing the number of realisations used in Monte Carlo error simulations often seems to be fairly arbitrary. In order to avoid a similar arbitrary selection in this project, a more substantial method was required. The method used was based on the simple idea that realisations are only required if they significantly effect results. For example, if using 500 realisations produces the same result as using 200 realisations, then there is no need to use the additional 300 realisations. To apply this method, it is first necessary to generate realisations using the same error distribution. As the number of error realisations is increased, the standard deviation of each grid cell over all the realisations is determined, and a grid of standard deviations is generated. The standard deviation of this grid of standard deviations is found after the addition of each realisation, and this value is compared with the previous value. If the change in the standard deviation is not significant, then the addition of the realisation does not introduce any new information. Therefore, any realisations after this point would be redundant. An example of this procedure is shown below in Fig. 4 which shows the standard deviations versus the number of realisations for error realisations generated using a standard deviation of 3 m. After 500 realisations the curve begins to level out, and adding additional realisations has little impact on the resulting standard deviation. Consequently, a Monte Carlo simulation using 500 realisations would be enough to produce statistically reliable results. It should be noted that as long as the standard deviations of the error realisations are the same, the size

stdv of the stdv (m)

0.8 0.6 0.4 0.2 0.0

0

200

400

600

800

1000

Number of Realisations Fig. 4. The standard deviation of the standard deviation error matrix for uncorrelated random error with a standard deviation of 3 m. The standard deviation of each grid cell is taken over the number of realisations specified. The standard deviation is then taken of the standard deviation matrix.

of the grid used does not matter since the resulting shape of the curve does not change. The only change that can be seen is a slight vertical displacement, where the smaller the grid the higher the standard deviations. 3.2.1.3. Correlated error. Up to this point it has been assumed that DEM error is uncorrelated from one elevation to the next. This is a worst case scenario since the errors in neighbouring elevations tend to be positively correlated. This error correlation is related to terrain characteristics and the photogrammetric and/or radargrammetric methods used to generate the DEM. In general, the rougher the terrain the more difficult it is to measure elevation values. Elevation values located in close proximity will have similar terrain characteristics, and hence exhibit similar error properties. For DEM products derived from inSAR, the errors are also related to coherence which is driven by stability of the dielectric geometry between passes. One should not forget that many DEM products in use today are derived from topographic maps. The error characteristics of these data are largely unknown. While it is likely that DEM errors are correlated, information on the structure of this correlation is not supplied with the accuracy information provided with DEM products. The DEM used in this study was no exception. To investigate the effect of correlated error in this study, we chose what might be called a first order model of spatial correlation, a Gaussian correlation function. We constructed a 21 ! 21 kernel matrix consisting of values from a folded normal distribution and convolved it with the DEM error grid. The folded normal, with a standard deviation of 5 m and a mean of 0, was centred in the matrix so that the centre value had the largest weight, with cells farther away from having smaller weights. This agrees with the assumption that the closer the points, the more correlated the errors. The normal distribution is considered folded because it uses the distance from the centre point (r2 Z x2 C y2) as opposed to the x and y components individually. A standard deviation of 5 m was chosen for the folded normal since it created a distribution that fit the 21 ! 21 kernel matrix well. However, without any knowledge of the correlation, any standard deviation value could have been used. Convolution with a folded normal kernel successfully correlates the error; however, it has the undesirable side-effect of reducing the standard deviation of the entire matrix. Since the stated accuracy of the DEM was G3 m, the error realisations of the DEM need to have a standard deviation of 3 m, regardless of its correlation. In order to maintain a standard deviation of 3 m after correlation, the initial uncorrelated matrices had to have a much higher standard deviation. The initial standard deviation required in order

L.D. Raaflaub, M.J. Collins / Environmental Modelling & Software 21 (2006) 710e732

to produce a resulting standard deviation of 3 m can be determined by sfinal sinitial Zpffiffiffiffiffiffiffiffiffiffiffiffiffi P 2; W

ð8Þ

where s is the standard deviation, and W are the weights of convolution matrix. The weights P of the convolution matrix are normalised to one ð WZ1Þ: Based on Eq. (8), the initial standard deviation of the uncorrelated error grids used so that the standard deviation of the error grids after correlation was 3 was determined to be 49.52. In order to ensure that the convolution with a 21 ! 21 folded normal did increase the correlation of an error grid, the autocorrelation function of both the uncorrelated and correlated error grids were calculated. The autocorrelation function of the uncorrelated error grid, as expected, only produced a spike at the centre of the distribution e confirming that the uncorrelated error grid is completely uncorrelated. The autocorrelation function of the correlated error grid, on the other hand, had a distinctive 2 dimensional bell shape curve. This indicated that the correlated error is indeed correlated with distance. For comparison, a horizontal cross section of the two autocorrelation functions is provided in Fig. 5. The method used to apply the Monte Carlo method of error propagation to correlated error is the same as that used for the uncorrelated error. The only difference is that the error grids need to be correlated before they are added to the original DEM. It should be noted that correlating the error does not effect the number of realisations required. A plot of the standard deviation of the standard deviation for correlated error is exactly the same as that shown in Fig. 4. Therefore, as with the uncorrelated error, 500 realisations of the correlated error were used to examine the effect of correlated error on the calculation of topographic parameters.

Autocorrelation

300000

Uncorrelated Correlated

250000 200000 150000 100000 50000 0 −50000 −3500 −2625 −1750 −875

0

875

1750

Distance from the Centre (m)

2625

3500

Fig. 5. Horizontal cross section of the autocorrelation function (the correlation function is radially symmetric). Distances are based on a 25-m grid spacing.

721

3.2.2. Pit removal In order to determine the topographic index for each grid cell in the watershed, both the slope and upslope area must be defined everywhere. In addition, flow must be routed out of the watershed and not accumulate in localised areas. Problems occur when the slope e and by extension the aspect e are undefined. This occurs in flat areas and, depending on the slope and aspect algorithm used, at closed depressions (pits). An undefined aspect becomes a problem when trying to calculate upslope area because there is no knowledge on how to route water out of a cell. Therefore, to get accurate values of upslope area, and to be able to create a topographic index grid that can be applied to a hydrological model, areas of undefined slope and aspect must be corrected. In the watershed used in this study, there were only a few single cell pits. Therefore, only the problem of single localised pits was dealt with, and not the more complex problems associated with larger pits and flat areas. Because of the design of the algorithms studied, the only algorithms that had problems with these pits were those based on the steepest neighbour (MDG, MAG and MDN). The problems of pit removal in these algorithms could be further restricted to the problem of finding a method to route water through a pit using the MDG method, since the other two algorithms can be ultimately reduced to a form of this algorithm. A simple way to handle pits is to manually remove them from the DEM before the calculation of topographic parameters. In this study this was not feasible because of the many error realisations. Adding error grids to the original DEMs to create error realisations had the effect of generating pits that were not present in the original DEM. Because of this, the pits had to be removed after the creation of the aspect grids for each error realisation and could not just be removed from the original DEM. The method employed to deal with pits essentially involved rerouting flow out of the pit. This required the initial creation of an aspect grid. From the aspect grid, those cells with undefined aspect values were selected. An example of a pit cell, depicted in red, is shown in the left hand figure in Fig. 6. Once a pit was found it was assigned an aspect that corresponded to the direction of the neighbour with the lowest elevation. The flow path from this neighbour was then traced to ensure that it flowed out of the watershed. If the flow path did not leave the watershed, the aspect of the neighbour was changed to the direction of its next lowest neighbour. This process was repeated until a flow path was found that flowed out of the watershed. Reassignment of aspect was restricted to cells that were not already being rerouted. Using this method, only a few aspect values typically required modification. An example of a flow path reroute is shown in the right hand figure in Fig. 6.

722

L.D. Raaflaub, M.J. Collins / Environmental Modelling & Software 21 (2006) 710e732

To estimate overall variability, we calculated the mean standard deviation of slope/aspect and the standard deviation of standard deviation for the whole watershed. We are using the standard deviation of slope and aspect as the basis for evaluating the sensitivity of the algorithms to elevation uncertainty in the DEM. Large standard deviation indicates high sensitivity and small standard deviation indicates low sensitivity. The standard deviation maps will show any spatial pattern in error sensitivity while the watershed average will indicate the overall average error sensitivity. In addition we calculate the standard deviation of the error sensitivity over the whole watershed. This statistic will summarise the overall variability of the error sensitivity. It is desirable for slope/ aspect algorithms to be robust against uncertainty, hence high quality algorithms will have low sensitivity to elevation uncertainty and small standard deviations. In addition we would like an algorithm that minimises the variability in error sensitivity. Hence, we want both the mean and the variability about the mean to be small. 3.3. Upslope area

Fig. 6. Flow directions of cells using MDG before pit correction (top) and after pit correction with the new flow path marked in light blue (bottom). Quiver lines represent the aspect direction and the red dot is the location of the pit. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Once the problem of an undefined aspect for pits has been addressed, the only problem remaining is that of an undefined slope. The slope for a flat area is, by definition, zero. Pits can also be viewed as flat because before there can be any outflow from them, their water level must theoretically reach a height level with that of their lowest neighbour. However, in either case the zero slope value would create an undefined topographic index. Therefore, to avoid any undefined topographic indices, all flat grid cells, and by extension all pits in the watershed were assigned a very small value. This very small value was arbitrarily selected to be 0.000001. 3.2.3. Calculation of topographic parameters We generated 500 DEM realisations each for uncorrelated and correlated errors. For each DEM we calculated slope and aspect using the algorithms outlined in Section 2.3. For each set of 500 slope and aspect maps we calculated two maps, a mean slope/aspect map and a standard deviation slope/aspect map. In addition, we used the standard deviation maps to estimate overall variability of slope/aspect error for the whole watershed.

Aspect values of the watershed are used to calculate the upslope area for each grid cell in the watershed. The determination of upslope area essentially involves in tracing the flow paths of water through a watershed to the drainage outlet. The algorithms available to calculate the upslope area compliment the algorithms used to calculate the slope and aspect. An upslope area grid was determined for each realisation and for each of the eight slope/aspect algorithms. Since the ultimate goal was to find the topographic index value for each cell and since the topographic index uses the upslope area in the form of the specific catchment area, the results for upslope area were expressed in terms of the specific catchment area (units are m instead of m2). The only difference between the specific catchment area and the upslope area is that the specific catchment area is the upslope area divided by the grid size (25 m in this study). 3.4. The topographic index Topographic index grids were calculated for each of the eight slope and aspect algorithms and for every error realisation of the watershed. The value of the topographic index for each cell in a realisation was found using the slope and corresponding upslope area grids. The topographic index for each grid cell was then found using Eq. (6). 3.4.1. Edge cells Each of the slope and aspect algorithms require knowledge of a cell’s neighbours. This means that all of these algorithms experience problems at the edges of

723

L.D. Raaflaub, M.J. Collins / Environmental Modelling & Software 21 (2006) 710e732

a DEM. Rather than invent new algorithms to deal with edge cells we ignored the edge cells which had the effect of ‘shrinking’ the watershed by one grid cell along its edge. Since all the information in the watershed was required for hydrological modelling, the watershed DEM used to determine the slope and aspect was expanded by one grid cell all around. By doing this, none of the information available for the watershed itself was lost.

Table 1 Mean ( s) and standard deviation (ss) of the standard deviation of slope for the entire watershed Error

Algorithm

 s

ss

Uncorrelated

MDG MAG MDN TPP FCN FDN ENU OODS

0.1325 0.1013 0.0785 0.1495 0.0817 0.0588 0.0483 0.0512

0.0144 0.0110 0.0105 0.0250 0.0053 0.0030 0.0022 0.0024

Correlated

MDG MAG MDN TPP FCN FDN ENU OODS

0.0180 0.0179 0.0164 0.0182 0.0176 0.0174 0.0175 0.0175

0.0008 0.0008 0.0086 0.0006 0.0006 0.0006 0.0006 0.0006

4. Experimental results We found that the MDG and ENU algorithms represent two extremes in their results with MDG showing the greatest variability and the ENU showing the least variability. This is consistent with our expectations since MDG is based on a selection process that uses only one neighbour to calculate slope and aspect, and is the method most often used in hydrology. ENU, in contrast, uses all eight neighbours in the calculation of slope and aspect. This method has also been shown to produce slope and aspect values that more closely approximate those of the terrain (Skidmore, 1989; Guth, 1995). Hence we present mean and standard deviation maps for these two algorithms only. In addition we present mean and standard deviation maps for uncorrelated error only. The filtering procedure we used to correlate the error naturally smoothes the error and reduces its effect. Therefore, the correlated error maps are similar to the uncorrelated error maps, and differ only in the magnitude of the error and in the extent of high error regions. The error map based on correlated error is, however, provided for topographic index because the calculation of it requires all of the other topographic parameters determined in the study. Consequently, it nicely summarises the cumulative effect of the error. 4.1. Slope The overall error sensitivity of the slope estimate for the eight algorithms is presented in Table 1. All eight algorithms are, as expected, much less sensitive to correlated error than they are to uncorrelated error. This is consistent with our expectations since all algorithms are based on local neighbourhood operations, and correlating the error reduces the variation of elevation in local neighbourhoods. In other words, the error in a local neighbourhood will be more of an elevation bias and will hence will have less of an impact on the resulting slope value. Another interesting observation regarding correlated error is that correlating the error reduces the variation between algorithms. The mean standard deviation for uncorrelated error ranges between 0.1495 and 0.0483, while for correlated error the mean standard deviation only ranges between 0.0182 and 0.0164. Therefore, all slope algorithms are equally

sensitive to correlated error. Again, the reason is that the correlated error results in an elevation bias. Because none of the slope algorithms use the absolute elevation values, this bias does not effect their calculation. While the error sensitivity in slope for the eight algorithms are fairly uniform for correlated error, the differences between algorithms are much larger for uncorrelated error. A general trend is that the sensitivity to error decreases as the number of neighbours used in the algorithm increases. The overall variability (i.e. sensitivity to error) in the steepest neighbour algorithms is, in general, larger than for the four neighbour algorithms, which is in turn larger than those of the eight neighbour algorithms. This relationship is consistent with our expectations since the more neighbours involved in the algorithm, the less likely the slope values will be affected by one, or a few, highly erroneous elevation values in neighbouring cells. Since the MDN algorithm calculates slope based on an average of all downslope neighbours, it is also reasonable that it is less sensitive to error than the other steepest neighbour algorithms and even the FCN algorithm. An apparent variation to this trend is the TPP method, which had the largest error sensitivity of all the algorithms. However, this algorithm uses the fewest neighbours and is most sensitive to the orientation of the grid. Even though the steepest neighbour algorithms use only one (or more than one as is the case in MDN) neighbour to calculate slope, they do take into account the other eight neighbours in the selection process. The TPP method, on the other hand, uses only two specific neighbours, ignoring the other six. This limitation can make this method more vulnerable to highly erroneous elevation values. Based on the significant differences in slope error between the FCN and FDN methods, along with the differences in slope error between the ENU and OODS

724

L.D. Raaflaub, M.J. Collins / Environmental Modelling & Software 21 (2006) 710e732

methods, it appears that this watershed suffers from a grid orientation bias. That is, over the entire watershed, slopes are more dominant in the diagonal directions than they are in the cardinal directions, and consequently, the diagonal direction is less sensitive to error than the cardinal direction. This is evident by the error sensitivity in the FDN slope algorithm being smaller than that in the FCN slope algorithm, and also in the error sensitivity in the ENU algorithm being smaller than that in the OODS algorithm. In the ENU algorithm all eight neighbours are weighted the same, where as in the OODS algorithm the diagonal neighbours are weighted less than the cardinal neighbours. This essentially results in, for this particular watershed, the more erroneous neighbours having a more significant role in the determination of slope, resulting in a higher error sensitivity. The significance of the differences between the overall error in the FCN and FDN slope algorithms, as well as between the ENU and OODS slope algorithms, was statistically confirmed using a student’s t-test. While there has been little to no research into the error sensitivity of slope algorithms, there has been some researches into the ability of these algorithms to accurately estimate the ‘true’ values of slope (Guth,

1995; Jones, 1998; Skidmore, 1989; Veregin, 1997; Quinn et al., 1991). These studies reported similar results to those found in this study e the greater the number of neighbours used in a slope algorithm the more accurately the slope values agreed with the actual value of slope. The steepest neighbour algorithms tended to produce poorer results, while the ENU and FCN algorithms tended to produce more accurate results. The spatial structure of slope error sensitivity in the watershed may be found in the mean and standard deviation maps for the MDG and the ENU methods as shown in Fig. 7. These maps allow to estimate the effect of slope value on algorithm sensitivity to error and to estimate bias in the slope values. From the mean slope maps it is apparent that the MDG method, on average, produces higher slope values than the ENU method. The mean value over the entire grid for the former was 0.446, while for the latter it was 0.406. This is a significant difference (confirmed using a student’s t-test) and similar to results obtained by Guth (1995), who found that the MDG method, and by extension the MAG method, tended to give significantly steeper slopes, especially in steep terrain. Though the mean slope values at each cell are different between the two methods, the overall slope pattern is similar.

Fig. 7. Slope error maps for uncorrelated error. The mean and standard deviation are determined for each cell over the 500 realisations. Mean slope using MDG (top left). Standard deviation of slope using MDG (top right). Mean slope using ENU (bottom left). Standard deviation of slope using ENU (bottom right).

725

L.D. Raaflaub, M.J. Collins / Environmental Modelling & Software 21 (2006) 710e732

Although the spatial pattern in the mean slope maps for the ENU and MDG methods is similar, there are significant differences in the standard deviation maps. The standard deviation map for the ENU method looks like random noise, hence the variability in error sensitivity in slope for the ENU method is generally uniform across the watershed. An exception is a small region near the upper right hand edge of the watershed where the standard deviation are quite small indicating a low sensitivity to error. This region is fairly flat and has low slope, which suggests that areas of low slope are less sensitive to error. This confirms the results presented by Veregin (1997). For the MDG method the error sensitivity appears to be correlated with slope value, however, this correlation is very slight. As with the ENU method, areas of low slope are less sensitive to error. However, in general, error sensitivity for the MDG method is much higher than that for the ENU method across the entire watershed. These results are confirmed by Fig. 8, which is a plot of the standard deviation of slope as a function of the mean slope. For both methods there appears to be a threshold mean slope below which the standard deviation (i.e. error sensitivity) drops off. It is clear from this figure that the overall error sensitivity is much higher for the MDG method and variability of this error sensitivity is also very high in comparison to the ENU method. 4.2. Aspect The overall error sensitivity of the aspect estimate for the eight algorithms is presented in Table 2. The results are similar to those for slope in several respects. First, the sensitivity to correlated error is much smaller than that for uncorrelated error. For slope the results for correlated error were about an order of magnitude smaller than that for uncorrelated error; for aspect the correlated error is around 20% of the uncorrelated error. The error sensitivity of the algorithm groups is also similar: the eight neighbour algorithms have the lowest error followed by the four neighbour algorithms followed by the steepest neighbour and finally TPP.

Fig. 8. The standard deviation of slope as a function of the mean slope for each cell in the watershed for uncorrelated DEM error.

Table 2 Mean ( s) and standard deviation (ss) of the standard deviation of aspect in degrees for the entire watershed Error

Algorithm

 s

ss

Uncorrelated

MDG MAG MDN TPP FCN FDN ENU OODS

26.64 26.90 46.89 34.36 16.58 11.53 9.26 9.81

15.29 16.08 37.71 31.69 17.94 13.74 11.74 11.78

Correlated

MDG MAG MDN TPP FCN FDN ENU OODS

5.27 5.61 10.84 3.35 3.18 3.32 3.20 3.18

9.48 10.19 21.92 4.94 4.45 5.59 5.24 4.32

The difference is the performance of the MDN algorithm which had the greatest sensitivity to error: about double that for the other steepest neighbour algorithms and 30% higher than the TPP result. This came as a surprise as the slope sensitivity for the MDN algorithm was the lowest of the steepest neighbour group. The variability of the error (ss) mirrors the result for the error itself: the higher the error, the greater the variability. Maps of the mean and standard deviation of the aspect estimates across the watershed for the MDG and ENU methods are presented in Fig. 9. From the mean aspect maps it can be seen that there is little difference in the aspect values between the two algorithms. The only difference between them is that the ENU values are slightly smoother. This is consistent with our expectations as the MDG method is restricted to aspect values that are in 45  intervals, while the aspect values for the ENU method can take any value, resulting in smoother variation across the landscape. The standard deviation maps for both algorithms (Fig. 9) are also similar. This is not unexpected, because both algorithms have similar overall error values. The error sensitivity is fairly constant across the watershed for both algorithms. Error tends to be a little higher for MDG than for ENU in the lower slope areas and along the drainage network, but this variation is not as large as that seen in slope errors. For both algorithms the area of highest aspect error is located in the same location as the area of lowest error for slope, i.e. the upper right hand edge. As shown in Fig. 7, this is a low slope area. It is understandable that areas of low slope are also areas of high aspect error since a small error in slope can dramatically effect the aspect. Fig. 10 (top) shows the standard deviation of aspect as a function of mean aspect and mean slope. There is no apparent relationship between aspect error sensitivity

726

L.D. Raaflaub, M.J. Collins / Environmental Modelling & Software 21 (2006) 710e732

Fig. 9. Aspect error maps for uncorrelated error. The mean and standard deviation are determined for each cell over the 500 realisations. Aspect values are expressed in degrees. Mean aspect using MDG (top left). Absolute standard deviation of aspect using MDG (top right). Mean aspect using ENU (bottom left). Absolute standard deviation of aspect using ENU (bottom right).

and the value of the aspect estimate. Fig. 10 (bottom) shows a strong relationship between aspect error sensitivity and slope value. As the slope increases, the aspect error sensitivity also decreases. The aspect error sensitivity is quite high at low slopes (flat areas) and decreases quite rapidly to a threshold after which it decreases more gradually. This relationship is much stronger for the ENU method where the variability of error sensitivity against slope is very small. For the MDG method the variability of aspect error sensitivity against slope is high as was the case for the variability of slope error sensitivity. Hence, as with slope, the error sensitivity for aspect grows as the number of neighbours participating in the estimate falls. 4.3. Upslope area The overall error sensitivity in upslope area is given in Table 3. As was found for slope errors, all eight algorithms are less sensitive, in general, to correlated error than they are to uncorrelated error. The exceptions are the MDN and TPP algorithms where the sensitivity to correlated error is significantly higher. The reasons for this are similar to those given for slope. However, unlike the errors in slope, correlating the DEM error does not reduce the variation in error between the

algorithms. This difference from the slope results can be attributed to the difference between global and neighbourhood operations. Slope, as calculated in this study, is restricted to neighbourhood operations, and, hence, treats correlated error as a bias. Upslope area, on the other hand, is more of a global operation, and hence is affected by values outside the region of correlation. The overall error between the eight slope/aspect algorithms, for both uncorrelated and correlated DEM error, appears to follow two patterns. These two patterns correspond to the two distinctive types of algorithms: those based on the steepest neighbour, and those that are not. For algorithms in the steepest neighbour group, it is apparent that the smaller the dependence on a single steepest neighbour, the smaller the sensitivity to DEM error. The MDG algorithm, which is based entirely on a single steepest neighbour has the larger error, while the MDN algorithm, which uses all downslope neighbours, has the lowest error. Based on the way the upslope area is calculated for these two algorithms, this result is understandable. The MDG algorithm only allows single cell flow e flow from one cell flow completely into a single neighbouring cell. Since the flow is never divided, it is difficult for this method to ‘recover’ from erroneous elevation values. If the flow is misdirected for a specific cell in a given realisation by even one grid cell, there is little to no

727

L.D. Raaflaub, M.J. Collins / Environmental Modelling & Software 21 (2006) 710e732

Table 3 Mean and standard deviation of the standard deviation of upslope area for the entire watershed

Fig. 10. The standard deviation of aspect as a function of the mean aspect (top) and mean slope (bottom) for each cell in the watershed for uncorrelated DEM error.

chance that the flow will follow the same path as that in another realisation. Since upslope area of one grid cell is connected to the upslope area of the other cells in the watershed, errors compound downslope. For the MDN algorithm flow is divided between all downslope neighbours. It is therefore possible for flow to diverge and converge as it travels downslope. Because of this, the variation in upslope area for each grid cell will be less than that observed using the MDG algorithm. The reason that the MAG algorithm is less sensitive to error than the MDG algorithm, even though it is also based on a single neighbour, is that it finds the neighbour with the maximum slope, not just the neighbour with the maximum downslope. Therefore, this method has a greater ability to choose the ‘correct’ flow direction. However, the MAG method is still limited by the same factors affecting the MDG method, and therefore it has a sensitivity to error that is larger than that of the MDN method. The five remaining slope/aspect algorithms (TPP, FCN, FDN, ENU and OODS) exhibit an error pattern that is the opposite of that seen for the steepest neighbour based algorithms and of the slope/aspect

Error

Algorithm

Mean standard deviation (m)

Standard deviation of standard deviation (m)

Uncorrelated

MDG MAG MDN TPP FCN FDN ENU OODS

4314.49 797.12 151.79 87.22 180.34 226.99 253.50 264.83

18212.98 1613.58 464.93 147.99 526.34 732.47 741.65 793.60

Correlated

MDG MAG MDN TPP FCN FDN ENU OODS

1099.26 615.43 281.52 98.85 133.28 175.00 151.80 143.80

8465.67 2189.31 1391.18 369.43 689.04 774.27 657.34 667.55

results in general. For these algorithms, the error in upslope area increases as the number of neighbours used in the algorithm increases. All of these methods calculate upslope area in a similar manner using plane fitting. Therefore, the only difference between the algorithms is the number of neighbours used to determine aspect. None of these methods placed any restrictions on aspect values, i.e. aspect was not limited to 45  bins. Since the error in aspect decreased as the number of neighbours increased, the results obtained for upslope area are a little baffling. In many cases, the goal of estimating upslope area is the extraction of a drainage network. Comparing the two mean upslope area plots (Fig. 11) it is clear that the MDG method produces a visible drainage network with a definable outflow location, while the ENU method, except for a few disconnected lines, does not. Similarly poor connectivity is observed in all algorithms that are not based on the steepest neighbour principle. Furthermore, as the number of neighbours used in an algorithm decreases, the more disconnected the drainage network becomes. Or, alternatively, the fewer neighbours used in an algorithm, the more uniform the upslope area values are across the watershed. If the upslope area is relatively uniform across the watershed, then it is less likely that there will be any variation in a single cell’s value between realisations. If, on the other hand, there are distinctive features that have widely different values, there is a greater opportunity for variation over realisations. It is therefore understandable that, for slope and aspect algorithms not based on the steepest neighbour, the overall error in upslope area will increase as the number of neighbours increases. In other words, the fewer the neighbours, the more uniform the upslope area, and the more uniform the upslope area, the smaller the variation

728

L.D. Raaflaub, M.J. Collins / Environmental Modelling & Software 21 (2006) 710e732

Fig. 11. Upslope area error maps for uncorrelated error. Shown are the mean (top left) and standard deviation (top right) of the upslope area using MDG and mean (bottom left) and standard deviation (bottom right) of the upslope area using ENU. The mean and standard deviation are determined for each cell over the 500 realisations. Upslope area values are expressed in metres.

outlet. Upslope area errors are generally larger for the MDG method as opposed to the ENU method whose errors are not only small but apparently unrelated to the upslope area value. An unusual phenomenon can be observed in Fig. 12 for the upslope area calculated using the MDG method. Three distinct curves are visible, each with a different maximum error value. These different curves likely correspond to data from different geographic locations

Standard Deviation

between realisations. For those algorithms not based upon the steepest neighbour principle, this explains the observed rise in upslope area error as the number of neighbours increased. Hence we have results that are somewhat conflicting. The overall error sensitivity of the steepest neighbour algorithms (as a group) is poorer than the fitted plane algorithms, and yet the latter group is unable to produce a drainage network. In fact, the best performer in terms of error sensitivity, the TPP algorithm, is the poorest performer in terms of drainage network production. However, notwithstanding the performance of the TPP algorithm, the best performance for both error sensitivity and drainage network production was the MDN algorithm. One would expect that the error sensitivity of upslope area would increase with the value of upslope area (as is visible in Fig. 11) where the drainage network is the location of high upslope area values. This relationship is more apparent in Fig. 12 where the standard deviation of upslope area for each cell in the watershed is plotted as a function of the mean. The error for the MDG method does have an upward trajectory against the upslope area value as expected, although it increases and then decreases. The largest error is at the drainage

2

x 10

5

1.5 1 MDG ENU

0.5 0

0

0.5

1

1.5

2

Mean Upslope Area (m)

2.5

3

x 10

5

Fig. 12. The standard deviation of upslope area as a function of the mean upslope area for each cell in the watershed for uncorrelated DEM error.

L.D. Raaflaub, M.J. Collins / Environmental Modelling & Software 21 (2006) 710e732

in the watershed, i.e. along the different branches of the drainage network. As observed before, flatter areas tend to be more sensitive to error. It stands to reason therefore, that drainage networks from flatter areas will in general experience less upslope area error sensitivity than drainage networks from steeper areas. 4.4. Topographic index The overall affect of uncorrelated and correlated DEM error on the eight algorithms used to calculate the topographic index are presented in Table 4. The patterns observed for errors in the topographic index between algorithms and between uncorrelated and correlated error are a combination of those observed for both slope and upslope area. Correlated DEM error reduces both the amount and variability of error in the topographic index for all eight algorithms. As with slope, the topographic index error decreases with an increase in the number of neighbours utilised in the algorithm. However, unlike for slope, the steepest neighbour algorithms (MDG, MAG and MDN) have the greatest error sensitivity. This observation includes the TPP method, which (again, unlike slope) has a smaller error than all of the steepest neighbour algorithms. The reason for this difference is related to the large error in upslope area found in the steepest neighbour algorithms. Another similarity between slope and topographic index error is the apparent grid orientation bias. The explanations for the observed patterns in topographic index error are the same as those given for the errors in slope and upslope area. Topographic index error maps for the MDG and ENU methods are presented in Fig. 13. As with the overall error discussed above, the mean and standard Table 4 Mean and standard deviation of the standard deviation of topographic index for the entire watershed Error

Algorithm

Mean standard deviation

Standard deviation of standard deviation

Uncorrelated

MDG MAG MDN TPP FCN FDN ENU OODS

1.9297 1.2218 1.2392 0.7970 0.7022 0.5841 0.5783 0.6112

1.1937 0.4562 0.8833 0.2647 0.3602 0.3909 0.4060 0.4129

MDG MAG MDN TPP FCN FDN ENU OODS

0.4559 0.4774 0.3975 0.2883 0.2551 0.2330 0.2342 0.2360

0.4658 0.4638 0.4658 0.3789 0.3529 0.3146 0.3294 0.3395

Correlated

729

deviation maps have attributes similar to those seen in the error maps for slope and upslope area. By comparing the two mean topographic index maps it is apparent that the MDG method, on average, produces higher values than the ENU method. While a drainage network is visible in the mean plots for both methods, the one produced using the MDG method is more highly defined. These observations were also made for the mean maps for slope and upslope area. The standard deviation maps for the MDG and ENU algorithms have appearances that have already been described by the error maps of slope and upslope area. For both algorithms, the drainage network is visible in the error maps, but is more prominent for the MDG algorithm. Error for both methods is greatest along the drainage network and in areas of low slope. Topographic index error is, in general, larger for the MDG method than the ENU method across the entire watershed. To compare the error maps generated using uncorrelated and correlated DEM error, the topographic index error maps resulting from correlated error were created for the MDG and ENU algorithms, and are shown in Fig. 14. For both algorithms there is very little difference visible between the mean plots for the uncorrelated error and the mean plots for the correlated error. This is reflected in the small average differences in the mean topographic index values between the two types of error simulations. For the MDG algorithm the difference was 0.63, and for the ENU algorithm it was an even smaller 0.33. The difference between error types is, however, much more significant when the standard deviations are examined. This is because correlated error significantly reduces the error sensitivity of both algorithms. While features visible in the topographic index error plots generated using uncorrelated error are visible in the error plots generated using correlated error, correlated error does tend to make the error more uniform across the watershed. Therefore, while correlated DEM error does not effect the value of topographic index, it does reduce the error sensitivity.

5. Conclusions This research studied the effect of uncorrelated and correlated DEM error on the determination of slope and aspect in particular to determine which topographic algorithms had the greatest sensitivity to DEM error. Conclusions that can be made from the results presented in the previous section are as follows. 5.1. Correlated versus uncorrelated DEM error Both slope and aspect error sensitivity was significantly smaller in the presence of correlated DEM error than for uncorrelated DEM error. All algorithms

730

L.D. Raaflaub, M.J. Collins / Environmental Modelling & Software 21 (2006) 710e732

Fig. 13. Topographic index error maps for uncorrelated error. Shown are the mean (top left) and standard deviation (top right) using MDG and mean (bottom left) and standard deviation (bottom right) using ENU. The mean and standard deviation are determined for each cell over the 500 realisations.

produced slope estimates that were equally sensitive to correlated error, while the aspect error sensitivity mirrored that for uncorrelated error.

the steepest neighbour principle are, in general, more dependent on the slope value. 5.3. Aspect

5.2. Slope The sensitivity of the slope estimate to DEM error decreased as the number of neighbours used in the algorithm increased. This is true for both the steepest neighbour algorithms and those algorithms that use a fixed number of neighbours. Furthermore, the steepest neighbour algorithms were more sensitive to DEM error than the fixed neighbour methods. When the eight algorithms are ranked from highest to lowest error sensitivity the result is as follows: TPP, MDG, MAG, FCN, MDN, FDN, OODS and ENU. For the two algorithms studied in detail (ENU and MDG), the sensitivity of the slope error was smallest in flat regions where slope was small. Outside these regions the ENU algorithm had an error sensitivity that was fairly uniform. In contrast, the error sensitivity of the MDG method was dependent upon the value of slope. We believe that these results are representative of the two different categories of slope algorithms and that, compared to those algorithms that use a fixed number of neighbours, the error sensitivity of algorithms based on

The sensitivity of the aspect estimate to DEM error also decreased as the number of neighbours used in the algorithm increased. The variation in the magnitude of aspect error across the watershed is independent of the aspect value itself. Rather, it is related to the value of slope. In particular, areas of low slope tend to correspond to areas of high aspect error. Overall the ENU algorithm had the lowest sensitivity to error. There were definite performance advantages in using eight neighbours over four neighbours and using a eight or four neighbour algorithm over a steepest neighbour method. The TPP method had the poorest performance overall, although the MDN algorithm had poorest performance for aspect. The hydrological reasoning behind the steepest neighbour algorithms does lead to algorithms whose estimates of slope and aspect are robust against the uncertainties in elevation that are commonplace in digital elevation models. Users requiring robust estimates of slope and aspect are advised to use eight neighbour algorithms.

L.D. Raaflaub, M.J. Collins / Environmental Modelling & Software 21 (2006) 710e732

731

Fig. 14. Topographic index error maps for correlated error. Shown are the mean (top left) and standard deviation (top right) using MDG and mean (bottom left) and standard deviation (bottom right) using ENU. The mean and standard deviation are determined for each cell over the 500 realisations.

5.4. Upslope area For the upslope area derived from slope and aspect algorithms based on the steepest neighbour principle, it was observed that the smaller the dependance on a single downhill neighbour, the smaller the error sensitivity. Conversely, for those algorithms not based on the steepest neighbour, the error sensitivity of upslope area increases with an increase in the number of neighbours used in the algorithm. The error sensitivity of the upslope area is greatest for those algorithms that produce well-defined drainage networks. This is because it is along the drainage networks that the largest upslope area errors occur. The steepest neighbour algorithms are the most adept at producing well-defined networks, and consequently have the greatest error sensitivity. The eight algorithms can be ranked from highest to lowest amount of upslope area error as a result of DEM error as follows: MDG, MAG, MDN, OODS, ENU, FDN, FCN and TPP. 5.5. Topographic index Errors in the topographic index are largest for the algorithms based on the steepest neighbour principle. Within this category, a similar error pattern is observed as was seen for both the slope and upslope area: the

greater the dependence on a single neighbour the greater the error sensitivity. The topographic index error pattern for the algorithms that use a fixed number of neighbours is, again, similar to that observed for slope: the more neighbours used in an algorithm, the less sensitive the result is to DEM error. Errors in topographic index are greatest along the drainage network and in areas of low slope. As discussed previously, however, the steepest neighbour algorithms produce the most clearly defined drainage networks, and consequently have a greater overall error. Thus, when the eight algorithms can be ranked from largest to smallest amounts of error, the steepest neighbour algorithms come out on top: MDG, MDN, MAG, TPP, FCN, OODS, FDN and ENU. Some final caveats on the use of our results and future directions. We have studied a single watershed in the foothills of eastern slopes of the North American Rocky Mountains. As we have stated, the spatial structure of both the elevations and their errors are influenced by the particular landscape under study, so similar analyses of other landscapes would be useful. In this study we have not attempted to preserve the spatial structure of the elevations in each Monte Carlo realisation, nor have we attempted to impose any spatial structure on the errors. The reason, as we have stated, is that we have no knowledge of these spatial structures and we wish to limit the number of assumptions used in

732

L.D. Raaflaub, M.J. Collins / Environmental Modelling & Software 21 (2006) 710e732

the research. Building this knowledge base, while labour intensive as it implies extensive field work, would greatly benefit analyses of digital elevation data. Acknowledgements This research was funded largely through a GEOIDE NCE research grant to Collins. Additional funding came from an NSERC research grant to Collins. We further acknowledge the assistance of Dr. Caterina Valeo of the Department of Geomatics Engineering at The University of Calgary for many useful suggestions and the anonymous referee for helping to clarify several important issues. References Ackermann, F., 1992. High-quality digital terrain models e the SCOP program and derived products. EARSeL Advances in Remote Sensing 1 (3), 122e128. Band, L.E., 1986. Topographic partition of watersheds with digital elevation models. Water Resources Research 22 (1), 15e24. Beven, K.J., Kirkby, M.J., 1979. A physically based, variable contributing area model of basin hydrology. Hydrological Sciences Bulletin 24 (1), 43e69. Costa-Cabral, M.C., Burges, S.J., 1994. Digital elevation model networks (DEMON): a model of flow over hillslopes for computation of contributing and dispersal areas. Water Resources Research 30 (6), 1681e1692. Dozier, J., Strahler, A.H., 1983. Ground investigations in support of remote sensing. In: Colwell, R.N., Simonett, D.S. (Eds.), second ed. Manual of Remote Sensing, vol. 1. American Society of Photogrammetry, Virginia, pp. 959e986. Eyton, J.R., 1991. Rate-of-change maps. Cartography and Geographic Information Systems 18 (2), 87e103. Fairfield, J., Leymarie, P., 1991. Drainage networks from grid digital elevation models. Water Resources Research 27 (5), 709e717. Fisher, P.F., 1991. First experiments in viewshed uncertainty: the accuracy of the viewshed area. Photogrammetric Engineering and Remote Sensing 57 (10), 1321e1327. Freeman, T.G., 1991. Calculating catchment area with divergent flow based on a regular grid. Computers and Geosciences 17 (3), 413e422. Gallant, J.C., Wilson, J.P., 2000. Primary topographic attributes. In: Wilson, J.P., Gallant, J.C. (Eds.), Terrain Analysis: Principles and Applications. John Wiley & Sons, Inc., New York, pp. 51e85. Goodrich, D.C., Woolhiser, D.A., Keefer, T.O., 1991. Kinematic routing using finite elements on a triangular irregular network. Water Resources Research 27 (6), 995e1003. Guth, P.L., 1995. Slope and aspect calculations on gridded digital elevation models: examples from a geomorphometric toolbox for personal computers. Zeitschrift fu¨r Geomorphologie N. F., Supplementband 101, 31e52. Heuvelink, G.B.M., 1998. Error Propagation in Environmental Modelling with GIS. Taylor & Francis, London. Hodgson, M.E., 1995. What cell size does the computed slope/aspect angle represent. Photogrammetric Engineering and Remote Sensing 61 (5), 513e517. Hoffer, R.M., Fleming, M.D., Bartolucci, L.A., Davis, R.F., Nelson, S.M., 1979. Digital processing of landsat MSS and topographic data to improve capabilities for computerized mapping of forest cover types. LARS Technical Report 011579. Laboratory for Applications of Remote Sensing, Purdue University, West Lafayette, Indiana.

Holmes, K.W., Chadwick, O.A., Kyriakidis, P.C., 2000. Error in a USGS 30-meter digital elevation model and its impact on terrain modeling. Journal of Hydrology 233, 154e173. Holmgren, P., 1994. Multiple flow direction algorithms for runoff modelling in grid based elevation models: an empirical evaluation. Hydrological Processes 8, 327e334. Horn, B.K.P., 1981. Hill shading and the reflectance map. Proceedings of the IEEE 69 (1), 14e47. Hutchinson, M.F., 1989. A new procedure for gridding elevation and stream line data with automatic removal of spurious pits. Journal of Hydrology 106, 211e232. Jenson, S.K., Dominque, J.O., 1988. Extracting topographic structure from digital elevation data for geographic information system analysis. Photogrammetric Engineering and Remote Sensing 54 (11), 1593e1600. Jones, K.H., 1998. A comparison of algorithms used to compute hill slope as a property of the DEM. Computers and Geosciences 24 (4), 315e323. Martz, L.W., de Jong, E., 1988. CATCH: a FORTRAN program for measuring catchment area from digital elevation models. Computers and Geosciences 14 (5), 627e640. Moore, I.D., O’Loughlin, E.M., Burch, G.J., 1988. A contour-based topographic model for hydrological and ecological applications. Earth Surface Processes and Landforms 13, 305e320. O’Callaghan, J.F., Mark, D.M., 1984. The extraction of drainage networks form digital elevation data. Computer Vision, Graphics, and Image Processing 28, 323e344. O’Neill, M.P., Mark, D.M., 1987. On the frequency distribution of land slope. Earth Surface Processes and Landforms 12, 127e136. Papo, H.B., Gelbman, E., 1984. Digital terrain models for slopes and curvatures. Photogrammetric Engineering and Remote Sensing 50 (6), 695e701. Quinn, P., Beven, K., Chevallier, P., Planchon, O., 1991. The prediction of hillslope flow paths for distributed hydrological modelling using digital terrain models. Hydrological Processes 5, 59e79. Quinn, P.F., Beven, K.J., Lamb, R., 1995. The ln(a/tan b) index: how to calculate it and how to use it within the TOPMODEL framework. Hydrological Processes 9, 161e182. Ritter, P., 1987. A vector-based slope and aspect generation algorithm. Photogrammetric Engineering and Remote Sensing 53 (8), 1109e1111. Sharpnack, D.A., Akin, G., 1969. An algorithm for computing slope and aspect from elevations. Photogrammetric Engineering 35 (3), 247e248. 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), 323e334. Tarboton, D.G., 1997. A new method for the determination of flow directions and upslope areas in grid digital elevation models. Water Resources Research 33 (2), 309e319. Tribe, A., 1992. Automated recognition of valley lines and drainage networks from grid digital elevation models: a review and a new method. Journal of Hydrology 139, 263e293. Unwin, D., 1981. Introductory Spatial Analysis. Methuen, New York, pp. 212. Veregin, H., 1997. The effects of vertical error in digital elevation models on the determination of flow-path direction. Cartography and Geographic Information Systems 24 (2), 67e79. Wolock, D.M., McCable Jr., G.J., 1995. Comparison of single and multiple flow direction algorithms for computing topographic parameters in TOPMODEL. Water Resources Research 31 (5), 1315e1324. Zevenbergen, L.W., Thorne, C.R., 1987. Quantitative analysis of land surface topography. Earth Surface Processes and Landforms 12, 47e56.