PhD Thesis

Jo Wood's Home Info. Science Home GI Staff GI Research GI Teaching & Learning

Chapter Six - Results

This chapter contains an assessment of the DEM characterisation tools that have been developed in the previous chapters. This takes three forms; a theoretical evaluation of a technique's validity; a calibration of the visualisation tools; and an investigation into a technique's utility through application.

A theoretical evaluation has, in part, already been discussed in the previous three chapters. Since many of the techniques used involve the quadratic approximation of local surface patches, the implications of this form of modelling will be assessed in greater detail. The more abstracted measurements and their visualisation (for example, the lag diagrams discussed in section 4.5) require some form of calibration where the geomorphometric analysis yields measurements that are not intuitively obvious. This is most effectively carried out by applying such measures to surfaces with simple and known properties. These can then be applied to surfaces whose properties are not known. Finally the utility of the techniques and methodology proposed is illustrated with the analysis of 'real' topographic surfaces.

6.1 Control Surfaces.

The evaluation of characterisation tools involves using some control surfaces whose properties are known or can be systematically altered. This section outlines their derivation and utility.

6.1.1 Uncorrelated Gaussian Surfaces.

If one of the objectives of surface characterisation is to identify spatial structure, it is useful to make comparison with surfaces that have no spatial structure at all. For the purpose of this study, a GRASS module, r.gauss.surf was written (see Appendix) that generates an uncorrelated random surface with a Gaussian (normal) distribution. The resolution, size, mean and standard deviation can all be controlled from the program.

The Box-Muller method is used for generating the random deviates with Gaussian distribution (see Press et al, 1992, pp.288-290). This involves transforming two uniform random deviates between 0 and 1 (x1 and x2) as follows,

6.1.2 Object Surfaces.

Alternatively, a control may be used that has a high degree of spatial structure. Such surfaces are particularly useful for calibrating lag diagrams (for example, see the discussion on measures of spatial autocorrelation, section 4.5.1).

To create simple geometric objects within GRASS, a module r.xy (see Appendix) was written that creates two raster surfaces based on the coordinates of each cell location. An origin is set at the bottom left corner, with unit coordinate displacement,

These surfaces may be manipulated using map algebra (Tomlin, 1989; Shapiro and Westervelt, 1992), to create simple geometrical surface objects. For example to create a cylinder of height 20 units and of radius 50 units centred over a 201 x 201 cell raster, the following map algebra expression was issued using the GRASS module r.mapcalc,

To create a 20 x 20 x 20 unit cube in the top left corner of a 201 x 201 cell raster,

6.1.3 Fractal Surfaces

The third form of control surface used in evaluation allows a degree of spatial structure to be controlled. It was originally anticipated that this could be achieved by creating polynomial surfaces of a user-defined order. However, since much of the surface modelling involves quadratic fitting, there is a danger of producing circular evidence. That is, there may be a tendency for quadratic models to fit other polynomial surfaces more effectively than, for example, trigonometric or 'real' surfaces.

Fractal surface generation was adopted as a more realistic but controllable model of topographic variation. By changing the fractal dimension of the surface, the degree of spatial structure can be varied. There are numerous methods of generating fractal surfaces (see section 2.5), but the one adopted here uses the spectral synthesis approach described by Saupe (1988), pp.105-109. This technique involves selecting scaled (Gaussian) random Fourier coefficients and performing the inverse Fourier transform. It has the advantage over the more common midpoint displacement methods which produce characteristic artifacts at distances 2n units away from a local origin (Voss, 1988). More importantly for this work, this technique has been modified so that multiple surfaces may be realised with only selected Fourier coefficients (see r.frac.surf in the Appendix). The result is that the scale of fractal behaviour may be controlled as well as the fractal dimension itself. Figure 6.1 shows a fractal surface of dimension D=2.10 rendered with (a) 1/8th, (b) 1/4, (c) 1/2, and (d) all of the Fourier coefficients transformed.

Figure 6.1Figure 6.1 - Fractal surfaces created using Spectral Synthesis. Each image shows partial Fourier transformation of (a) 1/8, (b) 1/4, (c) 1/2, and (d) all, of the scaled Gaussian coefficients.

6.2 Terrain Models

The DEMs representing real terrain were all selected from the Ordnance Survey 1:50k coverage of the UK. Three areas for study were selected that show contrasting topographic variation. It should be remembered that the objective of this study is not to provide a comprehensive characterisation of the UK landscape, but to demonstrate the utility of the characterisation tools themselves.

6.2.1 Lake District

The following Ordnance Survey DEM tiles were used (with National Grid coordinates indicating extent of the coverage),

Each tile consists of 401 x 401 cells of size 50m x 50m, giving a 3600km2 coverage of the major part of the Lake District. Various 'sub-windows' were selected for the evaluation below (and previous chapters).

The area represented is a mountainous previously glaciated environment ranging in altitude from sea level to 977m (England's highest mountain). The central upland region (NY20) is dominated by Borrowdale Volcanics, with Skiddaw Slates to the north (NY22) and Silurian shales and sandstones to the south (SD28).

6.2.2 Dartmoor

The following Ordnance Survey tiles covering central Dartmoor were selected,

The area is entirely within the surface outcrop of the Dartmoor granitic intrusion giving rise to a characteristic upland 'moor and tor' topography. The extent of the selected DEM is shownin Figure 6.2 (15 x 14 km).

Figure 6.2Figure 6.2 - Dartmoor study area.

6.2.3 Peak District

The third area selected for study covers a wider 40 x 41km area around the south of the Peak District (40 x 41 km). The Derwent valley runs north-south across the map from the moorland in the High Peak area to the south of the region near Matlock (Figure 6.3). The geology is dominated by Carboniferous Limestone to the west of the Derwent and Millstone Grit to the east and north.

Figure 6.3Figure 6.3 - Peak District study area.

6.3 Assessment of Quadratic Approximation Tools

6.3.1 Non-spatial analysis of residuals.

Even at the 3 x 3 kernel size, a quadratic function cannot precisely model the 9 values used in its construction. The question arises therefore, does the approximation lose vital surface information, or worse, does it introduce systematic artifacts into the surface model?

To test this question, the pattern of residuals between the quadratic parameter f (elevation at centre of the kernel) and the measured elevation from the DEM was examined. Three surfaces were used for analysis, each with a different degree of spatial structure. Two fractal surfaces of 200 x 200 cells were created, one with a low fractal dimension (D=2.01), the other with a high fractal dimension (D=2.9). A third uncorrelated Gaussian surface was created for comparison. Summary statistics for all three layers are shown in Table 6. 1.

Table 6.2 shows the mean and standard deviation of the residuals for various kernel sizes. As expected, the standard deviation increases with larger kernels (where the quadratic function becomes increasingly overspecified). The nature of the relationship between residuals and kernel size is shown in Figures 6.4 and 6.5.

Figure 6.4Figure 6.4 - Kernel - residual size relationships for three control surfaces.

Figure 6.5Figure 6.5 - Kernel size - mean residual relationships for three control surfaces.

For the surfaces with the greatest spatial autocorrelation (fractal D=2.01), there appears a well defined (R2=0.996) linear relationship between kernel size and the magnitude of the residual between quadratic model and original elevation. This suggests that the degree of approximation of 'smooth' surfaces can be controlled entirely by kernel size. As the surface becomes 'rougher' there appears a well defined (R2=0.902) log-linear relationship between residuals and kernel size. The transition from linear to log-linear relationship suggests that for rougher surfaces there is a degree of 'diminishing returns' with increasing spatial approximation. That is, most approximation occurs at smaller spatial scales with little further generalisation introduced by larger kernels. This relationship is exemplified by the case of the uncorrelated surface where there is little control over residual magnitude beyond a kernel size of 7 x 7 cells.

The relationship between kernel size and residual mean (Figure 6.5) should indicate any systematic over- or under- estimates of elevation. The magnitude of maximum residual mean is less than 0.5% of the residual standard deviation, and so in practice, will have very little systematic effect on elevation estimates. It is worth noting however, that the degree of over-or under-estimation of elevation appears to be dependent the topographic form as well as kernel size. For example, a DEM representing a single valley system will tend to produce consistently higher overestimates of elevation as kernel size is increased. Likewise, a tendency to regress towards a mean elevation will produce increasingly large underestimates of elevation for a DEM representing a single peak.

6.3.2 Spatial analysis of residuals

As with the discussion of DEM uncertainty (see Chapter 3), global measures of elevation variation are not sufficient to understand the pattern of residuals, since they are likely to vary over space. Figure 6.6 shows selected 'residual maps' for the quadratic approximation of the two fractal surfaces at different scales.

Figure 6.6Figure 6.6 - Residual maps for four quadratic approximations. Colour scales are normalised so that pure red/blue is 1 standard deviation from mean. (a) 3x3 approximation of fractal D=2.01; (b) 25x25 approximation of D=2.01; (c) 3x3 approximation of fractal D=2.90; (d) 25x25 approximation of fractal D=2.90.

All four images in Figure 6.6 show some form of spatial clustering of residuals, that not unexpectedly corresponds with a change in kernel size. To describe the spatial pattern in more systematically, Moran's I lag diagrams (see section 4.5) were visualised. Figure 6.7 shows two selected examples corresponding to residuals produced by 25x25 and 13x13 cell kernels. Red indicates positive autocorrelation, blue negative. The magnitude of the statistic is additionally shown with selected profiles.

Figure 6.7Figure 6.7 - Spatial autocorrelation lag diagrams and profiles for two residual surfaces.

The clustering seen in Figure 6.6 is described by the high spatial autocorrelation at low lags. Spatial autocorrelation is reduced as lag increases, reaching a minimum value before becoming approximately uncorrelated. In both cases, the innermost uncorrelated (white) band corresponds to approximately half the kernel size, while the maximum negative value is at approximately 2/3 the kernel width. There is also a degree of anisotropy shown by the lag diagrams as 'diamond' rather than circular bands of equal autocorrelation. This is due to the square shape of the kernel that was used without any form of distance weighting. These findings suggest that filters that do not have any kind of distance weighting may well have tobe treated with caution, especially when processing surfaces with low spatial autocorrelation at the scale of the filter.

6.3.3 Quadratic approximation as a surface generalisation tool

The effectiveness of quadratic approximation as a generalisation process was examined by comparing the smoothing effect of mean convolution filtering (the most common form of surface smoothing) with that of quadratic approximation.

Figure 6.8 shows the visual effect of quadratic approximation as a generalisation process, applied to the central Lake District (35km x 35km). Kernel sizes range from 7x7 (350m) to 25x25 cells (1.25km).

Figure 6.8Figure 6.8 - Effect of quadratic approximation on terrain roughness. (a) original Lake District DEM; (b) 7x7 quadratic approximation kernel; (c) 13 x 13 kernel; (d) 25 x 25 kernel.

The effect of non-distance weighted quadratic smoothing on the frequency distribution properties of three surfaces was considered. The range and standard deviation of (i) an uncorrelated Gaussian surface, (ii) the Lake District DEM pictured in Figure 6.8; and (iii) the Leicestershire DEM were examined after filtering at all scales from 3x3 to 25x25 cells. Comparison was made with a non-distance weighted mean filter (using the GRASS program r.mfilter) at the same range of scales. Figure 6.10a-f show the relationship between elevation dispersion and kernel size. Note also, the non-linear smoothing effect of the Gaussian surface that was also observed in section 6.3.1.

As would be expected, the larger the filter, the greater the reduction in range and variance of elevation values. However quadratic smoothing retains the original global dispersion characteristics far more effectively than mean filtering. This is reinforced by visualising the smoothed surfaces. Figure 6.9 shows the south east corner of the Lake District DEM (Windermere/Ambleside) smoothed using both the 25x25 quadratic and mean filters. Quadratic approximation appears to preserve characteristics at the certain scales while smoothing detail at finer scales. Mean filtering would appear far more scale insensitive in its smoothing effects.

Figure 6.9Figure 6.9 - Comparison of (a) quadratic approximation and (b) mean filtering as a smoothing process. Figure 6.10Figure 6.10 - Effect of filtering on the range and standard deviation of surfaces.

6.4 Calibration of Spatial Dependence Measures.

Section 4.5 described how several tools were developed to visualise and measure the spatial structure of DEMs. This section evaluates their utility by 'calibrating' them with surfaces with known properties.

It was suggested that a lag diagram of Moran's I spatial autocorrelation index could reveal something of the relationship between spatial dependence and scale. To evaluate this, multiple realisations of a fractal surface were generated. Eight surfaces of fractal dimension D=2.10 were created, and eight of D=2.90. A lag diagram of Moran's I was generated for each surface. Figures 6.12 and 6.13 show the results for surfaces with D=2.10 and D=2.90 respectively. The bi-polar colour scheme used for all Moran diagrams is shown in Figure 6.11.

Figure 6.11Figure 6.11 - Bi-polar colour scheme used for Moran's I lag diagrams.

Figure 6.12Figure 6.12 - 8 realisations of a fractal surface (D=2.10) with corresponding Moran's I lag diagram shown below.

The notion of a distinction between texture and structure is a common one in image processing (eg Gonzalez and Woods, 1992) and is useful to consider here. All the surfaces shown in Figure 6.12 have similar (fine scale) textural characteristics (determined by the fractal dimension), but structurally (coarse scale) they vary. That is, due to the random generation of Fourier coefficients, a surface may describe on a low frequency peak, pit, ridge, channel or combination of features. What distinguishes texture from structure is simply the scale at which variation is considered. All the images in Figure 6.12 show similar fine scale variation, in that the central portions of all lag diagrams are similar. Variation tends to occur away fromthe centre seen as concentric zones of equal spatial autocorrelation. The most obvious feature of these larger lag measures, is the ability to represent anisotropy. The orientation of any coarse scale valley or ridge systems is picked out by elongated bands of equal spatial autocorrelation. In several cases such bands tend to 'twist' with increasing lag suggesting a change of orientation with scale. It should be noted that such variation would not be detected using the traditional (one-dimensional) variogram. Surfaces without a dominant ridge or valley system tend to produce lag diagrams with a greater proportion of positively autocorrelated measurements.

Figure 6.13Figure 6.13 - 8 realisations of a fractal surface (D=2.9) with corresponding Moran's I lag diagram shown below.

The eight sets of images shown in Figure 6.13 show a marked contrast with the surfaces of lower fractal dimension. Generally, the magnitude of either positive or negative autocorrelation is lower (images appear 'paler'). The characteristic diameter of high spatial autocorrelation is much smaller for these rougher surfaces. Coarse scale structure and anisotropy is still revealed, albeit less strongly.

It would appear that Moran's I lag diagrams are useful in detecting structural variation and anisotropy, and the relative dominance of textural and structural components of surface form. Although the terms texture and structure are used, these diagrams show the scale based continuum between the two rather than forcing an artificial dichotomy between them. To investigate the distinction between textural and structural variation further, a second set of fractal surfaces was produced. Using a single set of scaled Gaussian Fourier coefficients, selected coefficients were transformed. The effect was to generate a series of similar surfaces but each with increasingly high frequency detail added (see Figure 6.1 for an example).

Figure 6.14 shows two surfaces from the series, one with 1/8th of the Fourier coefficients transformed, the other with all coefficients. Surface texture is emphasised by combining elevation and local relief as a hue-intensity map. The Moran's I lag diagrams for both surfaces are shown below.

Figure 6.14Figure 6.14 - Fractal surface (D=2.1) with (a) 1/8 and (b) all Fourier coefficients. Their Moran's I lag diagrams are shown below as (c) and (d).

The two lag diagrams appear very similar, suggesting that they are not suitable visualisations of high frequency textural variation. Indeed, visualising the original surfaces (Figure 6.14a and b) is much more useful. This demonstrates that such lag diagrams are dominated by, and therefore most suitable for examining, structural variation. The change in texture between the two surfaces occurs only in the immediate (dark) central portion of each image.

To examine finer scale textural variation , co-occurrence textural measures were investigated (see also section 4.5.2). The co-occurrence matrices for the two surfaces shown in Figure 6.14 were calculated and textural measures calculated from them. A visualisation of the two matrices is shown in Figure 6.15 and the derived texture measures are recorded in Table 6.3.

Figure 6.15Figure 6.15 - Two co-occurrence matrices for fractal surface with (a) 1/8th, and (b) all Fourier coefficients transformed. Elevation is quantized into 100 categories, lag is one raster unit.

The global textural measures do appear to discriminate between the two surfaces although it is not clear which distinct texture properties each measure. The co-occurrence matrix itself also appears to vary between the surfaces with the rougher surface producing the more disperse matrix. However, it should be considered that both the matrix and the derived indices are measured at a single lag (unit offset in this case). It would be expected therefore, that they should show different patterns at different scales.

To determine any scale dependency in textural characteristics, lag diagrams were calculated for two of the measures - contrast and inverse distance moment. The others were excluded either because they are scaled variations of the same property or because they showed no scale dependency. Figure 6.16 shows the contrast lag diagram for the two surfaces. Darker shading indicates lower contrast. Contours at 500 unit intervals are included to emphasise surface variation.

Figure 6.16Figure 6.16 - Contrast lag diagrams for fractal surfaces with (a) 1/8th and (b) all Fourier coefficients.

Broadly, these two images appear similar, as they suffer from the same problem as Moran's I lag diagrams as representations of high frequency texture. They do however demonstrate the high degree of scale dependence of the measure, which ranges from order 10 at the lowest lags to order 10,000 at the largest lags. To investigate differences in texture the central portion ofthe two images was enlarged and the percentage difference in contrast between the two was calculated (something that is easily accomplished in a raster GIS environment). Figure 6.17a shows the scaled difference between the contrast values of the fractal surface and its smoother equivalent with 1/8th of the Fourier coefficients. Figure 6.17b shows a similar difference map between the fractal and is equivalent surfaces with 1/4 of the Fourier coefficients. Both images show lags of up to 12 raster cells (17 along diagonals).

Figure 6.17Figure 6.17 - Percentage difference in contrast measures between (a) fractal and fractal with 1/8th coefficients, and (b) fractal and fractal with 1/4 coefficients.

The difference map shows that proportionally, the greatest difference in texture occurs at the highest frequency (as would be expected). More significantly, a spatial resolution can be identified at which this difference becomes significant. A 10% difference occurs at a lag of approximately 4 raster units in Figure 6.17a and 2 raster units in Figure 6.17b. This corresponds with the doubling of Fourier coefficients between the two surfaces.

A similar proportional difference map was produced for the Inverse Distance Moment measure and is shown in Figure 6.18. The spatial distribution of difference appears contrasts with Figure 6.17, and does not appear to relate to scale in an obvious way. While demonstrating that the Inverse Distance Moment measures different texture properties to Contrast, it is not clear how its meaning should be interpreted.

Figure 6.18Figure 6.18 - Scaled differences between the Inverse Distance Moment lag maps of fractal and fractal with 1/8th transformed coefficients.

6.5 Terrain Characterisation.

This section describes how the multi-scale characterisation techniques developed can applied to DEMs representing 'real' terrain. The aim is not to provide a detailed geomorphological analysis, but to demonstrate how, with visual examples, these tools could be used as a mechanism for such analysis. In particular, it is the aim of this section to demonstrate that a multi-scaled characterisation reveals more (useful) information that traditional raster-based geomorphometric analysis.

Figure 6.19 shows the cross-sectional curvature derived from the Lake District DEM at kernel scales of 150m, 450m, 850m and 1.25km. The area represented in the figure is approximately 8 x 8 km with the Wasdale valley in the southwest and the Ennerdale valley running northwest - southeast along the north of the image.

It is clear that the pattern of local curvature varies considerably with scale. At the finest (150m) scale, a rather fragmented network of ridges and channels is revealed. Smaller, well dissected upland valley systems are relatively well defined, while the major valleys are notdelineated. Much of the surface is represented as poorly autocorrelated minor concavities and convexities. At a coarser scale, there is far greater spatial autocorrelation in measurement. Figure 6.19d for example shows an almost entirely connected valley network. This suggests that either the true surface varies in roughness with scale (ie not self-similar), or that there is a (random) noise or error component at the highest frequency modelled by the DEM.

Figure 6.19Figure 6.19 - Cross-sectional curvature maps of selected portion of the Lake District DEM. Kernel size set to (a) 3 x 3 (150m); (b) 9 x 9 (450m); (c) 17 x 17 (850m); and (d) 25 x 25 (1.25km).

The pattern of changing curvature with scale is not symmetrical for concavity and convexity measurements. The 'characteristic scale' beyond which networks are well connected appears to be finer for convexities than for concavities. At the 450m scale (Figure6.19b), ridges appearwell defined while channels systems are represented discontinuously. At the 1.25km scale (Figure 6.19d) channels are still broadly linear features, while surface convexities are in many cases 'spread' over peak features. This asymmetry in pattern is itself a useful diagnostic feature. It is indicative of a glaciated mountain environment previously dominated by erosional processes around ridges and depositional processes along major valley systems.

The animated sequence of changing cross-sectional curvature with scale from which the four images in Figure 6.19 were taken, demonstrates the inadequacy of traditional (single scaled) raster based measurements. A large valley systems such as Wasdale (southwest of the image) is clearly of some geomorphological importance, yet would not be 'detected' by raster processing at the DEM resolution of 50m. It would seem an unnecessary constraint to consider geomorphological surface variation at one scale alone if a single DEM can reveal pattern over a range of scales.

Figure 6.20Figure 6.20 - Cross-sectional curvature for Dartmoor DEM. (a) 9x9 kernel; (b) 21 x 21 kernel; (c) 33 x 33 kernel; (d) 45 x 35 kernel.

Figure 6.20 shows a similar set of cross-sectional curvature values for the Dartmoor region. In this case the kernel size ranges up to 45 x 45 cells (2.3km) over an area of 14km x 14km. As with the Lake District example, the network pattern varies considerably with scale. Again, as with the Lake District, the change in curvature with scale appears asymmetric, but this time, channels are more well defined at finer scales than ridges. It also appears that there is greater spatial variation in convexity with scale than with concavity. This implies the existence of well defined 'V-shaped' channels that are expressed over a range of scales.

The degree to which surface measurements vary with scale was investigated interactively using d.param.scale (see Appendix). An indication of scale dependency in morphometric feature classification is shown in Figure 6.21.

Figure 6.21Figure 6.21 - Dartmoor modal feature classification. (a) morphometric features (pit -black, channel-blue, pass-green, ridge-yellow, peak-red, multimodal-grey); (b) classification entropy (dark = high entropy).

Figure 6.21a displays the most frequent feature classification of each DEM cell based on kernel sizes ranging from 9 x 9 to 45 x 45 cells. Where a cell has more than one mode, it is coloured grey. The relationship with relief is emphasised by combing this classification withshaded relief using an hue-intensity mapping. Figure 6.21b shows the confidence with which classification can be made by showing the scaled entropy. The lighter the image, the greater the consistency of feature classification (see section 5.3). These may be combined in a single image using an hue-intensity mapping (Figure 6.22)

Figure 6.22Figure 6.22 - Dartmoor feature classification certainty (see text).

Modal classifications are not all expressed at the same scale. That is, the size of each feature is not consistent over the surface. Features become progressively less well defined towards their edges, as the influence of a wider neighbourhood over a central cell's classification becomes increasingly dominant.

Figure 6.23 shows the scale dependency of feature classification for the Lake District DEM at four different kernel sizes. The slope and curvature tolerance values selected were 20 degrees and 2 respectively.

Figure 6.23Figure 6.23 - Feature classification of the Lake District DEM using kernel sizes of (a) 5 x 5; (b) 11 x 11; (c)17 x 17; (d) 25 x 25.

The certainty of feature classification is shown in Figures 6.24 and 6.25 using the mode and entropy method discussed in section 5.3.

Figure 6.24Figure 6.24 - Lake District feature classification (a) mode, and (b) entropy.

Figure 6.25Figure 6.25 - Lake District feature classification certainty.

One final example of feature classification is provided by the Peak District DEM shown in Figure 6.26. In this case a 15 x 15 cell kernel was passed over the 780 x 800 cell DEM. The distribution of peaks in particular is not uniform of the entire surface, with a higher density of peaks seen to the west and south of the image. This corresponds to the Carboniferous Limestone outcrop which contrasts with the flat topped moorland of Millstone Grit to the north and east. At this scale, the topological relationships between surface features may be of importance. Figure 6.27 shows the Delaunay triangulation of the thinned point features shown in Figure 6.26 (pits, passes and peaks) (see r.feat.thin and v.delaunay in the Appendix).

Figure 6.26Figure 6.26 - Peak District surface features.

Figure 6.27Figure 6.27 - Peak District feature triangulation.