Christopher D. Lloyd and Peter M. Atkinson
Department of Geography, University of Southampton, Highfield,
Southampton, SO17 1BJ, UK
Email: cdl195@soton.ac.uk
Information on the scale or frequency of spatial variation in properties such as land-form is of value in a wide variety of contexts including classification of land-form types and as an input to environmental modelling applications. This paper utilises this information to demonstrate how producers of digital terrain data sets may ascertain the best approach to employ and the nature and configuration of data that would be required to fulfil a particular user's requirements in terms of information and accuracy. The approach presented is applicable whether data are sampled on the ground or by remote sensing. The research centres around an examination of Ordnance Survey(R) Land-Form PROFILE (TM) contour data and the particular focus is the potential of the application of geostatistics to the improvement of the National Height Dataset (NHD). The research illustrates that to apply geostatistics it was necessary to ensure homogeneity of spatial variation across the region in concern. That is, it was necessary to classify spatial variation. Once classification of spatial variation was achieved it was possible to quantify the variation (using the variogram) and identify dominant forms of spatial variation for particular regions to aid the design of optimal sampling strategies. The results demonstrate that significant gains in efficiency can be obtained by adapting (i) the geostatistical approach and (ii) classifying the spatial variation prior to the application of geostatistics.
Key words: variogram, scale, sampling, kriging variance, DTM.
Geostatistics are of widespread potential value for sampling elevation data and the post-processing and analysis of Digital Terrain Models (DTMs). The variogram and related structure functions may be used to characterise spatial variation in landform for classification, estimation, simulation and the design of optimal sampling strategies. As a tool for the characterisation of landform the variogram complements elevation derivatives such as slope and aspect as well as other statistics. This paper examines Ordnance Survey Land-Form PROFILE contour data, the particular focus is the assessment of the potential of the application of geostatistics to the improvement of the National Height Dataset (NHD).
Several researchers have explored the use of geostatistics in the characteristation of land-form (for example, Mulla 1988; Bian and Walsh 1993). In particular, the variogram, a central tool of geostatistics, has been assessed for the design of sampling strategies. Balce (1987) discussed several techiques for determining optimal sampling intervals. In one approach the slope of a linear model fitted to a variogram with both axes logged was used to estimate iteratively the density of samples required to resolve the spatial variation represented by the variogram. The variogram is used with the same aim, but through different means, in this paper.
Geostatistics provides a means by which the measured spatial variation in digital elevation data may be used to ascertain the sample spacing necessary to achieve a particular accuracy. In short, using geostatistics it is possible to maximise the efficiency of the sampling strategy to achieve the requirements of individual users and, therefore, define the producer's requirements.
A major problem with using the variogram to design optimal sampling strategies as described above is that the variogram must be assumed representative of the entire region of interest. This assumption is often made necessary by a relative sparcity of data. In practice, this assumption may be unjustified: the variogram is often non-stationary meaning that it is not constant across geographical space. Therefore, a conventional geostatistical approach to the improvement of the National Height Dataset would be of limited value. To improve upon existing Ordnance Survey products using geostatistics it is necessary to identify regions where the measured spatial variation is homogenous. Once this has been achieved a geostatistical approach has the potential to make better use of existing data and to inform sampling strategies to complement the NHD. In this paper, the objective is to classify the region of interest into separate sub-regions which may be considered stationary in terms of the variogram prior to the design of optimal sampling strategies for those regions.
Geostatistics are based around the principle of spatial dependence: the likelihood that observations close together in space will be more alike than those further apart. This property may be represented using tools such as the variogram. There are many detailed introductory texts (for example, Journel and Huijbregts 1978; Isaaks and Srivastava 1989; Goovaerts 1997) and only a brief outline of the basics will be presented here.
The basis of geostatistics is the random function (RF) model. In geostatistics, spatially referenced variables are viewed as the outcome of a RF model. Thus, a set of individual realisations is termed a regionalised variable (ReV) and ReV theory underpins geostatistics.
An essential requirement of any geostatistical analysis is that the variogram (or related function) represents spatial variation sufficiently well. That is, the variogram must be stationary. Several forms of stationarity have been identifed by geostatisticians (Journel and Huijbregts 1978; Myers 1989). Stationarity may be divided into three major classes: strict, second-order and intrinsic stationarity. A RF is called strict stationary if its multivariate cumulative distribution function (cdf) is invariant under translation (Journel and Huijbregts 1978; Deutsch and Journel 1992). Second-order stationarity requires that the expected value of a variable should not depend on its location and that the covariance exists and also that it should depend only on the lag, h (the vector by which data locations are separated). Intrinsic stationarity is weaker in that it requires only that the variance of the increments should exist and, again, that it should depend only on the lag.
The variogram (also known as the semivariogram) is the most commonly used measure of spatial variation in geostatistics. It relates the semivariance, half the expected squared difference between paired data values {Z(x),Z(x+h)}, to the lag h, by which data locations are separated:
The experimental or sample variogram is computed as the average semivariance for several lags:
The subscript v refers to the support of the observation: the size, geometry and orientation of the area over which measurements are made.
Spatial scale is an important consideration in the analysis of any form of spatial data, and digital terrain data are no exception. Tools such as the variogram which quantify spatial variation and spatial scale have, therefore, a direct role in the analysis and use of digital terrain data.
The variogram may be fitted with a mathematical model, which is usually one of the so-called `authorised models' (McBratney and Webster 1986). The coefficients of the fitted model may be used to assign optimal weights for mapping using kriging. Additionally, the variogram may be used to relate the sampling configuration to the estimation error and it is this use which forms the basis of the present research.
If there is a limit to variation and the variogram model reaches a point at which the semivariances cease to increase with lag then the variogram model is termed bounded and contains several coefficients. The nugget variance, c0 (5.325 in figure 1), is the intercept of the model with the y-axis. The n structured components, c1, ...,cn (the two structured components in figure 1 are 369.061 and 440.641), represent the variation accounted for by each component of the model. The nugget variance and the structured components together equal the sill, equal to the a priori variance. In figure 1 the sill equals 5.325 + 369.061 + 440.641 = 815.027. The point on the x-axis at which the sill is reached is called the range, a (892.948 in figure 1). The range is perhaps the most useful coefficient of the variogram and an approach to the approximation of the range for the purposes of classification is discussed in section 2.6.
Fig. 1. Variogram with a nugget variance, a Gaussian and a spherical component.
Kriging is known by the acronym BLUE: the best linear unbiased estimator. Its main advantage over other interpolation techniques is that it assigns weights to observations on the basis of the form of the variogram (or related function) rather than assigning some arbitrary weighting function (a good introduction to kriging is given in Oliver and Webster 1990). The most widely used variant of kriging is ordinary kriging (OK) (Deutsch and Journel 1992). OK estimates are weighted averages of the available n observations, {z(xi), i=1,2,3,...,n}:
To ensure that the estimate is unbiased the weights (where the weights are lambda) sum to one:
The estimation variance is the expected value of the squared difference between the estimate and the true value:
This may be expressed as:
Where g(V,xi) (g is used here to mean gamma) is the average semivariance between the block to be estimated and the ith sample point, g(xi,xj) is the semivariance between all pairs of data locations and g(V,V) is the within-block variance (Webster and Oliver 1990).
The estimation variance is mimimised with the addition of a Lagrange parameter (Psi):
with the constraint that the weights sum to one.
A by-product of kriging is the kriging variance. The kriging variance is a measure of confidence in estimates and is a function of the form of the variogram (or other related function) and the sampling configuration only; it is independent of the data values (Goovaerts 1997).
The kriging variance is given as:
For a point support the average g(xi, V) becomes g(xi,x0) and the within-block variance, g(V,V) is zero.
The kriging variance provides a means by which to assess the suitability of a particular data configuration to meet the accuracy requirements of the user in a particular situation (Atkinson 1996). Where the aim is the design of an optimal sampling strategy it is obviously necessary first to obtain a representative variogram. Observations are frequently obtained along a transect to provide data for this purpose.
Once a representative variogram has been obtained and a model has been fitted the coefficients of the model can be used to compute the maximum kriging variance for different sample spacings (McBratney et al. 1981). It is then possible to ascertain the largest sample spacing that could be used to achieve the desired kriging variance. This approach is relevant whether measurements are obtained by ground survey or by remote sensing (for the latter see Atkinson and Curran 1995; Atkinson 1997)
The range was approximated through a process consisting of three steps: (i) a third-order polyonomial was fitted to the variogram, (ii) the first derivative of the polynomial was obtained and (iii) the smallest of two roots was extracted from the derivative function and this was used to approximate the range (with code published by Press et al. 1986).
The study site is a region of Dorset, South-West England (Ordnance Survey tile ST92SW). The data set comprises a five kilometre square section of the Ordnance Survey Land-Form PROFILE contour and spot height data set. The region is bounded by the National Grid References (NGRs) South-West: 120000, 390000; North-East: 125000, 395000. Land-Form PROFILE data have contours at a 5 m vertical interval and 10 m in some mountainous regions.
The data set was converted to ASCII x,y,z format to enable its analysis with existing geostatistical software. The data set analysed comprised contour nodes and spot heights. The nature of the data means that spatial variation is smoothed but the dominant form of spatial variation appeared to be captured well by the variogram.
Fig. 2. Land-Form PROFILE
(TM) contours and spot heights for region ST92SW.
Based on Ordnance Survey Land-Form PROFILE data with the permission of The
Controller of Her Majesty's Stationery Office, (C) Crown copyright,
Southampton University, License no. ED281247.
Fig. 3. Omnidirectional variogram of region ST92SW.
To assess the problem of change in spatial variation across the region of analysis variograms were computed for twenty five 1 km square sub-regions. The variograms (figure 4) clearly illustrate the heterogenous nature of spatial variation across the region. Since the global variogram masks this variation a reasonable goal is to ascertain regions for which the variogram is representative.
Fig. 4. Omnidirectional variograms of 1 km sub-regions. Fitted models are
indicated by solid lines.
The development of an approach for the classification of spatial variation using the variogram is an important step in using the variogram to design optimal sampling strategies. An approach was developed by which the variogram was computed for a moving window (500 m by 500 m, incrementing by 50 m) and various coefficients, including an approximation of the range (computed as discussed in section 2.6), were obtained. Variograms were computed using code based on that supplied in GSLIB.
Ranges above the cut-off (which were a very small proportion of the total) were extracted and replaced with weighted averages of their neighbours. A 5 by 5 median filter was then used to smooth the ranges to reduce the effect of local noise. The maximum slope of the locally fitted first-order polynomial was used in conjunction with the range to classify spatial variation because the two variables describe the dominant scales of spatial variation within each window. A K-means clustering algorithm was used for this purpose. Maps of the maximum slope and the approximated range (smoothed with 5 by 5 median filter) are shown in figures 5 and 6 respectively.
Fig. 5. Map of maximum slope of local first-order polynomial, scale is in
0.1m/m. The anomalous areas contrasting with the otherwise smooth variation
are missing values where a polynomial could not be fitted due to insufficient
numbers of data or other factors. The missing values can be replaced with
weighted averages.
Fig. 6. Map of approximated ranges, scale is in metres.
The map of classes for which variograms were computed is given in figure 7.
Fig. 7. Map of classes. Class 1 = white; class 2 = light grey; class 3 =
dark grey; class 4 = black.
Fig. 8. Class-based variograms. Class 1 to 4, top to bottom, left to right.
| Nugget | Gaussian | ||
| Class | c0 | a | c |
| 1 | 0.08 | 87.018 | 79.211 |
| 2 | 0.38 | 103.331 | 379.994 |
| 3 | 0.42 | 130.895 | 423.886 |
| 4 | 0.58 | 124.707 | 576.359 |
The coefficients in table 1 indicate that both the range and the structured component vary between different classes. The nugget variances are approximately zero in all cases due to the nature of the data (that is, the contours represent smoothed spatial variation).
The maximum kriging variances for sample spacings from 2 to 40 m were computed for both the global and class-based variogram models. The relevant figures are given in table 2. The increases in maximum kriging variance with increase in sample spacing are more clearly illustrated in figures 9 and 10, which are plots for the global and the class-based variogram model coefficients respectively.
| Spacing (m) | Global | C1 | C2 | C3 | C4 |
| 2 | 0.600 | 0.100 | 0.475 | 0.525 | 0.725 |
| 4 | 0.600 | 0.100 | 0.475 | 0.525 | 0.725 |
| 6 | 0.601 | 0.100 | 0.477 | 0.525 | 0.726 |
| 8 | 0.604 | 0.100 | 0.481 | 0.527 | 0.729 |
| 10 | 0.610 | 0.101 | 0.491 | 0.532 | 0.736 |
| 12 | 0.622 | 0.103 | 0.509 | 0.539 | 0.749 |
| 14 | 0.641 | 0.106 | 0.538 | 0.552 | 0.770 |
| 16 | 0.670 | 0.111 | 0.582 | 0.571 | 0.802 |
| 18 | 0.712 | 0.118 | 0.646 | 0.599 | 0.848 |
| 20 | 0.771 | 0.128 | 0.734 | 0.638 | 0.912 |
| 22 | 0.849 | 0.140 | 0.852 | 0.690 | 0.997 |
| 24 | 0.951 | 0.157 | 1.006 | 0.758 | 1.109 |
| 26 | 1.082 | 0.179 | 1.201 | 0.845 | 1.252 |
| 28 | 1.245 | 0.205 | 1.444 | 0.953 | 1.430 |
| 30 | 1.445 | 0.238 | 1.743 | 1.087 | 1.649 |
| 32 | 1.687 | 0.278 | 2.102 | 1.249 | 1.914 |
| 34 | 1.977 | 0.325 | 2.530 | 1.442 | 2.231 |
| 36 | 2.319 | 0.381 | 3.033 | 1.671 | 2.606 |
| 38 | 2.718 | 0.447 | 3.619 | 1.939 | 3.044 |
| 40 | 3.181 | 0.522 | 4.294 | 2.249 | 3.551 |
It should be noted that the nugget variance of the variogram model contributes directly to the kriging variance. Since the nugget variances fitted to the variograms were, in effect, arbitrary, the maximum values of kriging variance implied for a sample spacing of zero should be considered approximately zero in these cases. With this in mind, figure 9 indicates that the maximum kriging variance up to a sample spacing of perhaps 14 m is relatively constant (around 0.6 m2). This suggests that there would no advantage in sampling on a grid finer than 14 m as the accuracy would be almost the same as with a 14 m grid. However, as observed above, the global variogram was computed for data which exhibited heterogenous spatial variation. The estimates of the required sample spacing indicated by the global variogram model coefficients, therefore, may be inappropriate locally.
Fig. 9. Maximum kriging variance against sample spacing for global variogram
model coefficients.
The same process was followed using the coefficients of the models fitted to the four class-based variograms. The four graphs in figure 10 exhibit similar forms to the graph in figure 9, but obvious differences are apparent. The graph for class 1 is, relatively, rather flat in appearance. The maximum kriging variance is relatively constant up to a sample spacing of about 18 m for which the maximum kriging variance is only 0.018 m2 greater than for a sample spacing of 2 m. The cut-off for the graph of class 2 is clearly smaller; the maximum kriging variance is constant to a sample spacing of about 10 to 12 m.
If the stated accuracy was 0.75 m2 the graph based on the global variogram model coefficents would imply a necessary sample spacing of between 18 and 20 m. For the class 1 variogram model coefficents the corresponding sample spacing would be >40m. For classes 2, 3 and 4 the necessary sample spacings would be estimated as (2) approximately 20 m, (3) between 22 and 24 m and (4) approximately 12 m. Allowing for the fitted nugget variances the differences between the classes are clear. The global variogram model coefficients would lead to over-sampling by 100 % (in terms of the sample spacing) for the region covered by class 1 whilst the regions of class 2 and 3 would be slightly under-sampled. Class 4 would, in contrast, by under-sampled by about 80 %.
Fig. 10. Maximum kriging variance against sample spacing for classified regions.
Class 1 to 4, top to bottom, left to right.
It is clear that the kriging variance is subject to specific limitations. It serves as a general guide to support sampling but is heavily dependent on how well the variogram represents spatial variation. The classification or segmentation of spatial variation are approaches that may help to ensure that the variogram is representative of spatial variation and thereby facilitate improvements in the use of the technique. However, it is necessary both to test the efficiency of the kriging variance and to utlise other techniques to obtain a guide to the sampling strategy required (whether in terms of irregularly distributed samples or the spatial resolution of DTMs derived from remote sensing).
The fact that the kriging variance does not take into account the data values is a particular weakness. Where neighbouring data have large differences the potential error is greater than when neighbouring data are very similar. However, the kriging variance would treat two identical data configurations in the same way irrespective of the data values. Clearly, to solve this problem measures of uncertainty which are data dependent but which still take into account the form of spatial variation and the data configuration would need to be utilised. Goovaerts (1997) outlines an approach to ascertaining uncertainty through the use of local probability distributions. Future work will assess local uncertainty through the conditional variance and other measures.
The method presented may serve as a guide to the acquisition of data whether measurements are obtained by ground survey or by remote sensing The modelling of spatial variation, if such variation is classified or segmented effectively, can provide a perspective which may support existing means of assessing the best approach to fulfilling the requirements of particular users of digital terrain data. It is likely that a segmentation (region growing) approach will provide superior results to those obtained using an aspatial classification algorithm. Future work will focus on the automated segmentation of spatial variation. Finally, it was noted that the use of the kriging variance is problematical as it does not take the data values into account. Therefore, another primary aim of future work must be the examination of other measures of local uncertainty to better use the variogram (or related functions) in the maximisation of sampling efficiency.
Atkinson, P. M. Optimal sampling strategies for raster-based geographical information systems. Global Ecology and Biogeography Letters, 5 (1996), 271-280.
Atkinson, P. M. Selecting the spatial resolution of airborne MSS imagery for small-scale agricultural mapping. International Journal of Remote Sensing, 18 (1997), 1903-1917.
Atkinson, P. M. and Curran, P. J. Defining an optimal size of support for remote sensing investigations. IEEE Transactions on Geoscience and Remote Sensing, 33 (1995), 768-776.
Balce, A. E. Determination of optimum sampling interval in grid digital elevation models (DEM) data acquisition. Photogrammetric Engineering and Remote Sensing, 53 (1987), 323-330.
Bian, L. and Walsh, S. J. Scale dependencies of vegetation and topography in a mountainous environment of Montana. Professional Geographer, 45 (1993), 1-11.
Deutsch, C. and Journel, A. G. GSLIB: Geostatistical software library and User's Guide. New York: Oxford University Press (1992).
Goovaerts, P. Geostatistics for Natural Resources Evaluation. New York: Oxford University Press (1997).
Herzfeld, U. C. and Higginson, C. A. Automated geostatistical seafloor classification - principles, parameters, feature vectors, and discrimination criteria. Computers and Geosciences, 22 (1996), 35-52.
Isaaks, E. H. and Srivastava, R. M. An Introduction to Applied Geostatistics. New York: Oxford University Press (1989).
Journel, A. G. and Huijbregts, C. J. Mining Geostatistics. London: Academic Press (1978).
McBratney, A. B. and Webster, R. The design of optimal sampling schemes for local estimation and mapping of regionalised variables. II. Program and examples. Computers and Geosciences, 7 (1981), 335-365.
McBratney, A. B. and Webster, R. Choosing functions for semi-variograms of soil properties and fitting them to sampling estimates. Journal of Soil Science, 37 (1986), 617-639.
McBratney, A. B., Webster, R. and Burgess, T. M. The design of optimal sampling schemes for local estimation and mapping of regionalised variables. I. Theory and method. Computers and Geosciences, 7 (1981), 331-334.
Miller, A. J. Algorithm AS 274: Least squares routines to supplement those of Gentleman. Applied Statistics, 41 (1992), 458-478.
Miranda, F. P., MacDonald, J. A. and Carr, J. R. Application of the semivariogram texture classifier (STC) for vegetation discrimination using SIR-B data of Borneo. International Journal of Remote Sensing, 13 (1992), 2349-2354.
Mulla, D. J. Using geostatistics and spectral analysis to study spatial patterns in the topography of southeastern Washington state, U.S.A. Earth Surface Processes and Landforms, 13 (1988), 389-405.
Myers, D. E. To be or not to be ... stationary? That is the question. Mathematical Geology, 21 (1989), 347-362.
Oliver, M. A. and Webster, R. Kriging: a method of interpolation for geographical information systems. International Journal of Geographical Information Systems, 4 (1990), 313-332.
Pebesma, E. J. Gstat User's Manual. Amsterdam: Netherlands Centre for Geo-Ecological Research (ICG), Faculty of Environmental Sciences, University of Amsterdam (1998). http://www.frw.uva.nl/~pebesma/gstat/.
Press, W. H., Flannery, B. P., Teukolsky, S. A. and Vetterling, W. T. Numerical Recipies: The Art of Scientific Computing. Cambridge: Cambridge University Press (1986).
Webster, R. and Oliver, M. A. Statistical Methods in Soil and Land Resource Survey. Oxford: Oxford University Press (1990).