1. Introduction

This R Markdown notebook illustrates kriging interpolation.

The code was written by Barry Boessenkool. Thank you Barry!

2. Toy interpolation

A single script to understand interpolation:

if(!requireNamespace("pacman", quietly=TRUE)) install.packages("pacman")
pacman::p_load(geoR, berryFunctions) # installs packages if needed
# berryFunctions: seqPal, distance, seqR, colPoints
# geoR: as.geodata, variofit, variog, krige.conv, krige.control
K_X <- c(4,4,5,5,6,6,6,7,7,8,9,9,9)
K_Y <- c(6,9,5,8,4,6,8,4,8,7,3,7,9)
K_Z <- c(5,9,2,6,3,5,9,4,8,8,3,6,7)
K_COLS <- seqPal()
K_MDIST <- max(distance(K_X,K_Y))
K_RESOLUTION <- median(distance(K_X,K_Y))/10
K_GEODATA <- as.geodata(cbind(K_X,K_Y,K_Z))
K_VARIOGRAM <- variofit(variog(K_GEODATA, max.dist=K_MDIST))
## variog: computing omnidirectional variogram
## variofit: covariance model used is matern 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim
## Warning in variofit(variog(K_GEODATA, max.dist = K_MDIST)): initial values not
## provided - running the default search
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq kappa
## initial.value "6.8"   "2.83" "0"   "0.5"
## status        "est"   "est"  "est" "fix"
## loss value: 220.5968039392
K_GRID <- expand.grid(seqR(K_X, by=K_RESOLUTION, extend=0.1),
seqR(K_Y, by=K_RESOLUTION, extend=0.1))
K_OK <- krige.conv(K_GEODATA, locations=K_GRID,
krige=krige.control(type.krige="OK", obj.model=K_VARIOGRAM))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
geoR:::image.kriging(K_OK, col=K_COLS)
colPoints(K_X, K_Y, K_Z, col=K_COLS, zlab="Kriging", cex=1)
points(K_X, K_Y)

# optional: interactive map:
pacman::p_load(raster, mapview) # installs packages if needed
K_RASTER <- raster::rasterFromXYZ(cbind(K_GRID, temp=K_OK$predict),
crs="+proj=longlat +datum=WGS84 +no_defs")
mapview::mapview(K_RASTER, col.regions=seqPal(256))

3. Real-world interpolation

The data used is a shapefile with sampling precipitation data. It can be downloaded from Barry repository at Github.

Use the following link: https://github.com/brry/course/blob/master/data/PrecBrandenburg.zip

Reading precipitation sampled data

# Shapefile:
p <- sf::st_read("datos/PrecBrandenburg/niederschlag.shp", quiet=TRUE)
p
 # Plot prep
pcol <- colorRampPalette(c("red","yellow","blue"))(50)
clss <- berryFunctions::classify(p$P1, breaks=50)$index
# Plot
par(mar = c(1, 1, 1.2, 1)) # for sf plots / colPointsLegend sf:::plot.sf(p, col=pcol[clss], max.plot=1) # P1: Precipitation # kriging coordinates
(cent <- sf::st_centroid(p))
## Warning in st_centroid.sf(p): st_centroid assumes attributes are constant over
## geometries of x
#berryFunctions::colPoints(cent$x, cent$y, p$P1, add=T, cex=0.7, legargs=list(y1=0.8,y2=1), col=pcol)
#points(cent$x, cent$y, cex=0.7)

Variogram

 library(geoR)
# Semivariance:
geoprec <- as.geodata(cbind(cent$x,cent$y,p$P1)) 
vario <- variog(geoprec, max.dist=130000)
## variog: computing omnidirectional variogram
## variog: computing omnidirectional variogram
fit <- variofit(vario)
## variofit: covariance model used is matern 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim
## Warning in variofit(vario): initial values not provided - running the default
## search
## variofit: searching for best initial value ... selected values:
##               sigmasq   phi        tausq kappa
## initial.value "1326.72" "19999.93" "0"   "0.5"
## status        "est"     "est"      "est" "fix"
## loss value: 107266266.76371

Variogram visualization

plot(vario) ; lines(fit)

# distance to closest other point:
d <- sapply(1:nrow(cent), function(i) min(berryFunctions::distance(
                          cent$x[i], cent$y[i], cent$x[-i], cent$y[-i])))
hist(d/1000, breaks=20, main="distance to closest gauge [km]")

mean(d) # 8 km - grid resolution in next slide
## [1] 8165.633

Kriging

# Kriging:
res <- 1000 # 1 km, since stations are 8 km apart on average 
grid <- expand.grid(seq(min(cent$x),max(cent$x),res),
                    seq(min(cent$y),max(cent$y),res))
krico <- krige.control(type.krige="OK", obj.model=fit)
krobj <- krige.conv(geoprec, locations=grid, krige=krico)
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
# Set values outside of Brandenburg to NA:
grid_sf <- sf::st_as_sf(grid, coords=1:2, crs=sf::st_crs(p))
isinp <- sapply(sf::st_within(grid_sf, p), length) > 0
krobj2 <- krobj
krobj2$predict[!isinp] <- NA

Kriging Visualization

geoR:::image.kriging(krobj2, col=pcol)

berryFunctions::colPoints(cent$x, cent$y, p$P1, col=pcol, zlab="Prec", cex=0.7,
     legargs=list(y1=0.1,y2=0.8, x1=0.78, x2=0.87, horizontal=F))

plot(p, add=T, col=NA, border=8)#; points(cent$x,cent$y, cex=0.8)
## Warning in plot.sf(p, add = T, col = NA, border = 8): ignoring all but the first
## attribute