Geostatistics — Variograms and Kriging (R)¶
gstat on the Meuse data¶
The R counterpart to kriging_python.ipynb. Where the Python notebook builds the variogram and kriging equations from scratch, this uses gstat (Pebesma) — the standard geostatistics package: variogram and fit.variogram for the spatial-dependence model, krige for ordinary and universal kriging, and krige.cv for leave-one-out validation. Same Meuse soil-zinc data, with sp/spplot maps.
options(repr.plot.width=6, repr.plot.height=4.5)
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths()))
suppressMessages({library(sp); library(gstat)})
data(meuse); coordinates(meuse) <- ~x+y
data(meuse.grid); coordinates(meuse.grid) <- ~x+y; gridded(meuse.grid) <- TRUE
cat(nrow(meuse), "soil samples;", length(meuse.grid), "grid cells; zinc", min(meuse$zinc), "-", max(meuse$zinc), "ppm\n")
print(spplot(meuse, "zinc", main="Meuse zinc samples (ppm)", col.regions=heat.colors(20, rev=TRUE), cex=0.8))
Warning message: "package 'gstat' was built under R version 4.6.1"
155 soil samples; 3103 grid cells; zinc 113 - 1839 ppm
1. The variogram¶
variogram(log(zinc) ~ 1, ...) computes the empirical semivariogram; fit.variogram fits an exponential model (nugget, partial sill, range).
v <- variogram(log(zinc) ~ 1, meuse)
vfit <- fit.variogram(v, vgm(0.6, "Exp", 300, 0.05))
print(vfit)
print(plot(v, vfit, main="Empirical and fitted variogram of log-zinc"))
model psill range 1 Nug 0.0000000 0.000 2 Exp 0.7186526 449.758
2. Ordinary kriging¶
krige(log(zinc) ~ 1, ...) predicts the surface and returns the kriging variance. We map both, and validate with leave-one-out cross-validation (krige.cv).
ok <- krige(log(zinc) ~ 1, meuse, meuse.grid, model=vfit, debug.level=0)
cv <- krige.cv(log(zinc) ~ 1, meuse, model=vfit, verbose=FALSE)
cat(sprintf("ordinary kriging LOO RMSE (log-zinc): %.3f\n", sqrt(mean(cv$residual^2))))
print(spplot(ok, "var1.pred", main="Ordinary kriging: predicted log-zinc", col.regions=heat.colors(40, rev=TRUE)))
print(spplot(ok, "var1.var", main="Kriging variance (uncertainty)", col.regions=cm.colors(40)))
ordinary kriging LOO RMSE (log-zinc): 0.393
3. Universal kriging with a river-distance trend¶
Adding sqrt(dist) (distance to the river) as a trend and kriging the residuals — the universal-kriging model — improves prediction, as leave-one-out cross-validation confirms.
vu <- variogram(log(zinc) ~ sqrt(dist), meuse)
vufit <- fit.variogram(vu, vgm(0.4, "Exp", 300, 0.05))
uk <- krige(log(zinc) ~ sqrt(dist), meuse, meuse.grid, model=vufit, debug.level=0)
cvu <- krige.cv(log(zinc) ~ sqrt(dist), meuse, model=vufit, verbose=FALSE)
cat(sprintf("universal kriging LOO RMSE: %.3f vs ordinary %.3f\n", sqrt(mean(cvu$residual^2)), sqrt(mean(cv$residual^2))))
print(spplot(uk, "var1.pred", main="Universal kriging (river-distance trend): predicted log-zinc", col.regions=heat.colors(40, rev=TRUE)))
universal kriging LOO RMSE: 0.378 vs ordinary 0.393
4. Summary¶
gstat reproduces the from-scratch geostatistics. variogram/fit.variogram estimate the spatial-dependence model, krige produces the ordinary- and universal-kriging surfaces with their kriging variances, and krige.cv confirms that adding the river-distance trend (universal kriging) lowers the leave-one-out error — matching kriging_python.ipynb. The spplot maps show the same high-zinc river corridor and the uncertainty growing away from the samples.
gstat (with sp/sf) is the field-standard geostatistics toolkit and the package mirror of the from-scratch kriging equations, which are themselves Gaussian-process regression (the PyMC GP in the Python notebook). This point-referenced subsection complements the areal (CAR/SAR) and spatial-econometric models on regions. Next in the arc: spatial point processes, then the spatiotemporal capstone.