Return to GeoComputation 99 Index

Designing optimal sampling configurations with ordinary and indicator kriging

Christopher D. Lloyd and Peter M. Atkinson
University of Southampton, Department of Geography, Highfield, Southampton SO17 1BJ, United Kingdom
E-mail: cdl195@soton.ac.uk

Abstract

The objective of this paper is to examine the applicability of two geostatistical approaches, ordinary kriging (OK) and ordinary indicator kriging (IK), to the design of optimal sampling strategies. This paper uses the OK variance and the conditional variance of the conditional cumulative distribution function (ccdf) derived through IK to assess local uncertainty in estimates. The mean OK variance and IK variance for a given grid spacing, using data sampled from a remotely sensed DTM, were used to ascertain the sample spacing required to achieve a particular accuracy. The estimates of optimal sample spacing were assessed with reference to the complete DTM. Judgement on the success of the two approaches was made on the basis of the correlation between the OK variance and the errors of OK estimation and between the IK variance and the errors of IK estimation, and the form of the relationships were discussed. The IK variance was found to represent the estimation errors more accurately than the OK variance, although OK estimates of elevation values were more accurate than those for IK. The IK is significantly more costly to implement than OK in terms of expenditure of time and effort and the implementation of the technique was demonstrated to be problematic in the presence of a low frequency trend. Despite these limitations IK was recommended for the design of optimal sampling strategies where the estimation of accuracy of a specified degree is particularly important and where the analysis relates to an area over which the spatial variability may be considered similar. Comparison and integration with a segmentation approach was suggested for future work, and the results of a related approach presented, to improve the implementation of IK in the context of digital elevation data.

1. Introduction

Optimal strategies for sampling digital elevation data have been designed using a variety of techniques including, for example, the tools of geostatistics (Lloyd and Atkinson, 1998) and fractal analysis (Balce, 1987). Geostatistics have been used for the analysis of landform in several contexts, such as the characterisation and classification of landforms (Mulla 1988; Bian and Walsh 1993) including seafloor landform (Herzfeld and Higginson 1996). Elevation has also been used as a co-variate for kriging (for example, Leenaers and others, 1990).

Kriging is, in basis, a form of weighted average estimator. The weights are assigned on the basis on the form of a model fitted to a function, such as the variogram, which represents spatial variability in the property of interest. The most frequently used form of kriging is ordinary kriging (OK). A by-product of OK is the OK variance, a measure of confidence in estimates. Previously, to determine an acceptable sample grid spacing, investigators have plotted the maximum OK variance for a range of sample spacings and used the plot to select a sample spacing that achieved a given precision of estimation (McBratney and others, 1981; McBratney and Webster, 1981; Lloyd and Atkinson, 1998).

Geostatistical approaches, such as kriging, rely on measures of spatial variability such as the variogram. Where the estimates and statistics derived are sensitive to small changes in the variogram, or other function, the use of such approaches may be problematic, as spatial variability of landform is rarely similar over even small areas. The results obtained may, therefore, be biased when the variogram is not representative of the local area. One solution to this problem is to determine sub-regions of a data set within which the spatial variability may be considered similar. Measures of spatial variability may then be obtained for each segment. This kind of approach has yielded reasonable results (Lloyd and Atkinson, 1998) where the OK variance has been used as a measure of confidence in estimates.

One major limitation of the OK variance is that it is not conditional on the data values; therefore, unless the spatial variability is very similar across the region in concern, the results may grossly over-estimate or under-estimate confidence in the estimates. An alternative geostatistical approach, with the potential to overcome some of the problems considered, is IK. The IK alternative to the OK variance is the IK (conditional) variance. Since the conditional variance of the conditional cumulative distribution function (ccdf) estimated using IK is conditional on the data values, the problem of data independence of the OK variance is overcome. In this paper, the mean OK variance and IK conditional variance are used in the same way as the maximum OK variance has been used to relate estimation error to sample spacing.

2. Geostatistics

The most widely used tools of geostatistics have been presented in a variety of good introductions (for example, Journal and Huijbregts, 1978; Isaaks and Srivastava, 1989; Goovaerts, 1997) and only a summary of the tools used in this paper will be given here.

The principal tool of most geostatistical analyses is the variogram, a function that relates half the average squared difference between paired data values to the distance (and direction, where anisotropy is considered) by which they are separated. A mathematical model may be fitted to the variogram and the coefficients of the model may be used to assign optimal weights for interpolation using kriging. As stated previously, OK is the most widely used variant of kriging. A by-product of OK is the OK variance (or standard error). It is a function of (1) the form of spatial variability of the data (modelled, for example, by the variogram) and (2) the spatial configuration of the observations. The OK variance serves as a measure of confidence in estimates although, since it is independent of the data values, its practical value is frequently severely limited unless the data are segmented into sub-sets within which the spatial variability is similar. For a given sampling configuration, the OK variance will be the same irrespective of the data values locally. If data are measured on a regular grid, the maximum OK variance will be identical for a given location; however, much of the data values vary locally (Goovaerts, 1997). A more suitable approach, especially where the variogram is not stationary, would be to use a measure that was conditional on the data values.

3. The indicator approach

The utility of conventional variograms and ordinary kriging for geostatistical analysis is limited by particular assumptions. The principal problems concern (i) normality (although the techniques are fairly robust) and (ii) the independence of the estimation variance on data values. The indicator approach is one means of overcoming both of these limitations. The indicator approach involves the designation of several thresholds, for example the deciles of the sample cumulative distribution (Goovaerts, 1997) that are used to transform the data to a set of binary variables. The indicator transform (as defined in Deutsch and Journel, 1998) I(x) (at location x) with cut-off cut_k for datum value z(x) may be given as:

(1)

Variograms are then computed for the k classes using the transformed data.

3.1 Indicator kriging

With IK, a least-squares estimate of the ccdf is made at each cut-off. The ccdf is constructed by putting together the IK estimates for each cut-off. The ccdf signifies a probabilistic model for the uncertainty around the unknown value z(x) (Deutsch and Journel, 1998). The "E-type" estimate is the mean of the ccdf, which is compared to OK estimates in this study. The IK (conditional) variance, used in this paper to measure uncertainty in estimates, is a measure of dispersion of the ccdf around its’ mean (Goovaerts, 1997).

4. Study area and data

The analysis was based on a subset of a digital terrain model (DTM) derived as part of the LANDMASS project (Ridley and others, 1997). The DTM represents an area in the Lake District, Northern England and falls primarily in the British Ordnance Survey^® tile NY11NE. The spatial resolution of the DTM was 5 m, and the cells were treated as discrete points (the sample spacing being, equal to the spatial resolution). The geographical limits of the subset DTM (Figure 1) defined in British National Grid References were: North 517742.68; South 515902.68; West: 315953.79; East 317353.79. The edges of the DTM were not regular but the maximum number of rows was 368 pixels and the maximum number of columns was 280 pixels.

Figure 1. LANDMASS DTM (region NY11NE) with values estimated, using OK, to a 5-m grid. The Northwest corner of the region is an area of water.

5. Methods and approaches

5.1 Variograms

In this paper, omnidirectional variograms only were computed. While there is clear directional variability in the DTM, this initial comparison of the OK and IK approaches was restricted to maintain simplicity, particularly in relation to IK, and to serve as a basis for future work.

An omnidirectional variogram was estimated for the entire region (Figure 1). All standard experimental variograms were estimated and all models were fitted in Gstat (Pebesma and Wesseling, 1998). The shape of the landform was a gradually sloping valley. This was apparent from the form of the variogram: semivariance increased without limit, a common indication of the presence of a low-frequency trend (Starks and Fang, 1982). A variogram of residuals from a second-order polynomial, for which the range and sill were more distinct, was used instead and kriging was carried out using the residuals. A nugget and two Gaussian structures (the coefficients are given in Figure 2) were fitted to the variogram. It was found necessary to increase the nugget beyond that originally fitted (to about 1.4% of the sill) as the use of Gaussian model components resulted initially in negative kriging variances. Problems with the kriging matrix are frequently encountered when Gaussian models are used (Deutsch and Journel, 1998).

The OK and IK approaches were compared by referring to estimates of the residuals and measures of confidence in these estimates. The use of ordinary least squares fitting to estimate a trend directly has been criticised although, in practical applications, this approach often provides a good result (Kitanidis, 1993; Kitanidis, 1997). Furthermore, accuracy's of estimates obtained in work relating to this paper were generally higher where the data were detrended in this manner prior to kriging than when the raw data were used at all stages. It would normally be necessary to integrate the uncertainty associated with the polynomial trend with the kriging variances, although in this analysis the residuals were treated directly as the variables of interest.

Figure 2. Omnidirectional variogram of residuals from a second-order polynomial trend. The poor fit of the model at the origin of the variogram is due to the necessity of increasing the nugget variance to prevent negative OK variances.

5.2 Ordinary kriging

The ordinary kriging functionality of GSLIB (Deutsch and Journel, 1998) was used. As noted previously, kriging was carried out using residuals from a second-order polynomial trend. For OK and IK, a search radius of 100 m was used with a minimum of four and a maximum of 16 observations used for estimation. The sole exception was that a minimum of 8 observations were used for IK in the segment-based analysis presented at the end of the analysis section.

5.3 Indicator variograms

Indicator variograms were computed using the gamv program provided in GSLIB (Deutsch and Journel, 1998). The cut-offs were based on the deciles of the data. Since the data are located on a regular grid, no declustering algorithm was used. As for the experimental variogram, the indicator variograms were estimated for the residuals from a second-order polynomial trend. The cut-offs and models fitted to the indicator variograms are given in Table 1. The models were fitted using a combination of fitting by eye and the weighted least squares functionality of Gstat.

Table 1. Indicator variogram cut-offs for residuals from a second-order polynomial trend and fitted models.

5.4 Indicator kriging

Indicator kriging was carried out using the ik3d code of GSLIB (Deutsch and Journel, 1998). The "E-type" estimates and conditional variances were derived using the GSLIB postik routine. For the main part of the analysis, linear models were used to extrapolate the tails of the ccdfs, although the models and parameters were varied for the latter part of the analysis.

6. Analysis

The implementation of IK proved problematic in this analysis because of the presence of a low-frequency trend — the gradually sloping side of a valley. Although the low frequency variation was reduced using the residuals from a fitted second-order polynomial, this function was insufficient to remove all variation of low frequency. For OK this did not present a problem, but for IK, the end result was a series of patches of estimates with the same value, giving in some areas a poor representation of the topographic form. IK decomposes the variable into a set of binary classes, but in the presence of a clear trend, this may mean that large sub-regions of a data set may comprise only members of one class. When the local ccdf is estimated, the available observations may all belong to a single class with the result that estimates will be identical across quite large areas, creating flat patches. This can be seen in the map of IK variances in Figure 3. This problem does not affect the illustration of the IK conditional variance as a useful statistic. The problem does require resolving through some means, such as segmentation of spatial variability, into sub-sets for which the spatial variability is either similar, or effectively be made similar, by detrending with a low-order polynomial (that is, less than an order of two). A brief example of such an approach will be given after the main analysis.

6.1 Analysis of subset DTM

Initially, the OK and IK were carried out using the subset DTM. Following this, a segmentation-based approach was used.

Figure 3. Map of IK variances.

Mean and root mean square (RMS) errors and estimation standard errors (square root of estimation variances) are given for OK in Table 2. Summary statistics for the differences between the errors of OK estimates and the OK estimation standard errors are given in Table 3. It should be noted that direct comparison of the error and estimation standard errors does not imply that one is directly equivalent to the other but simply that the latter may serve as a guide to the former and that some relationship between the two may be expected. It is clear that, globally, the OK standard error does not represent the estimation errors well. For the 10-m sample spacing, the mean and RMS OK standard errors are between two and three times larger than the errors themselves, which would indicate expenditure on sampling effort much greater than that necessary. The difference between the mean and RMS OK standard errors and the mean and RMS errors increases from a grid spacing of 10 m to 25 m and then reduces slightly for grid spacings of 30 m and 35 m (Table 3). Table 3 indicates a tendency for OK to under-estimate the error (since the mean error minus estimation standard error is negative), though the largest differences are for over-estimation of the error.

Table 2. Mean OK estimation errors and variances (standard errors) for different sample spacings. All units are in metres.

Table 3. Summary statistics for error minus OK estimation variance. All units are in metres.

It was noted previously that implementation of IK in this analysis proved problematic. This is reflected in the values given in Tables 4 and 5. The mean and RMS errors are much larger than those for OK (Tables 2 and 3). In terms of the mean and RMS IK standard errors and the mean and RMS errors, IK has estimated more successfully the errors with, for example, a difference of only 0.87 m in the mean error and mean IK standard error. The spread of differences between errors and estimated errors is greater for IK than it was for OK, with a more symmetrical distribution of values. The minima and maxima for the errors in estimation minus the IK estimates indicate again that the presence of a low frequency trend has resulted in estimates that do not match the true values well in areas that constitute only members of a single indicator class.

Table 4. Mean IK estimation errors and variances (standard errors) for different sample spacings. All units are in metres.

Table 5. Summary statistics for error minus IK estimation variance. All units are in metres.

The correlation between the OK variance and the errors was, as expected, very small. The r² figures typically being of the order of 0.001. In contrast the r² figures for the IK variance against IK errors were: 10 m, 0.779; 15 m, 0.800; 20 m, 0.774; 25 m, 0.748; 30 m, 0.737; 35 m, 0.715; however, a large proportion of the estimates clearly deviated from a linear relation (as illustrated in Figure 4, for IK with a sample spacing of 10 m). The fairly large coefficients of determination were primarily due to the fact that most of the estimation standard errors tended to increase in value with increase in error, though in a curvilinear manner. The bulk of errors of less than 10 m were represented well by the IK variance, although errors between 10 m and 40 m were frequently over-estimated while error over 40 m tended to be under-estimated. The observed deviations in the plot were due, in part, to the problems noted above in relation to the use of IK in the presence of a low-frequency trend, which has not been accounted for sufficiently by detrending. The strips of values parallel with the x-axis are indicative of ‘platforms’ of identical values formed for the reason elucidated above; therefore, as for the estimates, there were large sub-regions for which the conditional variance was identical and, locally, did not represent estimation error well.

Figure 4. Error in IK estimates against IK estimation standard error, 10-m sample spacing.

6.2 Analysis of DTM segment

OK and IK also were carried out for a segment of the DTM with a 10-m sample spacing. The segment was identified using a combination of the first principal component of a locally-estimated experimental variogram and a region-growing segmentation algorithm (a similar approach was presented in Lloyd and Atkinson, 1998).

The models and parameters used to extrapolate the tails of the ccdfs were varied to identify the most effective combinations in terms of accuracy of estimation and the form of the relationship between the IK variance and the estimation errors. The experimental variogram of the segment data was fitted with a nugget and two spherical components. The mean error of the OK estimates was 2.56 m. The OK variances bore no relationship with the errors of estimation as the OK variances were generally one of only two values for the whole range of errors. That is, OK variances for estimates 5 m away from two observations in the same row and OK variances for estimates equidistant from four observations.

Table 6 gives the models and parameters for extrapolation of the ccdf tails with the mean, minimum and maximum errors of estimation. The w parameters (Goovaerts, 1997) were varied iteratively to reduce the mean estimation error reduced, but in doing so the estimation errors and the IK variance became less correlated for the models and parameters in table 6 that resulted in an mean error of 2.62 m and the estimation errors and the IK variance were uncorrelated (r² of 0.0009). In the other case (mean error of 5.05 m) the relationship was clearly curvilinear (with r² of 0.6451) as shown in Figure 5.

Table 6. Models and parameters for extrapolation of the ccdf tails and associated summary statistics of the estimation errors.

Figure 5. Error in IK estimates against IK estimation standard error, 10 m sample spacing, for a segment of the DTM.

Figure 5, as Figure 4 indicates, the majority of errors of IK estimation were less than 10 m, and the IK variance represented errors of less than 10 m with a general tendency to over-estimate errors of greater than 10 m.

7. Discussion

Clearly, a balance between accuracy of estimation and the usefulness of the IK variance should be sought. The reduction in accuracy of IK estimates, as observed above, was not a primary concern as the aim of the paper was to use the IK variance as a guide to accuracy; however, it is obvious that since the IK variance serves as a measure of confidence in IK estimates, those estimates must be acceptable. Given the correct choice of models and parameters for extrapolation of the tails of the ccdf, this requirement may be fulfilled. The IK variance does provide a useful guide to uncertainty that may be used to inform sampling strategies. The IK variance has been shown to represent error, both globally (in terms of mean error and mean IK variance), and locally (point-by-point comparison of error and IK variance). The IK variance may, therefore, be considered appropriate for the design of new sampling configurations (on a grid basis informed by mean IK variance), and for the sampling of data to improve existing data sets (by using mapped IK variance as demonstrated in Figure 3). The use of a segmentation-based approach as presented here improved the results and indicates a way forward.

The implementation of IK can be a lengthy task. In this analysis, the OK estimates were more accurate than those provided by IK (given the models and parameters used); however the IK variance has been shown to provide frequently an accurate estimate of the errors of IK estimation, unlike the OK variance. The implication of this analysis is that IK is a potentially useful tool where the aim is to design optimal configurations for sampling elevation data. The approach seems likely to be beneficial only where (i) there are sufficient indicator cut-offs; (ii) there is little low frequency spatial variability, or the trend can be effectively removed, either directly or in combination with segmentation, and (iii) the spatial variability is similar across the region of concern, which again may be resolved using segmentation.

8. Conclusions

The results presented in this paper have identified strengths and weaknesses of the OK and IK, and the paper has identified several key issues:

1. OK is easier to implement than IK and provides more accurate estimates than IK (i) where there are insufficient indicator classes and (ii) in the presence of a low frequency trend, that cannot be removed by detrending using a polynomial of order two or less.
2. In the absence of some form of segmentation of spatial variability, both the OK standard errors and the IK standard errors may be of limited value in the design of optimal sampling configurations.
3. The IK variance was much more correlated with the errors of the IK estimates than were the OK variances with the errors of the OK estimates.
4. The evidence indicates that the segmentation of spatial variability prior to the use of IK, may provide a powerful tool in the design of optimal sampling configurations.

In the context of the sampling of digital elevation data with the aid of IK, the principle research aim should perhaps be the integration of a segmentation-based approach (for example, Lloyd and Atkinson, 1998) with the use of IK. Such an approach was presented at the end of this analysis. This would hopefully combine to at least partially resolve the problems of (i) stationarity; (ii) the assumption of a Gaussian distribution, and also (iii) the data independence of the OK variance.

Acknowledgments

The authors would like to thank the University of Southampton, England and the Ordnance Survey^® in Southampton, England, for financial support. Dr. David Holland of the Ordnance Survey is thanked for providing the data on which this project was based.

References

Balce, A.E., 1987. Determination of optimum sampling interval in grid digital elevation models (DEM) data acquisition. Photogrammetric Engineering and Remote Sensing, 53, pp. 323-330.

Bian, L. and S.J. Walsh, 1993. Scale dependencies of vegetation and topography in a mountainous environment of Montana. Professional Geographer, 45, pp. 1-1.

Deutsch, C. and A.G. Journel, 1998. GSLIB: Geostatistical Software Library and User’s Guide. Second Edition. New York: Oxford University Press.

Goovaerts, P., 1997. Geostatistics for Natural Resources Evaluation. New York: Oxford University Press.

Herzfeld, U.C. and C.A. Higginson, 1996. Automated geostatistical seafloor classification — principles, parameters, feature vectors and discrimination criteria. Computers and Geosciences, 22, pp. 35-52.

Isaaks, E.H. and R.M. Srivastava, 1989. An Introduction to Applied Geostatistics. New York: Oxford University Press.

Journel, A.G. and C.J. Huijbregts, 1978. Mining Geostatistics. London: Academic Press.

Kitanidis, P.K., 1993. Generalised covariance functions in estimation. Mathematical Geology, 25, pp. 525-540.

Kitanidis, P.K., 1997. Introduction to Geostatistics: Applications to hydrogeology. Cambridge: Cambridge University Press.

Leenaers, H., J.P. Okx, and P.A. Burrough, 1990. Comparison of spatial prediction methods for mapping floodplain soil pollution. Catena, 17, pp. 535-550.

Lloyd, C.D. and P.M. Atkinson, 1998. Scale and the spatial structure of landform: optimising sampling strategies with geostatistics. In Proceedings of the 3rd International Conference on GeoComputation, University of Bristol, 17th–19th September 1998. GeoComputation CD-ROM, Manchester. http://divcom.otago.ac.nz/sirc/geoc98/15/gc_15.htm.

McBratney, A.B. and R. Webster, R., 1981. The design of optimal sampling schemes for local estimation and mapping of regionalised variables. II. Program and examples. Computers and Geosciences, 7, pp. 335-365.

McBratney, A.B., R. Webster, and T.M. Burgess, 1981. The design of optimal sampling schemes for local estimation and mapping of regionalised variables. I. Theory and method. Computers and Geosciences, 7, pp. 331-334.

Mulla, D.J., 1988. 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, pp. 389-405.

Pebesma, E.J. and C.G. Wesseling, 1998. Gstat, a program for geostatistical modelling, prediction and simulation. Computers and Geosciences, 24, pp. 17-31.

Ridley, H.M., P.M. Atkinson, P. Aplin, J.-P. Muller, and I. Dowman, 1997. Evaluating the potential of the forthcoming commercial U.S. high-resolution satellite sensor imagery at the Ordnance Survey. Photogrammetric Engineering and Remote Sensing, 63, pp. 997-1,005.

Starks, T.H. and J.H. Fang, 1982. The effect of drift on the experimental semivariogram. Mathematical Geology, 14, pp. 309-319.