Download Kriging Example: Spatial Interpolation and Prediction and more Summaries Statistics in PDF only on Docsity!
Kriging Example
STAT 197 - Intro to Spatial Statistics
Kriging (Matheron 1963) is a spatial interpolation approach for making predictions at unsampled places using geostatistical data. This method was developed in the field of mining geology and is named after South African mining engineer Danie G. Krige. Suppose that we have observed data Z ( s 1 ) , ..., Z ( sn ), and wish to predict the value of Z at an arbitrary location s 0 ∈ D. The Ordinary Kriging estimator of Z ( s 0 ) is defined as the linear unbiased estimator
Z^ ˆ( s 0 ) =
∑^ n
i =
λiZ ( si )
that minimizes the mean squared prediction error defined as
E [( Z ˆ( s 0 ) − Z ( s 0 ))^2 ].
The Kriging weights are based on the sampled data’s predicted spatial structure. Weights are calculated by first fitting a variogram model to the observed data, which allows us to understand how the correla- tion between observation values varies with distance between places. Once the Kriging weights have been determined, they are applied to the known data values at the sampled locations to compute the predicted results at unsampled places. The Kriging weights describe the spatial correlation of the data, taking into consideration the geographical proximity and similarity of data points. Thus, observed locations that are associated and close to the predicted locations receive more weight than those that are uncorrelated and farther apart. Weights also consider the spatial arrangement of all observations, so clusters of observations in oversampled areas are weighted less heavily because they carry less information than individual sites.
Kriging Predictions of Zinc Concentrations
Data
The meuse data from the sp package contains zinc and other soil-heavy metal concentrations collected at locations in a flood plain of the river Meuse near Stein, The Netherlands. The meuse.grid contains prediction grid locations for the meuse dataset. We convert the meuse and meuse.grid data frames to sf objects using the st_as_sf() function and setting the coordinate reference system to EPSG 28992.
#install.packages('gstat')
library (sp) library (gstat) library (sf) library (ggplot2) library (viridis)
data (meuse) data (meuse.grid)
meuse <- st_as_sf (meuse, coords = c ("x", "y"), crs = 28992)
meuse.grid <- st_as_sf (meuse.grid, coords = c ("x", "y"), crs = 28992)
ggplot () + geom_sf (data = meuse, aes (color = meuse $ zinc)) + scale_color_viridis (name = "Zinc Content") + theme_bw ()
50.955°N
50.960°N
50.965°N
50.970°N
50.975°N
50.980°N
50.985°N
50.990°N
5.73°E 5.74°E 5.75°E 5.76°E
500
1000
1500
Zinc Content
Variogram Cloud
The variogram cloud shows half of all possible squared differences of data observation pairs against their separation distance h:
( Z ( s ) − Z ( s + h ))^2
The variogram() function of gstat can be used to calculate the variogram cloud. In variogram(), we set argument object to z ~ 1 if we wish to obtain the variogram for data z, or to a formula of a linear model with covariates if we wish the variogram for the residuals. This plot can be inspected to assess whether sample pairs closer to each other are more similar than pairs further apart.
distance
semivariance
Fitted Variogram
The sample variogram obtained is a function of distance h estimated at discrete lags (e.g., h (1) , h (2) , ..., h ( k )). Then, we can fit a variogram model to these estimated values,
{( h ( j ) , 2ˆ γ ( h ( j ))) : j = 1 , 2 , ..., k }.
The vgm() function of gstat generates a variogram model. This function has arguments model with the model type (e.g., Sph for spherical and Exp for exponential), nugget, psill (partial sill, which is the sill minus the nugget) and range. A list of available models can be printed with vgm() and visualized with show.vgms().
vgm ()
short long
1 Nug Nug (nugget)
2 Exp Exp (exponential)
3 Sph Sph (spherical)
4 Gau Gau (gaussian)
5 Exc Exclass (Exponential class/stable)
6 Mat Mat (Matern)
7 Ste Mat (Matern, M. Stein’s parameterization)
8 Cir Cir (circular)
9 Lin Lin (linear)
10 Bes Bes (bessel)
11 Pen Pen (pentaspherical)
12 Per Per (periodic)
13 Wav Wav (wave)
14 Hol Hol (hole)
15 Log Log (logarithmic)
16 Pow Pow (power)
17 Spl Spl (spline)
18 Leg Leg (Legendre)
19 Err Err (Measurement error)
20 Int Int (Intercept)
show.vgms (par.strip.text = list (cex = 0.75))
distance
semivariance
vgm(1,"Nug",0)
vgm(1,"Exp",1) vgm(1,"Sph",1)
vgm(1,"Gau",1) vgm(1,"Exc",1)
vgm(1,"Mat",1)
vgm(1,"Ste",1) vgm(1,"Cir",1) vgm(1,"Lin",0) vgm(1,"Bes",1) vgm(1,"Pen",1)
vgm(1,"Per",1)
vgm(1,"Wav",1) vgm(1,"Hol",1)
vgm(1,"Log",1) vgm(1,"Pow",1)
vgm(1,"Spl",1)
The fit.variogram() function fits a variogram model to a sample variogram. By inspecting the sample variogram of the zinc data, we may think that a good model for the variogram would be spherical with initial values for the partial sill equal to 0.5, range equal to 900, and nugget equal to 0.1. Plot below shows the sample variogram v calculated with variogram(), and the variogram generated with vgm() using a spherical model, and with partial sill, range, and nugget equal to our initial guess values. This plot allows us to assess our initial model.
vinitial <- vgm (psill = 0.5, model = "Sph", range = 900, nugget = 0.1) plot (v, vinitial, cutoff = 1000, cex = 1.5)
distance
semivariance
Kriging
Kriging uses the fitted variogram values fv to obtain predictions at unsampled locations. Here, we use the gstat() function of gstat to compute the Kriging predictions.
k <- gstat (formula = log (zinc) ~ 1, data = meuse, model = fv)
Then, we obtain predictions at the unsampled locations given by meuse.grid with the predict.gstat() function.
kpred <- predict (k, meuse.grid)
[using ordinary kriging]
The returned object of predict() contains the predictions (var1.pred) and their variance (var1.var) which allow us to quantify the uncertainty of these predictions. Below shows maps with the predictions and variance values. Note that the Kriging predictions are calculated using weights that depend on the variogram estimates, and they are therefore sensitive to misspecification of the variogram model.
ggplot () + geom_sf (data = kpred, aes (color = var1.pred)) + geom_sf (data = meuse) + scale_color_viridis (name = "log(zinc)") + theme_bw ()
50.955°N
50.960°N
50.965°N
50.970°N
50.975°N
50.980°N
50.985°N
50.990°N
5.72°E 5.73°E 5.74°E 5.75°E 5.76°E
log(zinc)
ggplot () + geom_sf (data = kpred, aes (color = var1.var)) + geom_sf (data = meuse) + scale_color_viridis (name = "variance") + theme_bw ()
