Return to GeoComputation 99 Index

Performance comparison of geostatistical algorithms for incorporating elevation into the mapping of precipitation

Abstract

This paper presents three geostatistical algorithms for incorporating a digital elevation model into the spatial prediction of rainfall: simple kriging with varying local means, kriging with an external drift, and colocated cokriging. The techniques are illustrated using annual and monthly rainfall observations at 36 climatic stations in a 5,000 km² region of Portugal. Cross validation is used to compare the prediction performances of the three geostatistical interpolation algorithms with the straightforward linear regression of rainfall against elevation and three univariate techniques: Thiessen polygon, inverse square distance, and ordinary kriging.

Larger prediction errors are obtained for the two algorithms (inverse square distance, Thiessen polygon) that ignore both elevation and rainfall records at surrounding stations. The three multivariate geostatistical algorithms outperform other interpolators, in particular linear regression, which stresses the importance of accounting for spatially dependent rainfall observations in addition to the colocated elevation. Last, ordinary kriging yields more accurate predictions than linear regression when the correlation between rainfall and elevation is moderate (less than 0.75 in the case study).

1. Introduction

Measured rainfall data are important to many problems in hydrologic analysis and designs. For example, the ability of obtaining high resolution estimates of spatial variability in rainfall fields becomes important for identification of locally intense storms that could lead to floods and especially flash floods. The accurate estimation of the spatial distribution of rainfall requires a very dense network of instruments, which entails large installation and operational costs. Also, vandalism or the failure of the observer to make the necessary visit to the gage may result in even lower sampling density; thus, it is necessary to estimate point rainfall at unrecorded locations from values at surrounding sites.

A number of methods have been proposed for the interpolation of rainfall data. The simplest approach consists of assigning to the unsampled location the record of the closest gauge (Thiessen, 1911). This method amounts to drawing around each gauge a polygon of influence with the boundaries at a distance halfway between gauge pairs, hence the name Thiessen polygon for the technique. In 1972, the U.S. National Weather Service developed another method whereby the unknown rainfall depth is estimated as a weighted average of surrounding values, the weights being reciprocal to the square distances from the unsampled location (Bedient and Huber, 1992, p.25). Like the Thiessen polygon method, the inverse square distance technique does not allow the hydrologist to consider factors, such as topography, that can affect the catch at a gage. The isohyetal method (McCuen, 1998, p.190) is designed to overcome this deficiency. The idea is to use the location and catch for each gage, as well as knowledge of the factors affecting these catches, to draw lines of equal rainfall depth (isohyets). The amount of rainfall at the unsampled location is then estimated by interpolation within the isohyets. A limitation of the technique is that an extensive gage network is required to draw isohyets accurately.

Geostatistics, which is based on the theory of regionalized variables (Journel and Huijbregts, 1978; Goovaerts, 1997, 1998a), is increasingly preferred because it allows one to capitalize on the spatial correlation between neighboring observations to predict attribute values at unsampled locations. Several authors (Tabios and Salas, 1985; Phillips et al., 1992) have shown that geostatistical prediction techniques (kriging) provide better estimates of rainfall than conventional methods. Recently, Dirks et al. (1998) found that the results depend on the sampling density and that, for high-resolution networks (e.g., 13 rain gauges over a 35 km² area), the kriging method does not show significantly greater predictive skill than simpler techniques, such as the inverse square distance method. Borga and Vizzaccaro (1997) found similar results when they compared kriging and multiquadratic surface fitting for various gauge densities. In fact, besides providing a measure of prediction error (kriging variance), a major advantage of kriging over simpler methods is that sparsely sampled observations of the primary attribute can be complemented by secondary attributes that are more densely sampled. For rainfall, secondary information can take the form of weather-radar observations. A multivariate extension of kriging, known as cokriging, has been used for merging raingage and radar-rainfall data (Creutin et al., 1988; Azimi-Zonooz et al., 1989). Raspa et al. (1997) used another geostatistical technique, kriging with an external drift, to combine both types of information. In this paper, another valuable and cheaper source of secondary information is considered: digital elevation model (DEM). A straightforward approach consists of estimating rainfall at a DEM grid cell through a regression of rainfall versus elevation (Daly et al., 1994). More sophisticated approaches, like cokriging, have been used to incorporate elevation into the mapping of rainfall (Hevesi et al., 1992 a,b).

In this paper, annual and monthly rainfall data from the Algarve region (Portugal) are interpolated using two types of techniques: 1) methods that use only rainfall data recorded at 36 stations (Thiessen polygon, inverse square distance, and ordinary kriging), and 2) algorithms that combine rainfall data with a digital elevation model (linear regression, simple kriging with varying local means, kriging with an external drift, and colocated ordinary cokriging). Prediction performances of the different algorithms are compared using cross validation and are related to the strength of the correlation between rainfall and elevation, and the pattern of spatial dependence of rainfall.

2. Setting the problem

The Algarve is the most southern region of Portugal, with an area of approximately 5,000 km². Figure 1 shows the location of 36 daily-read raingauge stations used in this study. The monthly and annual rainfall depths have been computed for the period of January 1970 to March 1995, and basic sample statistics (mean, standard deviation, minimum, maximum) are given in Table 1.

Another source of information is the elevation map shown in Figure 2. Each grid cell represents 1 km² and its elevation was computed as the average of the elevations at 4 discrete points within the cell. The two main Algarve's mountains dominate the relief: the Monchique (left) and the Caldeirão (right). Table 1 (right column) indicates that the correlation between rainfall and elevation ranges from 0.33 to 0.83, hence it seems worth accounting for this exhaustive secondary information into the mapping of rainfall.

3. Interpolation procedures

This section briefly introduces the different estimators used in the case study. Interested readers should refer to Goovaerts (1997) for a detailed presentation of the different kriging algorithms, and Deutsch and Journel (1998) for their implementation in the public-domain Geostatistical Software Library (Gslib).

3.1 Univariate estimation

Consider first the problem of estimating the rainfall depth z at an unsampled location u using only rainfall data. Let {z(u_i), i=1,...,n} be the set of rainfall data measured at n=36 locations u_i: The most straightforward approach is the Thiessen polygon method whereby the value of the closest observation is simply assigned to u:

Figure 3 (bottom graph) shows the map of annual rainfall interpolated at the nodes of a 1 km x 1 km grid corresponding to the resolution of the elevation model. The map displays the characteristic polygonal zones of influence around the 36 gages.

To avoid unrealistic patchy maps, the rainfall depth z can be estimated as a linear combination of several surrounding observations, with the weights being inversely proportional to the square distance between observations and u:

Figure 4 shows the map of annual rainfall produced using the inverse square distance method and n(u)=16 surrounding observations.

The basic idea behind the weighting scheme (2) is that observations that are close to each other on the ground tend to be more alike than those further apart, hence observations closer to u should receive a larger weight. Instead of the Euclidian distance, geostatistics uses the semivariogram as a measure of dissimilarity between observations. The experimental semivariogram is computed as half the average squared difference between the components of data pairs:

where N(h) is the number of pairs of data locations a vector h apart. The semivariogram is a function of both the distance and direction, and so it can account for direction-dependent variability (anisotropic spatial pattern).

Figure 5 shows the semivariogram of annual rainfall computed from the 36 data of Figure 3. Because of the lack of data only the omnidirectional semivariogram was computed, hence the spatial variability is assumed to be identical in all directions. Semivariogram values increase with the separation distance, reflecting our intuitive feeling that two rainfall data close to each other on the ground are more alike, thus, their squared difference is smaller than those further apart. The semivariogram reaches a maximum at 25 km before dipping and fluctuating around a sill value. The so-called "hole effect" typically reflects pseudo-periodic or cyclic phenomena (Journal and Huijbregts, 1978, p. 403). Here, the hole effect relates to the existence of two mountains 40 km apart (recall Figure 2) which creates two high-valued areas in the rainfall field.

Kriging is a generalized least-square regression technique that allows one to account for the spatial dependence between observations, as revealed by the semivariogram, into spatial prediction. Most of geostatistics is based on the concept of a random function, whereby the set of unknown values is regarded as a set of spatially dependent random variables. Each measurement z(u_i) is thus interpreted as a particular realization of a random variable Z(u_i). Interested readers should refer to textbooks such as Isaaks and Srivastava (1989, pp. 196-236) or Goovaerts (1997, pp. 59-74) for a detailed presentation of the theory of random functions. Like the inverse square distance method, geostatistical interpolation amounts at estimating the unknown rainfall depth z at the unsampled location u as a linear combination of neighboring observations:

The ordinary kriging weights are determined such as to minimize the estimation variance, Var{Z^* _OK(u) - Z(u)}, while ensuring the unbiasedness of the estimator, E{Z^* _OK(u) - Z(u)}=0. These weights are obtained by solving a system of linear equations which is known as ``ordinary kriging system":

The only information required by the kriging system (5) are semivariogram values for different lags, and these are readily derived once a semivariogram model has been fitted to experimental values. Figure 5 shows three different types of permissible models that are combined with a nugget-effect model for the fitting of the experimental semivariogram of annual rainfall:

1. Spherical model with range a:



2. Cubic model with range a:



3. Dampened hole effect model:

The spherical model is the most widely used semivariogram model and is characterized by a linear behaviour at the origin. The cubic model (Galli et al., 1984) allows one to account for the parabolic behaviour displayed by the rainfall semivariogram for the first lags and that is, to a large extent, imparted by elevation that varies smoothly across the study area. This model is preferred to the Gaussian model because it avoids numerical unstability in the kriging system (Wackernagel, 1998, p.120). The last type of function is more complex and used to model hole effect (Deutsch and Journel, 1998, p.26). The three models have been fitted using regression and are such that the weighted sum of squares (WSS) of differences between experimental and model semivariogram values is minimum:

The weighting scheme aims at giving more importance to the first lags and the ones computed from more data pairs. The cubic model has been retained because it yields the smallest WSS value while being more parsimonious (less parameters to estimate) than the dampened hole effect model. Figure 6 shows the rainfall map produced by ordinary kriging using the cubic semivariogram model and the 16 closest observations at each grid node. Like the inverse square distance method, the rainfall map is quite crude, which stresses the importance of accounting for more densely sampled information, such as the digital elevation model of Figure 2.

3.2 Accounting for elevation

Consider now the situation where the rainfall data {z(u_i), i=1,...,n} are supplemented by elevation data available at all estimation grid nodes and denoted y(u).

A straightforward approach consists of predicting the rainfall as a function of the co-located elevation, e.g., using a linear relation such as:

where the two regression coefficients a₀^*, a₁^* are estimated from the set of collocated rainfall and elevation data {(z(u_i), y(u_i)), i=1,...,n}. For example, the relation between annual rainfall and elevation was modeled as z(u) = 324.1 + 0.922 y(u) R² = 0.62, leading to the rainfall map shown in Figure 7. A major shortcoming of this type of regression is that the rainfall at a particular grid node u is derived only from the elevation at u, regardless of the records at the surrounding raingauges u_i. Such an approach amounts at assuming that the residual values r(u_i) = z(u_i) - f[y(u_i)] are spatially uncorrelated. Spatial correlation of the residuals or of the rainfall observations can be taken into account using the three types of geostatistical algorithms described below.

Simple kriging with varying local means (SKlm) amounts at replacing the known stationary mean in the simple kriging estimate by known varying means m_SK^* (u) derived from the secondary information (Goovaerts, 1997, pp. 190-191):

If the local means are derived using a relation of type (10), the SKlm estimate can be rewritten as the sum of the regression estimate and the SK estimate of the residual value at u:

where the weights are obtained by solving the simple kriging system:

where C_R(h) is the covariance function of the residual R(u) = Z(u) - m(u), not that of Z(u) itself. If the residuals are uncorrelated, C_R(h) = 0 for all h, hence all kriging weights in equation (12) are zero and the SKlm estimate is but the value provided by linear regression. Although the kriging system (13) is expressed in terms of covariances, the common practice consists of estimating and modeling the semivariogram, then retrieving the covariance C_R(h) by substracting the semivariogram value from the variance C_R(0). For example, Figure 8 shows the semivariogram of residuals for annual rainfall with the model fitted. This model was used to generate the rainfall map on Figure 9. The impact of elevation on the rainfall map is less pronounced than for the map generated using linear regression that was but a rescaling of the elevation model (Figure 7).

Like the SKlm approach, kriging with an external drift (KED) uses the secondary information to derive the local mean of the primary attribute z, then performs simple kriging on the corresponding residuals:

where m_KED^* (u) = a₀^*(u) + a₁^*(u) y(u). Estimates (11) and (14) differ by the definition of the local mean or trend. The trend coefficients a₀^* and a₁^* are derived once and independently of the kriging system in the SKlm approach, whereas in the KED approach the regression coefficients a₀^*(u) and a₁^*(u) are implicitly estimated through the kriging system within each search neighborhood. In other words, the relation between elevation and rainfall is assessed locally, which allows one to account for changes in correlation across the study area. The usual and equivalent expression for the KED estimate (Wackernagel, 1998, pp. 199-201; Goovaerts, 1997, pp. 194-198) is:

The kriging weights are the solution of the following system of (n(u) + 2) linear equations:

with 2 Lagrange parameters accounting for the constraints on the weights. Figure 10 shows the map of KED estimates of annual rainfall, which looks very similar to the SKlm map.

Another approach for incorporating secondary information is cokriging, a multivariate extension of kriging (Goovaerts, 1997, pp. 203-248). When the secondary variable is known everywhere and varies smoothly across the study area (e.g., elevation) there is little loss in retaining in the cokriging system, only the secondary datum co-located with the location u being estimated (Xu et al., 1992; Goovaerts, 1998b). Indeed the co-located elevation y(u) tends to screen the influence of further away elevation data. The cokriging estimate is then:

where m_Z and m_Y are the global means of the rainfall and elevation data, see Table 1. The main difference between cokriging and the two previous geostatistical algorithms lies in how the elevation value is handled. Whereas that datum directly influences the cokriging estimate, in the SKlm and KED approaches, elevation provides information only about the primary trend at location u. The cokriging weights are solutions of the following system of (n(u) + 2) linear equations:

where C_ZY(u_i - u) is the cross covariance between primary and secondary variables at locations u_i and u, respectively. Again, the common practice consists of estimating and modeling the (cross) semivariogram, then retrieving the (cross) covariance. Figure 11 shows the experimental semivariograms of elevation and annual rainfall, and their cross semivariogram computed as:

A linear model of coregionalization was fitted using an iterative procedure developed by Goulard (1989). The three semivariograms are modeled as a linear combination of the same set of basic models (here a nugget effect and a cubic model with range 29.18 km) so that the WSS criterion is minimum under the constraints of positive semi-definiteness of the matrix of sills, see Goovaerts (1998a) for more details. Figure 12 shows the map of cokriging estimates of annual rainfall. Unlike the three previous techniques, the details of the elevation map do not appear in the rainfall map that actually shows more similarities with the ordinary kriging map of Figure 6.

4. Evaluation of the different interpolators

The performances of the seven interpolators were assessed and compared using cross validation (Isaaks and Srivastava, 1989, pp. 351-368). The idea consists of removing temporarily one rainfall observation at a time from the data set and ``re-estimate" this value from remaining data using the alternative algorithms. The comparison criterion is the mean square error (MSE) of prediction, which measures the average square difference between the true rainfall z(u_i) and its estimate z^*(u_i):

where n=36 for the Algarve data set. The value of this criterion should be close to zero if the algorithm is accurate. For linear regression, the MSE was simply computed as the average square residual value for the linear model fitted using all 36 observations, which means that the prediction error would tend to be underestimated for this method.

Figure 13 shows the mean square errors of prediction produced by each of the seven interpolation algorithms for the monthly (1 to 12) and annual (13) rainfall. Results are expressed as proportions of the prediction error of the linear regression approach. The conclusions are as follows:

1. Larger predictions errors are obtained for the three algorithms that ignore elevation, with the worst results produced by Thiessen polygon. It is noteworthy that for several months, and on average over the year, ordinary kriging yields smaller prediction errors than linear regression of rainfall against elevation.
2. Except for the period from June to September, which is characterized by low rainfall amounts (Table 1), multivariate geostatistical algorithms perform better than linear regression, which disregards the information provided by surrounding climatic stations.
3. Among geostatistical algorithms that account for elevation, simple kriging with varying local means and kriging with an external drift yield slightly better results than the more complex ordinary cokriging.

To identify the factors that might be responsible for these relative prediction performances, ratios of MSE values have been plotted against parameters, such as the relative nugget effect of the rainfall semivariogram, the relative nugget effect of the residual semivariogram, and the correlation coefficient between rainfall and elevation. The first scattergram (Figure 14, left top graph) shows that the benefit of using ordinary kriging instead of the inverse square distance method (i.e., a larger MSE ratio Inv/OK) decreases as the spatial dependence between observations weakens, which is indicated by a larger relative nugget effect for the rainfall semivariogram. Similarly, the benefit of using a geostatistical approach (SKlm) instead of linear regression to account for elevation decreases as the spatial dependence between residual data weakens (larger nugget effect for the residual semivariogram), see Figure 14 (right top graph).

The two bottom graphs of Figure 14 show the impact of the strength of the correlation between elevation and rainfall on the relative performances of ordinary kriging versus cokriging (left graph) and ordinary kriging versus linear regression (right graph). Note that the smallest correlation coefficients observed during the summer have not been taken into account because they also correspond to negligible amounts of rainfall, see Table 1. The gain of ordinary cokriging versus kriging increases as elevation brings more information on rainfall (higher correlation coefficient). The other graph indicates that ordinary kriging performs better than linear regression (ratio smaller than 1) when the correlation between elevation and rainfall is moderate, i.e., less than 0.75 and the information from surrounding stations is worth integration. For larger correlation coefficients, the elevation at the estimation grid node screens the influence of rainfall data at surrounding sites, so spatial information is less valuable.

5. Conclusions

Our results confirm previous findings that for low-resolution networks, geostatistical interpolation outperforms techniques, such as the inverse square distance or Thiessen polygon, and ignore the pattern of spatial dependence that is usually observed for rainfall data. Prediction can be further improved if correlated secondary information, such as a digital elevation model, is taken into account. This paper has reviewed different ways to incorporate such exhaustive secondary information, and cross validation has shown that prediction performances can vary greatly among algorithms.

The most straightforward approach consists of deriving the rainfall value directly from the colocated elevation through a (non)linear regression. By so doing, one ignores the information provided by surrounding climatic stations that are critical when the correlation between the two variables is not too strong and when the residuals are spatially correlated. An easy way to account for both elevation and spatial correlation is to interpolate the regression residuals using geostatistics, that is, to use simple kriging with varying local means. A more sophisticated approach consists of assessing the correlation between elevation and rainfall within each search neighborhood (kriging with an external drift). The last technique is cokriging, which interpolates the rainfall as a linear combination of surrounding rainfall observations and the colocated elevation. This approach is the most demanding because three semivariograms must be inferred and jointly modeled, a task that is alleviated by the recent development of automatic fitting procedures. Cokriged maps show less details than the SKlm and KED maps that are greatly influenced by the pattern of the DEM. In this case study, the additional complexity of cokriging does not pay off.

Finally, this paper has introduced the little known cubic semivariogram model that allows the modeling of parabolic behaviour typically observed at the origin of the semivariograms of elevation and correlated rainfall, while avoiding unstability problems frequently encountered when using the Gaussian semivariogram model in kriging.

Acknowledgments

The author thanks Mr. Nuno de Santos Loureiro for the Algarve data set.

References

Azimi-Zonooz, A., W.F. Krajewski, D.S. Bowles and D.J. Seo, 1989. Spatial rainfall estimation by linear and non-linear cokriging of radar-rainfall and raingage data, Stochastic Hydrol. Hydraul., 3, pp. 51-67.

Bedient, P.B. and W.C. Huber, 1992. Hydrology and Floodplain Analysis: Second Edition, Addison-Wesley, New York, NY.

Borga, M. and A. Vizzaccaro, 1997. On the interpolation of hydrologic variables: formal equivalence of multiquadratic surface fitting and kriging, J. Hydrol., 195(1-4), pp. 160-171.

Creutin, J.D., G. Delrieu and T. Lebel, 1988. Rain measurement by raingage-radar combination: a geostatistical approach, J. Atm. and Oceanic Tech., 5(1), pp. 102-115.

Daly, C., R.P. Neilson and D.L. Phillips, 1994. A statistical topographic model for mapping climatological precipitation over mountainous terrain, J. Appl. Meteor., 33(2), pp. 140-158.

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

Dirks, K.N., J.E. Hay, C.D. Stow and D. Harris, 1998. High-resolution studies of rainfall on Norfolk Island Part II: interpolation of rainfall data, J. Hydrol., 208(3-4), pp. 187-193.

Galli, A., F. Gerdil-Neuillet and C. Dadou, 1984. Factorial kriging analysis: a substitute to spectral analysis of magnetic data. In: G. Verly, M. David, A. G. Journel, A. Maréchal (Editors). Geostatistics for Natural Resources Characterization, Reidel, Dordrecht, pp. 543-557.

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

Goovaerts, P., 1998a. Geostatistical tools for characterizing the spatial variability of microbiological and physico-chemical soil properties, Biol. Fertil. Soils, 27(4), pp. 315-334.

Goovaerts, P., 1998b. Ordinary cokriging revisited, Math. Geol., 30(1), pp. 21-42.

Goulard, M., 1989. Inference in a coregionalization model. In: M. Armstrong, Ed. Geostatistics, Kluwer, Dordrecht, pp. 397-408.

Hevesi, J.A., J.D. Istok and A.L. Flint, 1992a. Precipitation estimation in mountainous terrain using multivariate geostatistics. Part I: Structural analysis, J. Appl. Meteor., 31, pp. 661-676.

Hevesi, J.A., A.L. Flint and J.D. Istok, 1992b. Precipitation estimation in mountainous terrain using multivariate geostatistics. Part II: Isohyetal maps, J. Appl. Meteor., 31, pp. 677-688.

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

Journel, A.G. and C.J. Huijbregts, 1978. Mining Geostatistics, Academic Press, New York, NY.

McCuen, R.H., 1998. Hydrologic Analysis and Design, 2nd edition, Prentice Hall, NJ.

Phillips, D.L., J. Dolph and D. Marks, 1992. A comparison of geostatistical procedures for spatial analysis of precipitations in mountainous terrain, Agric. and Forest Meteor., 58, pp. 119-141.

Raspa, G., M. Tucci and R. Bruno, 1997. Reconstruction of rainfall fields by combining ground raingauges data with radar maps using external drift method. In: E.Y. Baafi, N.A. Schofield (Eds.). Geostatistics Wollongong '96, Kluwer Academic Publishers, Dordrecht, pp. 941-950.

Tabios, G.Q. and J.D. Salas, 1985. A comparative analysis of techniques for spatial interpolation of precipitation, Water Resources Bulletin, 21(3), pp. 365-380.

Thiessen, A.H., 1911. Precipitation averages for large areas, Monthly Weather Review, 39(7), pp. 1,082-1,084.

Wackernagel, H., 1998. Multivariate Geostatistics: Second, completely revised Edition, Springer-Verlag, Berlin, Germany.

Xu, W., T.T. Tran, R.M. Srivastava and A.G. Journel, 1992. Integrating seismic data in reservoir modeling: the collocated cokriging alternative, Society of Petroleum Engineers, paper no. 24,742.