Index | Search | About | Help
This chapter reviews and discusses the research that provides a context for this study. It is not meant to be an exhaustive review of the literature, but more an indication of the theoretical development of, and problems associated with, characterising surface form using Digital Elevation Models. More comprehensive reviews, collections and bibliographies include Pike (1993), Goudie (1990), (geomorphometry); Petrie and Kennie(1990), Weibel and Heller (1991), (Digital Elevation Models); and Wilkinson et al (1995), Unwin (1989), (fractals in the geosciences).
The form of surface model used for this entire study is that of the Digital Elevation Model (DEM). Although the term is used inconsistently in the literature (Burrough, 1986; Weibel and Heller, 1991) it is strictly defined here consistent with the terms of Burrough (1986), as a regular gridded matrix representation of the continuous variation of relief over space. Georeferencing of elements that comprise the surface are implicitly defined by the ordering of elevation values within the matrix. No additional or explicit representation of surface form is recorded. The contrast should be noted between such surface models and the Digital Terrain Model (DTM) which includes some additional explicit representation of surface form. Examples of DTMs include the Triangulated Irregular Network (TIN) (Peucker, 1978), digital contours with 'form lines' (Ordnance Survey, 1992), and the 'richline model' of Douglas (1988) that uses ridge, valley and form lines to define an elevation model.
The attraction of a simple matrix of elevation values (amenable to processing with procedural computer languages) was one of the reasons for its uptake in the early 1970s as a model suitable for landscape analysis (eg Evans, 1972). More recently, its continued widespread use as a model of surface form may be attributed to its easy integration within a GIS environment (Weibel and Heller, 1990,1991).
There are a number of conceptual problems that must be addressed when considering DEMs as models of surface form. Firstly, despite being a model of continuous surface form, a DEM is a set of discrete elevation measurements of what is usually an undifferentiatable surface. The fidelity with which the DEM models the true surface will depend on surface roughness and DEM resolution. Yet, evidence from fractal research (see section 2.5) suggests that there will always be detail at a finer scale than that measured at the DEM resolution. This suggests that all DEMs implicitly model at a certain scale implied by the grid cell resolution. Although the scale dependency of measurement has been recognised and even modelled by many authors (eg Garbrecht and Martz, 1993; Ackerman, 1993; Hodgson, 1995), it still remains hidden in many aspects of DEM analysis.
A second conceptual problem must be addressed in considering what each elevation value within a gridded matrix represents. It is possible to conceive of the matrix as a collection of elevation samples at point locations. If this assumption is made, the model is not one of continuous surface form without invoking an additional assumption about the spatial relationship between recorded elevation values. Although the interpolation between grid cells is rarely addressed by those using DEM analysis, the form of implicit interpolation can have significant effects on analytical results. Fisher (1993) showed that DEM derived viewshed areas could vary drastically depending on how adjacent elevation values were interpolated to form a continuous surface. Kumler (1994) suggested that a linear interpolation between adjacent DEM cells produced the smallest residuals when compared with TIN models of the same area. He also recognised that different interpolation procedures gave rise to different elevation estimates. This problem has led to a consistent (quadratic) form of cell interpolation to be adopted for this study where a continuous model of surface form is required.
The measurement of shape has received much attention in many branches of science ranging from, biology and palaeontology (for example, Thompson, 1917), to mathematics and imageprocessing (Serra, 1982). The systematic measurement of topography from cartographic sources can be traced back at least as far as the mid-nineteenth century (Cayley, 1859). Significantly, the development of a science of surface measurement has been accompanied by the development of surface storage and representation methods. Cayley's analysis of "slope and contour lines" was only possible (and useful) because of the introduction of accurate contour representations of relief in the nineteenth century (Imhof, 1982). This section considers the science of geomorphometry, or the science "which treats the geometry of the landscape" (Chorley et al, 1957).
A substantial part of twentieth century geomorphology, at least since Tricart (1947), has been devoted to the measurement and quantification of topographic form, although as Pike (1993) has noted, this has proved less successful than the characterisation of geomorphological process. The notable exception has been in the quantification of channel form at the drainage basin scale (Horton, 1945; Strahler, 1964; Shreve, 1966). The quantification of network form has allowed the empirical relationships between measurements to be investigated, often with claims of universality (eg Horton, 1945; Schumm, 1956; Hack, 1957, Shreve, 1974). A criticism of this approach has been that the indices used (eg stream order) or the empirical relationships observed (eg basin mainstream length ratio) are insensitive to topographic variation (Abrahams, 1985; Knighton, 1984). Further, such empirical observations may be as much a function of scale of sampling and measurement as they are of geomorphological variation (eg Church and Mark, 1980; see also section 2.5 below).
The drainage network forms only a part of an overall description of topography. Several authors have attempted to place network quantification within a wider context of topographic form. The early work of Cayley (1859) and Maxwell (1870) considered the topological links between channels, ridges and local extrema. Mark (1979) considered the topological properties of ridge networks. Werner (1988) attempted to formalise the relationship between the channel network and interdigitating ridge network. Wolf (1989,1991) provided a similar formalisation using graph theory. However these topological descriptions do not consider the geometrical properties of surface form, nor do they provide descriptions of continuous variation over a surface.
Evans (1972) made the distinction between attempts to quantify specific geomorphological features (specific geomorphometry) and general geomorphometry, "the measurement and analysis of those characteristics of landform which are applicable to any continuous rough surface" (Evans, 1972, p.18). The origins of this approach arise from the cartometric analysis of map data (eg Glock, 1932; Johnson, 1933; Wood and Snell, 1957, 1960) and the derivation of descriptive indices from them. Mark (1975a) suggests that what unifies this approach, is an attempt to quantify the concept of terrain "roughness". He suggests that the planimetric variation of roughness embodies two primary scales, namely the grain and texture. Grain refers to the longest significant wavelength of a topographic surface, texture, the shortest. Altimetric range has been characterised by various measures of relief, including local relief (altitude range over a local area, Smith, 1935); available relief (vertical difference between a former upland surface and degraded channel, Glock, 1932); and alternative measure of available relief (vertical difference between adjacent ridge and channel, Dury, 1951).
A number of significant problems arise from this form of quantification. Firstly, many of these early measures involved arbitrary or ambiguous definitions, such as 'significant wavelength', or 'former position of upland surface', making comparisons of measures derived by different authors difficult. Secondly there is significant dimensional overlap between measures, resulting in a redundancy in description. Attempts to remove descriptive redundancy using factor analysis (eg Lewis, 1969; King, 1968) or correlation (Doornkamp, 1968) have resolved this problem in part, although inadequate operational definition of indices make the interpretation of components difficult (Evans, 1972). A more systematic parameterisation of altimetric and planimetric variation was suggested and used by Evans (1972, 1979, 1980, 1984). He suggested that a unifying framework is provided by taking the first and second derivatives of altitude (slope, aspect, profile convexity and plan convexity). Dimensionally these measures, or morphometric parameters, are orthogonal and process related.
Evans (1979,1980) analysed the frequency distributions of morphometric parameters derived from DEMs of glaciated and mountainous terrains, and found they characteristically possessedhigh mean and variability in slope, high variability in altitude, and a positive skew in profile and plan convexities. Frequency distribution based characterisations were carried out by Speight (1971) who analysed slope distributions and Ohmori and Hirano (1984) who identified characteristic hypsometric distributions. Evans (1984) applied factor analysis to the 4 momental measures (mean, standard deviation, skewness and kurtosis) of the altitude and its four derivatives (slope, aspect, profile and plan convexity). He identified 10 'key variables' from this analysis that allowed him to characterise local variations in terrain. Herdegen and Beran (1982) and Pike (1988) both used slope and curvature to characterise different landforms by constructing landform 'signatures' (Pike, 1988) from these parameters.
The process of landform classification is not one that is considered by this work, but it should be recognised that a successful surface characterisation provides the information with which to classify landform. Pike (1988) suggests that a successful surface parameterisation is necessary for a flexible terrain taxonomy. Geomorphometric classification of terrain has tended to be either into 'homogeneous regions' (eg Speight, 1968, 1973, 1976, Dikau, 1989) or the identification of specific geomorphological features (eg valley heads by Tribe, 1990). Richards (1981) highlights a number of problems with geomorphometry that apply to both approaches. In particular, the problem of scale of both spatial extent and resolution make single objective classifications of landscape unfeasible.
The notion of scale in geomorphological analysis is an important one, and forms the basis of this study. Evans (1979) recognised that many of the morphometric parameters measured from a DEM would vary with the size of grid used to model the surface. He distinguished this (sampling) effect from the size of area measured, which he found to be far less scale dependent. Frank et al (1986) found that scale dependency in operational definition of some surface features was sufficiently high to preclude objective definition altogether. In a manual measurement context, Gerrard and Robinson (1971) found high levels of variation in slope measurement with different sampling interval. Garbrecht and Martz (1993) analysed the effect of DEM grid size on automated drainage network identification concluding that the grid resolution should be no larger that 5% of the drainage basin area, if effective characterisation is to take place. What becomes clear from these (and many other) studies, is that it is unwiseto ignore the scale based effects of sampling interval when characterising surfaces.
Probably the most widely examined and cited geomorphological feature extraction processes from DEMs concerns the derivation of hydrological and fluvial features from surface models. Consequently, it deserves special attention when reviewing existing surface characterisation work. The routing of water over the surface of a landscape represents a fundamental geomorphological process that is intimately tied to its form. The subdivision of the continuous surface into discrete hydrological units provides an important step in the geomorphological treatment of an elevation model. The literature on this specific process alone is large, Pike (1993) identifying over 100 refereed articles. Consequently, rather than provide a comprehensive (and repetitive) treatment of all papers in the field, this section attempts a categorisation of the techniques, identifying the major works in each category.
Prior to the use of automated hydrological characterisation from DEMs, features would either have to be measured directly in the field, or derived from secondary sources such as paper maps. Much of the rich literature on the pattern of drainage networks is based upon such extraction (eg. see Abrahams, 1984 for a review of techniques). If paper maps are used, drainage channels may be measured directly from the 'blue line' network, or inferred from contour crenellations. The accuracy of the blue line network clearly depends on the scale of the map source and quality of original surveying, but also the dynamism of the network itself. Networks with ephemeral streams may have particular symbolic representations (eg USGS 1:24000 topographic maps), or they may not be distinguished from permanent streams (eg, Ordnance Survey 1:50000 Landranger maps). The Ordnance Survey instructs its surveyors to represent the network at a constant hydrological state, at the "normal winter level" (Ovenden and Gregory, 1980).
If an alternative measure of drainage density is required, contour data may be examined so that channel form may be extracted. Many authors have questioned the accuracy of contour maps as sources of such information (eg. Drummond, 1974; Dunkerley, 1977; Mark, 1983) both because the accuracy of the contours themselves may be questionable and because there can be varying interpretations of the same contour data. Tribe (1990) noted that when two experienced geomorphologists were given the task of identifying the locations and areal extent of valley heads from a contour map, one consistently delimited a smaller head extent than the other.
From the multitude of literature on ridge and channel extraction, it is possible to characterise all the methods according to five dichotomous classification criteria. These are summarised in Table 2.1. Criteria have been selected that are general enough to group any method of network extraction from any type or dimension of surface model. They could, for example, equally be used to classify the extraction of a least cost path from a cost surface, or migration path from a space-time 'volume'. This classification procedure was first described by Wood (1990) and is analogous to Lam's (1983) classification of spatial interpolation methods. In fact, the relationship between feature extraction and spatial interpolation is not an arbitrary one. Both involve the characterisation of as much useful source material as possible, from which a target is derived. The essential difference between spatial interpolation and feature extraction is that the former involves extracting a target that is the same quantity as the source (eg elevation from elevation) whereas the latter involves a further transformation (eg. from elevation to morphometric feature type).
(i) Topological / Geometrical
An analogy can be drawn between this dichotomous classification of feature extraction techniques and that of Lam's (1983) point/area interpolation dichotomy. Both define the metric used for the source and target. For interpolation the source and target are either 0-dimensional point space or 2-dimensional area space. For drainage networkextraction, the target is either n-dimensional space or a more abstract topological 'space'.
It has long been recognised that surface models contain important topological information that characterises a surface. Early work by Cayley (1859) and Maxwell (1870) described how any contour map describing surface form contains a set of unique topological relationships between summits (local maxima), immits (local minima), bars (lowest point dividing two immits), and passes (lowest point dividing two summits). From the topological connectivity of these point locations, the line features of watercourses and watersheds, and the areal features of hills and dales, could all be delimited. The exact distribution in space of these features is not fully defined, but rather the number of links and their degree of connectedness.
This topological characterisation has been developed more rigorously by Warntz (1966), Pfaltz (1976) and Wolf (1984, 1989, 1992). Wolf used the more standard classification of surface topology (described in Chapter 5) by considering bars and passes as having the same topological form. More significantly, he modelled these connectivity relationships using graph theory. Thus the topology of any surface could be described using a weighted surface network of pits, peaks and passes connected by ridges and channels. This has been used as the basis for generalisation (Mackaness and Beard, 1993; Wolf, 1989) and stream network identification (Wolf, 1992).
The topological characterisation of derived features as a concise description of surface form has been examined by many authors. Werner (1988) looked at the topological relationship between ridge and channel pattern within a drainage basin and derived a number of consistent relations. The topological characteristics of channel networks themselves have been examined in great detail, the most widespread approach incorporating the Random Topology Model (Shreve, 1966).
Network topology may not provide sufficient information to characterise hydrology for many applications. It must also be recognised that GIS, the platforms for mostautomated extraction techniques, are not currently well equipped for processing topological surface data. Thus a distinction can be made between the techniques described above and the majority of extraction routines that are concerned with the geometry of hydrological features. Spatial location of surface features is required for much hydrological analysis (eg Lammers and Band, 1990), conversion to other data structures (Kumler, 1994) and realistic cartographic representation.
Conversion between topological and geometrical representation of drainage networks remains one of the outstanding problems of feature extraction. Wolf's surface network approach has limited application because it is difficult to attach spatial location to the processed weighted networks. Likewise, the fragmentation of networks produced by many of the geometric techniques (eg Peucker and Douglas, 1974; Jenson, 1985; Skidmore, 1990), makes the identification of topological relationships difficult. Two categories of solution to this problem have been adopted. Hutchinson (1989) describes a method of interpolating elevation using a `drainage enforcement algorithm' to force hydrological connectivity. This is done by identifying peaks and passes and forcing topological connectivity via channels that contain no pits. The alternative approach adopted by many authors (eg Band, 1986; Wood, 1990b) is to post-process the derived network to force topological connectivity. This may be in the form of line thinning (Skidmore, 1990), line joining and elimination (Wood, 1990b), combination of external data sources such as remotely sensed imagery (Riley, 1993).
(ii) Global / Local
Lam (1983) distinguishes between spatial interpolation routines that use some local neighbourhood in the interpolation of new values, and those that use the entire source. The same distinction may be made for the extraction of drainage features. As with spatial interpolation, the distinction is less a dichotomous one, more two ends of a continuum. For the purposes of classification, a distinction will be made here between three degrees of operation. Local extraction routines use a fixed window size that is less than the size of the entire surface mode. Quasi-local methods use an adaptivewindow that is in most senses local, but may change size according to the characteristics of the surface mode. Global routines require information from the entire surface model for extraction.
Peucker and Douglas (1974) describe how a local 2 by 2 window may be passed over a raster while flagging the highest cell in each window. Cells that remain unflagged after the entire DEM has been scanned are described as channel cells. This procedure (which may be conveniently reversed to detect ridge cells) is based on the assumption that all channel cells will have at least once neighbouring cell of a higher elevation. This, and similar methods are referred to as hillslope elimination methods by Wood (1990a). The majority of local routines are based on more precise morphometric properties referred to as concavity detection methods by Wood (1990a). An alternate method suggested by Peucker and Douglas (1974), Jensen (1985), Band (1986) all use a fixed local 3 by 3 neighbourhood to identify 2-dimensional profiles from which channel shape (and other morphometric features) may be identified. Fixed local windows have been used ranging from 2 by 2 (Douglas and Peucker, 1974) to 7 by 7 (Zhang et al., 1990).
Many authors recognise that using a fixed size window only allows channel features of a certain scale to be identified. Since this scale is arbitrarily defined by the resolution of the elevation model, it may not be always appropriate for feature extraction. Quasi-local methods allow the size of window used for sampling source data to vary according to the properties of the surface model. This variation can take several forms. Jensen and Domingue (1988) and Skidmore (1990) both describe methods similar to those of Jensen (1985) but where window size is allowed to expand in areas where profiles are ambiguous (eg. flat regions). Such methods have the potential to identify cross-sectional form at a variety of scales as the window expands. An alternative set of quasi-local techniques was adopted by Wood (1990b), where a variety of scales of surface form are examined using a range of window sizes.
Global feature extraction methods require the entire surface model to be examinedbefore any characterisation can be made. As with global spatial interpolation, every elevation value can potentially influence the outcome of characterisation. Collins (1975) used such a method by ordering all elevation values from highest to lowest before processing. This method, although computationally efficient disregards much of the natural spatial autocorrelation of real surfaces during extraction. Douglas (1986) shows that the method described by Collins may also produce incorrect results in certain circumstances. By far the most widely used group of methods of global hydrological characterisation is the flow magnitude method first described by Mark (1983a) and subsequently modified by many others (eg, Marks et al., 1983; O'Callaghan and Mark, 1984; Band and Wood, 1988; Jensen and Domingue, 1988; Bennett and Armstrong, 1989; Band, 1989). This method (described in more detail below) involves traversing the entire surface model accumulating predicted surface flow. Drainage channels can be identified as made up of locations that have exceeded an (arbitrarily) defined flow magnitude threshold.
(iii) Exact / Approximate
This dichotomous distinction is similar but not identical to the exact/approximate distinction of spatial interpolation made by Lam (1983). Lam's distinction refers to the similarity between coincident source and target values. An exact interpolator will honour all source values such that spatially coincident source and targets will have identical values. Approximate interpolators may result in deviation between the two. Hydrological feature extraction can be regarded in a similar context. Given that all methods use some kind of morphometric characterisation, it possible to distinguish between those that honour elevation values precisely and those that make some approximation about local elevation values.
Approximate characterisation of elevation may be for a number of reasons. Several methods apply some preprocessing of the elevation model before extraction such as mean filtering (eg. Mark, 1983a; Band, 1986). Alternatively, the elevation may bemodified to ensure hydrological connectivity. Hutchinson (1989) forced interpolation of contours to eliminate any apparent pits modeled by contour lines. Morris and Flavin (1984) describe how planimetric position of contour lines may be altered so that they are consistent with a surveyed blue line network. Thus a deviation is introduced between original elevation values and those used for the extraction of hydrological features.
Two of the outstanding problems with feature extraction concern the characteristion of form of flat regions and 'noisy' elevation models. Several authors have suggested that a tolerance value should be associated with each elevation before using it to characterise morphometry (eg. Marks et al., 1983; Skidmore, 1990). This is usually expressed as the vertical deviation required to define a real change in elevation. A third way of introducing approximate characterisation of surface form concerns those methods that model local neighbourhoods with some form of mathematical function or 'patch'. Evans (1979) suggested a bi-variate quadratic approximation of 3 by 3 neighbourhoods could be used to identify drainage features. Due to the overspecification of the six coefficients of the quadratic by the nine elevation values of the local neighbourhood, an approximate (but optimal) model of sampled elevations is produced. Mitasova and Mitas (1993) describe how a bivariate spline function can be used to detect channel features. Tension parameters are introduced that control the degree to which modelled values fit the original source data. Zhang et al (1990) apply a discrete cosine transformation to elevation, forming an approximate mathematical model of sampled points, before computing derivatives to find surface concavities.
Exact models honour all sampled points exactly during the extraction process. Zevenbergen and Thorne (1987) use a partial quartic expression to model a 3 by 3 local neighbourhood. Since this expression contains 9 coefficients, the model is specified precisely with no redundancy or deviation between sampled and modelled data. Other methods use the elevation values directly without patch modelling, pre-processing or using tolerance values (eg Jensen, 1985).
(iv) Direct / Indirect
In extracting hydrological features from an elevation model, two procedures may be adopted. Either the morphometry of the target can be identified and measured (eg. channel cross-section), or the target features may be associated with some other set of properties which can be in turn related to morphometry. There is a loose analogy here with the distinction between deterministic and stochastic interpolation methods. If the source and target values are both morphometric, it is possible to invoke a deterministic relationship between the two. If an indirect, intermediate rule must be invoked to relate the two, there is greater likelihood of indeterminacy.
Mark (1983a) defines a drainage channel as being made up of "points at which fluvial processes are sufficiently dominant over slope processes" (Mark, 1983a, p.169). If a measure of fluvial and slope processes can be made, so drainage channels can be (indirectly) identified. This indirect method of identifying drainage channels by the predicted volume of water flowing over their surface is the basis for many hydrological characterisations of surfaces (eg. O'Callaghan and Mark, 1984; Band and Wood, 1988; Band, 1989). Such methods are desirable when morphometry alone is not sufficient to characterise hydrology because other factors which may vary spatially have importance in hydrological modelling (eg Bathurst, 1986). Alternatively it may be that features may not be well defined morphometrically. For example drainage divides are relatively unambiguously defined in terms of their hydrological function, but may not be well expressed as morphometric ridge features. Conversely, channels may have a strong morphometric expression but have a widely varying hydrological role. Heavily dissected badlands in a semi-arid environment are particularly susceptible to such variation (eg loess plateaux described by Zonghu, 1986; Wood, 1990).
(v) Systematic / Recursive
The final dichotomy describes the way in which the feature extraction is applied spatially over the surface model. Systematic methods proceed in some orderly way thatis entirely independent of the characteristics of the source. Recursive methods are those which traverse the source in a pattern determined by the source itself. The parallel can be drawn with the gradual/abrupt distinction of Lam (1983). Gradual interpolation applies the same rules over the entire source whereas abrupt interpolation can involve the application of different rules at different points determined by 'barriers' in the source.
It should be noted that although the term recursive is used here, it does not necessarily follow that recursive algorithms are always used to implement recursive techniques. One of the most common recursive methods drainage network identification is to locate a drainage outflow point, recursively move upstream to the drainage divide and fall back to the outflow accumulating surface flow. This method has been implemented using recursive function calling by Marks et al (1983) and Band (1989), but non-recursively by O'Callaghan and Mark (1984).
A Classification of Hydrological Feature Extraction Techniques
| Author | Topol/Geom | Loc/Quasi/Glob | Exact/Approx | Direct/Indirect | Sys/Recursive |
| Cayley (1859) | T | Q | D | R | |
| Maxwell (1870) | T | Q | D | R | |
| Pfaltz (1976) | T | Q | D | R | |
| Wolf(1984,88) | T | Q | D | R | |
| Band (1986) | G | L | E | D | S |
| Jenson (1985) | G | L | E | D | S |
| Peucker & Douglas (1974)1 | G | L | E | D | S |
| Dikau (1989) | G | L | A | D | S |
| Evans (1979) | G | L | A | D | S |
| Zhang et al (1990) | G | L | A | D | S |
| Skidmore (1990) | G | Q | E | D | S |
| Speight (1968) | G | Q | E | D | S |
| Tribe (1990) | G | Q | E | D | S |
| Palacios-Velez et al (1986) | G | Q | E | D | S |
| Collins (1975) | G | G | E | D | R |
| Marks et al (1983) | G | G | E | D | R |
| Band (1989) | G | G | E | I | R |
| Band & Wood (1988) | G | G | E | I | R |
| Mark (1983a) | G | G | E | I | S |
| O'Callaghan & Mark (1984) | G | G | E | I | S |
| Jensen & Domingue (1988) | G | G | E | I | S |
| Bennett & Armstrong (1989) | G | L/G | E | D/I | S |
| Wood (1990b) | G | L/G | E | D/I | S |
Where no entries are given, the method can be considered in either grouping. Where two alternatives are given, the method uses both groupings to produce the final extraction.
Table 2.1 A classification of some hydrological feature extraction techniques
A consideration of the extraction of information from gridded elevation models cannot ignore the very similar process of information extraction from gridded digital imagery. In particular the image processing techniques of texture analysis and pattern recognition are considered here. This section will not consider the process of image segmentation (partition of homogeneous regions) or mathematical morphology (shape definition) which are more closely associated with specific geomorphology.
The analysis of texture within a digital image is closely allied to the geomorphometricmeasurement of roughness. There appear to be no formal definitions of texture (Gonzalez and Woods, 1992) although many authors make the distinction between statistical texture, and structural components of an image (eg.Haralick, 1979; Lu and Fu, 1978; Tomita et al, 1982). Statistical texture is the description of the highest frequency of variation within an image, while structure is the arrangement of 'textural primitives' that make up an image (Gonzalez and Woods, 1992). It can be argued that the distinction is simply one of scale, such that the texture at one resolution becomes the structure at another. A large branch of image processing is concerned with spectral analysis of images where a (Fourier) transformation of an image allows periodicities at a range of scales to be examined (eg. Bajcsy, 1973).
Early measurements of statistical texture include analysis of frequency distributions (Darling and Joseph, 1968), and measures of autocorrelation (Kaizer, 1955) but, according to Haralick et al (1973), failed to adequately define or characterise texture. They proposed a method of texture analysis involving the computation of the grey-tone spatial dependence matrix (co-occurrence matrix) from a digital image. This two-dimensional matrix is, in effect, a measure of the spatial autocorrelation properties of an image. From this matrix, various summary statistics can be calculated measuring distinct textural properties such as entropy and linearity. Weszka et al (1976) assessed the effectiveness of these measures by comparison with Fourier techniques, concluding that in general, spectral analysis was less discriminating than the co-occurrence matrix measures.
It should be noted that the co-occurrence matrix is a resolution dependent measure. The degree of spatial association between grey levels in many images will depend in part on their distance of separation (a property measured by the variogram or correlogram (eg. Goodchild, 1986)). Such scale dependency has been measured using Fourier transforms, but the trigonometrical modelling of surface form makes this technique most suitable for images that comprise continuous and regular periodicities (Gonzalez and Woods, 1992). More recent developments in the use of wavelet transforms (Graps, 1995) may prove more promising in that they allow discrete non-periodic scale dependencies to be measured. Gallant and Hutchinson (1996) use such wavelet transforms to reconstruct DEMs over a (controlled) range of scales.
Independently of the geostatistical tradition of characterising spatial dependency, there has arisen in the last 20 years, an alternative set of statistical tools that are able to describe aspects of scale dependency in spatial objects. Mandelbrot's persuasive treatises on fractals (Mandelbrot, 1975, 1977, 1982) have led to a prolific literature (over 100,000 articles according to Raffy, 1995) on scale-based characterisation of spatial and non-spatial phenomena. While much of the fractal literature is not strictly relevant to geomorphometry, there are a number of important ideas that have arisen from this field of study. It is ironic that one of the most persuasive aspects of Mandelbrot's work, according to Feder (1988), are the images created by Richard Voss (Voss, 1985) of "landscapes [that] look natural - one must believe that fractals somehow capture the essence of the surface topography of the earth" (Feder, 1988, p.3). In this context it is important to examine the relevant fractal literature in order to assess the contribution the theoretical advance has made to geomorphometry.
The essence of any fractal object is its self-similarity and can be defined as "a shape made of parts similar to the whole in some way" (Mandelbrot, 1986). More specifically, fractals can be defined in four contexts - deterministic or stochastic and self-similar or self-affine. Deterministic fractals can usually be defined by some mathematical function that is recursively applied over a wide range of scales. The result is a (usually geometric) object that contains scaled copies of itself. If these copies are identical bar their isotropic scaling and rotation, the fractal is termed self-similar (Mandelbrot, 1982). If the scaling process is itself anisotropic but linear, the fractal is termed self-affine (Loevejoy and Schertzer, 1995). Examples of deterministic self-similar fractals include the well known 19th century mathematical "monsters", the Koch Curve and Sierpinski Carpet (Mandelbrot, 1982). Although not directly relevant to geomorphometric characterisation, such fractals have been shown to be useful in developing spatial data structures to mimic surface spatial autocorrelation (Goodchild and Grandfield, 1983). Deterministic self-affine fractals have been largely used in image compression (Barnsely, 1988), botanical simulation (Govaerts and Verstraete, 1995) and atmospheric modelling (Lovejoy and Schertzer, 1995). More useful for the analysis anddescription of terrain is the concept of statistical self similarity and self-affinity. Here measurable statistical parameters of an object are repeated at all scales. Thus, although a statistically self similar object will not contain exact copies of itself at all scales, the object would contain no geometric indication of its scale (Goodchild and Mark, 1987).
The fractal behaviour of surfaces (statistical self similarity) can be detected with many of the geostatistical techniques used for interpolation. In particular the variogram can be used to identify self-similar scaling since the log of variance plotted against the log of lag will yield a linear relationship for self-similar form (Rees, 1995). It is this behaviour that has one of the most important implications for spatial science in general and surface characterisation in particular. The so-called Steinhaus Paradox (Goodchild and Mark, 1987) is a recognition that the measurement of length increases with increased accuracy (Steinhaus, 1960). A similar observation known as the Schwartz area paradox (Feder, 1988) has been noted in the scale-dependent measurement of surface area of a volume by triangulation. Any truly self-similar object would have indeterminate measurable geometry since at whatever scale of measurement used, there will always be a finer level of detail that can only be approximated (underestimated) by that measurement. Richardson (1961) showed empirically that this is often the case when measuring cartographic lines. Schertzer and Lovejoy (1995) have shown this is also the case for measurements such as vegetation cover or albedo from increasingly higher resolution satellite imagery. The importance of this observation is that it suggests that any measurement of length (or area or volume) of a fractal object, must be accompanied by the scale at which the measurement was made if it is to have meaning. The corollary of this is that any raster based model of a surface, which by definition is recorded at a specific resolution scale, is only a measure of part of a fractal object. It should be with caution, that any analysis carried out at a particular raster resolution is extrapolated to any other scale.
The degree to which length measurements change with scale was investigated by Richardson (1961) by plotting the log of line length against log of measurement sampling interval. The resultant straight line relationship indicated a consistent power 'law'. Similar empirical results were found for time series data by Hurst et al. (1965) who observed through 'R/S analysis' (Rees, 1995), a constant power law relationship between measurement variance and scale.Mandelbrot (1967, 1977) placed such observations in a fractal framework, by defining a fractional dimension D = log(n/n0) / log(l/l0) where l and l0 are two different length step sizes, and n and n0 are the number of these steps that measure the length of a line. For any truly self-similar object the value D will be a constant regardless of the two scales chosen for measurement. In practice D would be estimated by taking measurement at a variety of scales and taking the best fit regression. This may be measured directly using the so-called walking-dividers method (eg. Schelberg et al, 1982), from Richardson plots, where gradient of the best fit line = 1-D (Goodchild and Mark, 1987), power spectra (Burrough, 1981; Pentland 1984; Huang and Turcotte, 1990), or the variogram (Mark and Aronson, 1984, Polidori, 1995).
A number of methods of measuring the fractal dimension of surfaces have been widely used. It has been suggested that a truly fractal surface of dimension D will contain profiles of dimension D-1 (Voss, 1988, p.30). These profiles may be measured in the XY plane (ie contours) (eg. Shelberg et al, 1983), or with vertical profiles (Rees, 1995). A number of other techniques rely on box or cell counting (eg Voss, 1988, Goodchild, 1982). Alternatively, a 2-dimensional extension of the walking-dividers method has been proposed by Clarke (1985) and modified to look at local variation by de Jong (1995). A comparative assessment of these techniques was carried out by Roy et al (1987).
Given that fractal dimension can be measured, it is possible to ask (and possibly answer) the more important question as to whether topography is fractal. Implicit in the argument of Wood and Snell (1957) is that the relief measured over a range of sampling scales can be used to predict relief characteristics at other scales. Mandelbrot (1967) provides some evidence that the British coastline is fractal suggesting a characteristic range of D=1.2 to 1.3. This implicitly suggests the terrain it bounds may also be so. Rather indirectly, he also suggests that landscapes are fractal since fractal simulations "look right" (Mandelbrot, 1982). Goodchild (1980) points out that although the west coast exhibits a relatively constant value of D (1.25), the east coast varies between 1.15 and 1.31. Burrough (1981) examined a wide variety of environmental data including topography identifying a characteristic (profile) dimension of 1.5. Clarke (1986) found the fractal dimension of the USGS DEM of the San Gabriel mountains, California to be 2.19. Huang and Turcotte found the mean fractal dimension from USGS 1degree DEMs of Arizona to be 2.59 (Huang and Turcotte, 1989) and Oregon to be 2.586 with a standard deviation of 0.15 (Huang and Turcotte, 1990). Deitler and Zhang (1992) found the characteristic fractal dimension of Swiss alpine topography to be 2.43. Rees (1995) who examined topographic profiles of glaciated regions (Snowdonia, UK) and ice masses (Svalbard) using variograms, found D to vary between 1.1 (Snowdonia) through 1.3 (lowland UK terrain) to 1.4 for large ice masses.
What becomes apparent when comparing empirical results is that a single fractal dimension is not a very useful measure. It may fail to distinguish between contrasting topographies (eg. Huang and Turcotte, 1989, 1990; Fioravanti, 1995). The measure itself can depend on the procedure used to derive it (Roy et al, 1987). There is some degree of disagreement over what constitutes a 'typical' topographic value (2.3 - Shelberg et al (1983); 2.5 - Outcalt and Melton (1992); 2.2 - Voss (1988)). Most significantly, most empirical work suggests that landscapes do not possess a single fractal dimension, but a variety of values that change with scale. Of the 17 surfaces examined by Mark and Aronson (1984), 15 were said to show at least 2 characteristic values. At short lags (<1km) values were typically around 2.4, at larger lags, D appeared to be much higher (up to 2.8). Roy et al (1987) attribute this in part to sampling bias and possible anisotropic effects. Hakanson (1978), although not measuring fractal dimension explicitly noted a smoothing in lake shoreline as measurement size decreases. Outcalt and Melton (1992) noted a systematic spatial variation in D within DEMs.
Systematic variation in values of D can be interpreted in two ways. Polidori (1995) suggests a typical log-log variogram measured from a DEM will have three 'zones'. At the medium lags, self-similar (linear) behaviour will be evident, from which a fractal dimension could be derived. At larger lags, the variogram slope is unpredictable as the sample size rapidly decreases. At small lags, variance (and so fractal dimension) decreases more rapidly as the generalisation effects of DEM resolution begin to dominate. Alternative causes are attributable to geomorphic process. It is attractive to assume that since geomorphic process varies with scale, so might the surface form produced by such process. Mark and Aronson (1984) suggest the scales at which fractal dimension changes represent the change in dominance from one form of geomorphic process to another. Certainly, there is additional empirical evidence forcharacteristic scales or breakpoints in fractal dimension. For example Dietler (1995) found Swiss alpine topography to be effectively made up of self-affine units of size 5km. Gilbert (1989) noted a general linear trend, but with local non linear stretches over limited scales.
While the overwhelming weight of evidence suggests that topography is at best fractal over a limited range of scales, the concept forces the recognition of two important points. Firstly, the scale at which surfaces are modelled and analysed will in part determine surface characterisation. Secondly, the very absence of scale independence is itself a useful diagnostic, be it through systematic change in fractal behaviour, or as characteristic scales of change. More recent development in the field of multifractals is an attempt to recognise and model scale and attribute dependencies. Lovejoy and Schertzer (1995) suggest that traditional (mono)fractal analysis is inappropriate for modelling what they regard as multifractal topography. They suggest that fractal dimension systematically decreases with altitude. Again this has an appeal to geomorphologists who acknowledge different upland and lowland process and form. Multifractal characterisation replaces a single scale invariant value with a scale and attribute dependent function D(q). Fioravanti (1995) uses such D(q) functions to characterise the texture of SAR images and shows the technique to be far more discriminating than either a single fractal dimension or other traditional texture measures (such as the co-occurrence matrix measures). Allied to the use of multifractal analysis is the use of wavelet transforms to specifically map discrete scale dependencies within the frequency domain (Graps, 1995).
