John C. Gallant and Michael F. Hutchinson
Grid resolution of digital elevation data profoundly influences both the spatial pattern and the frequency distribution of derived topographic attributes, such as slope and specific catchment area, and therefore influences environmental models built on such attributes. A model of topography that explicitly incorporates scale is a prerequisite for an understanding of the origins and nature of the scale dependence of topographic attributes. We propose a model of surface topography that directly represents features in the landscape at a range of scales. This model allows generalisation and refinement of elevation models and also characterises the morphology - shapes, sizes and orientations - of features in the surface. Features are derived from a grid DEMusing a correlation detection process with successive refinement from broad scales to finer scales.
Topographic analysis (or terrain analysis) is the quantitative analysis of topographic surfaces with the aim of studying surface and near-surface processes. A number of topographic attributes can be calculated, including specific catchment area, slope, aspect and plan curvature (contour curvature). Slope and specific catchment area are key variables in hydrology and are used to predict spatial patterns of soil water content and erosion (Beven and Kirkby, 1979; Moore et al., 1991; Moore and Wilson, 1992; Moore etal., 1993d; Moore, 1995; Wilson and Gallant, 1996). Solar radiation estimation is based on slope and aspect, modified by topographic shadowing (Moore et al., 1993c; Gallant and Wilson, 1996a). The spatial distribution of soil physical and chemical properties can be modelled within uniform geological settings using a combination of topographic attributes (Moore et al., 1993a; Gessler et al., 1995). Vegetation distribution, which responds to water, light and nutrient availability, can be modelled using combinations of topographic attributes which capture much of the landscape-scale variability of these parameters (Moore et al., 1993c). In short, topographic analysis provides the basis for a wide range of landscape-scale environmental models which are used to address both research and managementquestions.
It is now widely recognised that topographic analysis results are sensitive to the resolution of the source generalised. This affects all topographic attributes but in varying ways. The resolution-dependence of slope and specific catchment area have been the most intensively studied because of their regular application in hydrological modelling (Moore et al., 1993b; Zhang and Montgomery, 1994; Quinn et al., 1995). In these and similar studies the primary question to be addressed has been "What grid resolution should I use for a particular modelling exercise?", and some useful answers have emerged. But it is worth taking a step back and asking why topographic attributes are scale dependent and what that can tell us about the topographic surface. What are the characteristics of topographic surfaces that induce the observed scale dependence of topographic attributes? Is the scale dependence consistent across scales? Do topographic features have similar shapes at different scales? Is there a different kind of organisation of the landscape at different scales? In this paper we will explore some approaches to studying scale dependence of topographic surfaces.
Our first observation is that grid resolution is not a particularly appropriate representation of scale. When we subsample an elevation grid to obtain another grid at coarser resolution, we are not only removing fine scale features of the surface (the intended change) but also changing the number of square cells into which the surface is divided. If grid resolution is used to study scale dependence of topographic attributes, the analysis is complicated by the different number of samples obtained from each resolution. Furthermore, specific catchment area is generally computed by accumulating cell areas from adjacent cells, and this network of connections is changed when the grid resolution ischanged. The minimum catchment area resolvable using the usual flow accumulation algorithms is also dependent on grid size, so grid resolution introduces a number of complicating artifacts to the analysis of scale dependence. If we wish to study the scale properties of a topographic surface, rather than the effect of grid resolution per se , it would be best to use a method that dealt with scale directly.
One technique for studying scale effects is spectral analysis, which provides information on relative amounts of variation at different wavelengths or spatial frequencies (Pike and Rozema, 1975; Gallant et al., 1994): wavelength is approximately equivalent to spatial scale. Spectral analysis as applied to topographic data is open to criticism on several fronts. Firstly, it assumes stationarity of the signal - the mean, variance and higher order moments should all be independent of location. This is clearly not true for topographic data, and in fact the spatial variation of such parameters is of considerable interest. The second difficulty with Fourier transforms is that they use oscillatory waveforms (sine functions) as the basis functions into which the sample data is decomposed. If the sine functions are not a good representation of the fundamental shapes occurring in the landscape then good localisation in scale is not possible. Non-sinusoidal shapes produce harmonics at shorter wavelengths which can overwhelm the contribution of smaller-amplitude features at those wavelengths. This also implies that a single feature in the landscape (if such a thing exists) might be represented in the Fourier transform by a substantial number of sinusoidal components, and conversely a single component in the Fourier transform contains contributions from a number of surface features. A basis function which better represents the fundamental shapes in the landscape would provide more meaningful representation of scale.
If we wish to overcome the first problem of non-stationarity we might use a small window which is moved across the data set to perform a spectral analysis at each window location. The main difficulty with this approach is that resolution must be traded off against location resolution, since a large window is needed to obtain good resolution in frequency (or wavelength) while a small window is needed to obtain good resolution in location. A more practical alternative is the wavelet transform which offers the best of both worlds by effectively using a window length that changes with wavelength. Short wavelengths use small windows thus giving good location resolution for fine-scale phenomena, while long wavelengths use large windows giving good wavelength resolution for broad-scale phenomena (Chui, 1992). Wavelet analysis uses a single basis function that is localised in both position and frequency and can translated and dilated to cover the entire position-frequency plane. This contrasts with Fourier analysis which uses many basis functions each at a different frequency. Wavelet analysis could produce a power spectrum at every point in the data set, which overcomes the non-stationary/non-local problem of the Fourier transform but results in a very large data set which then has to be processed further to obtain interpretable results. Furthermore, wavelet analysis still does not address the second problem of using oscillatory basis functions. Although there is a wide choice of basis functions, the wavelet transform requires the basis function to have zero mean which implies a degree of oscillation.
The representation we propose consists of a superposition of features at various scales. A surface can be constructed by introducing broad-scale features first and refining the surface by adding finer features onto the broader features. The features for such a representation should have at least the following properties:
The feature chosen to develop the method has an elliptical plan form and a smooth polynomial profile form as shown in Figure 1.
The parameters of a feature are its location, length, width, orientation and height: note that height can be negative. This is the simplest feature with sufficient flexibility to represent topographic surfaces. Figure 2 shows an idealised small catchment synthesised using 9 of these features to demonstrate that plausible topographic surfaces can be constructed using such features.
Synthesis is the "forward" problem and is reasonably straightforward. How do we solve the "inverse" problem: is it possible to take an existing surface (typically represented as a grid DEM) and decompose it into a set of features that closely approximate the original surface?
Positive wavelet analysis (Watson and Jones, 1993) is similar to wavelet analysis in that it is based upon a single function which is translated and dilated, but it replaces the oscillatory wavelet with a positive pulse. Because this positive wavelet has a non-zero mean the inverse wavelet transform does not exist - it is not possible to reconstruct the original data from the wavelet coefficients computed from the forward transform. However, by using a process of correlation detection, individual features in the form of translated and dilated copies of the basis function can be detected in a surface which can then be used to reconstruct an approximation of the original surface. It is thus perfectly suited to the problem of decomposing surfaces using features.
The correlation detection algorithm finds the location, scale, shape
and orientation of the feature giving the largest correlation value
which is the real 2-dimensional wavelet transform of the function
g(u, v) using the wavelet F(). The feature is located at x and y, L is
its length, Lw is its width and theta is its orientation; note that the
width w expresses the width of the ellipse as a fraction of length L and
takes values between 0 and 1. The scale of a feature is L w^½ which is
proportional to the square root of its area L^2 w pi/4.
The amplitude of a detected feature is
C is a constant dependent on the wavelet shape
The detected feature is subtracted from the surface and another feature is detected from the residual surface. The parameters of features detected from the surface are distorted by the effect of other overlapping features, so the parameters of overlapping features are jointly optimised using non-linear least-squares (Marquardt, 1963; Press et al., 1989) with sum of squared differences between data and modelled elevations as the objective function.
The iterative decomposition detects the largest magnitude features first (which in topographic data tend to also be at large (i.e. broad) scales) followed by progressively smaller features. This leads to a relatively sparse representation of the surface being analysed provided the shape of the positive wavelet is a good match to the shape of the features in the sample. It is likely that some of the small scale features will be corrections for errors in the shapes of large scale components rather than real features in the data. If this is the case, some additional parameters could be included to modify the profile shape of the features to improve the fit to the surface.
In practice, scanning the entire 5-dimensional correlation surface T(x, y, L, w, theta) for every new feature is computationally very expensive. A simpler search has been developed which searches over a small subset of the space that is highly likely to contain the maximum correlation value.
Figure 3 shows contours of a 1 km x 1 km section from a 20 m DEM of the Brindabella Range near Canberra, Australia. The main features are the hill in the northwest corner, the general easterly slope on the eastern half and the drainage line running to the east in the southern half of the square.
Figure 4 shows a reconstruction of the surface using the first four features detected by the correlation detection algorithm. Note that these four features have captured the dominant broad-scale characteristics of this piece of landscape. Further refinement of the surface to 20 features is shown in Figure 5. The shape of the surface has been essentially captured but some fine scale detail is still missing. A decomposition to an average accuracy of 1 m (largest difference about 3 m) requires about 50 features for this patch.
The patch of Figure 3 is too small to undertake a very meaningful scale analysis. The full 5.8 km x 3.7 km DEM from which the patch was taken has been decomposed to 690 features using a target average residual of 1 m. The surface was then reconstructed at four levels of generalisation by excluding features smaller than 200, 500, 1000 and 2000 m scale. All surfaces were reconstructed as 20 m resolution DEMs. The 2000 m scale surface contains only 12 features, the 1000 m surface has 37, the 500 m surface has 107 and the 200 m surface has 471 features. The contours of the original and the three most generalised surfaces are shown in Figure 6. The surface reconstructed using features larger than 200 m scale (not shown) differs from the original surface only at quite fine scales.
The five surfaces were then analysed using the TAPES-G terrain analysis package (Moore, 1995; Moore et al., 1993b; Gallant and Wilson, 1996b) to compute contributing area using a multiple flow direction algorithm (Freeman, 1991; Quinn et al., 1991).The spatial pattern of contributing area from the original DEM and the 1000 m generalised surface are shown in Figure 7 with darker shades representing larger contributing areas.
Figure 7 shows that much of the drainage network remains intact when features smaller than 1000 m scale (less than 1km^2 or 2500 grid cells) are removed, but that there is a much lower proportion of small contributing area values due to the absence of fine-scale division of the landscape into small catchments.
Figure 8 shows the frequency distributions for all 5 surfaces (original and four levels of generalisation) which also shows a progressive decrease of the proportion of small contributing area values as the surface becomes increasingly generalised. The most generalised surface shows a significantly different frequency distribution because it fails to separate two catchments on the eastern side of the DEM (see Figure 6. Previous analysis of the effect of scale on contributing area distributions using changes in grid resolution (Moore et al., 1993b; Zhang and Montgomery, 1994) have not been able to show the effect of scale at small contributing area because of the way minimum contributing area increases systematically with grid size.
The positive wavelet decomposition presented here is a useful tool for analysis of scale dependence in topography. It explicitly identifies features at a range of scales, allowing generalisation of a topographic surface to allow detailed study of the effect of scale on topographic attributes without introducing artifacts due to changes in grid resolution. The shapes and orientations of features identified in the landscape may also be useful for characterising landforms and delineating regions of contrasting surface structure.