Mapping Temperature with Kriging

Study Site

I have a data set with 110 weather stations across the central Cascade Mountains of Oregon. Each station records min/max temperature and precipitation on a daily time-step. For this example I will consider the annual means of the minimum temperature (T_min) records. Fig 1 shows the positions of the data, with contours indicating the hygric continentality (a measurement of the climatic regime) of the site. The Eastern side of the area has a dry, thermally variable, continental climate, whilst the Western side has a wetter, more thermally stable maritime climate. This climatic gradient, along with the considerable altitudinal variation present (Fig 2), makes the site a challenging and interesting example for meteorological interpolation.

Fig 1: Weather station locations; black outline is the region of interest, contours represent hygric continentality

We also have a DEM (digital elevation model) of the region, derived from radar interferometry (Fig 2). The elevation is recorded in metres. The projection is UTM zone 10, WGS84 datum with units of km. The highest peak is Mount Jefferson, at approximately 590,4950. To the south we have Three Fingered Jack, at approximately 590,4925. South of three fingered Jack is Mount Washington (590,4915). The small round peak to the East of these at about 610,4919 is Black Bute. The whole region is covered in coniferous forest, and is very beautiful.

Fig 2: DEM of the region of interest outlined above. Axes in Km, Elevation in metres.

Interpolation

Ordinary Kriging

To fill the gaps between the stations, we require an interpolation. Kriging allows us to make optimal (in the least squares sense) prediction of the temperature at unsampled points. we proceed by selecting a point in space for which no observation exists. We then predict the datum value from a weighted linear combination of the observed data set. How do we decide how much weight to put on each sample location? Logically, we might think it is a good idea to weight data points closest to the prediction location more heavily than those further away, having correctly observed that things close together tend to be similar.

The problem then is, what kind of law do we follow to decrease the weighting with distance? Simple methods define an arbitrary law, such as an inverse square relationship, or a linear falloff. Whilst this may be sufficient for certain known physical relationships, it does not adequately reflect the complexity of many natural phenomena.

Kriging solves this problem by using a mapping of the statistical difference between data points onto their physical separation distances to generate a set of weights. Such a mapping is known as the semi-variogram (fig 3). Distances separating each of the observations are calculated, and a measure of the statistical difference is calculated for each pair. In order for this relationship to make sense, we require that the mean and variance are invariant across the region of interest, as the presence of large scale trends in the data severely compromises the ability to estimate a stable spatial variance function.

Fig 3: Variogram of T_min data; similarity decreases with separation distance. Points indicate the lag averages (empirical semivariogram), whilst the red line indicates the variogram model. Distances are in km.

By measuring the distance between the point to be estimated, and all the observation locations, we can then build a set of weights from our variogram model which results in a statistically optimal estimation of the value of the parameter of interest at that location. Such an estimate is known as the best linear unbiased estimate (BLUE). By repeating the procedure at every prediction location, we can fill our area of interest with predictions to produce a map (fig 4).

Fig 4: Map of estimated mean annual T_min (°C) via ordinary Kriging.

External Drift Kriging

So far so good then; ordinary kriging allowed us to make planar enforced (that is, spatially exhaustive) estimates of mean annual T_min for our region of interest. A quick comparison of the original data plot (Fig 1) with our estimate seems reasonable: The north-east of the region is warmer than the south-west. We can also rationalise this trend with respect to the topography; in fig 2 we see the mountain range runs through the west of the region from north to south. The highest elevations are in the south-west of the region around Mount Jefferson and Three Fingered Jack. As we expect, the mountain tops are colder than the low lying regions.

Due to the mountainous characteristics of the region, however, we are unlikely to make good predictions of temperature without informing the interpolation of the effect of altitude on temperature. It so happens, in this region, we lose around 3.5 °C per km rise in elevation (fig 5). Ordinary kriging assumes that the global mean (and variance) of the region is constant; a property known as first order stationarity. As we see from fig 5, our region has a variable local mean due to the presence of a trend or "drift". Our data are therefore non-stationary, which complicates matters a little.

Fig5: The relationship between elevation and T_min.

We can allow for this in our predictions using a more exotic interpolation technique called "External Drift Kriging". Here the local mean temperature is inferred from the elevation trend in the data, and small scale variation is incorporated by interpolating the residual variation (fig 6). The results here appear much more reasonable in terms of the known physical relationships, which are imposed on the local mean by the elevation covariate. This partitioning of the data into a deterministic trend component, and a residual 'noise' component is the rationale of non-stationary geostatistics.

Fig 6: Map of estimated mean annual T_min (°C) via external drift Kriging.

References

Excellent introductory material can be found in Isaaks and Srivastava (1989), and Webster and Oliver (2001). Goovaerts (1997) presents a thorough and complete account. Deutsch and Journel (1998) provide software and users information for all algorithms, including advanced non-stationary methods and 3D kriging; technical issues are also well covered, with example data and exercises. Noel Cressie (1993) presents an exhaustive and very mathematical text on all aspects of spatial statistics. Finally, mulivariate implementations and non-stationary methods are introduced and developed in Wackernagel (2003).

• Cressie, N.A.C. (1993). Statistics for Spatial Data (Revised Edition). John Wiley and Sons LTD.
• Deutch, C.V. and Journel, A.G. (1998). GSLIB: Geostatistical Software Library and User's Guide (2nd Ed). Oxford University Press.
• Goovaerts, P. (1997). Geostatistics for Natural Resources Evaluation. Oxford University Press.
• Goovaerts, P., 2000. Geostatistical approaches for incorporating elevation into the spatial interpolation of rainfall. Journal of Hydrology, 228(1-2): 113-129.
• Hudson, G., and H. Wackernagel. 1994. Mapping temperature using kriging with external drift - theory and an example from Scotland. International Journal of Climatology 14:77-91.
• Isaaks, E.H. and Srivastava, E.H. (1989). An Introduction to Applied Geostatistics. Oxford University Press.
• Wackernagel, H. (2003). Multivariate Geostatistics; An Inroduction with Applications (3rd Ed). Springer.
• Webster, R. and Oliver, M.A. (2001). Geostatistics for Environmental Scientists. John Wiley and Sons LTD.