0.1 libraries

library(sp)
library(gstat)
library(sf)
library(mapview)
library("stars")  
library(ggplot2)
library(viridis)

meuse = point data (cadmium, zinc, copper, etc.)

meuse.grid = grid for interpolation

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)
mapview(meuse, zcol = "zinc",  map.types = "CartoDB.Voyager")
mapview(meuse.grid,  map.types = "CartoDB.Voyager")
##this command should be memorised

Cloud variogram(many dots): Shows raw variability

vc <- variogram(log(zinc) ~ 1, meuse, cloud = TRUE)
plot(vc)

Sample variogram (bin-averaged): Used to fit a variogram model

v <- variogram(log(zinc) ~ 1, data = meuse)
plot(v)

** Variogram models**

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)
vinitial <- vgm(psill = 0.5, model = "Sph",
                range = 900, nugget = 0.1)
vinitial
##   model psill range
## 1   Nug   0.1     0
## 2   Sph   0.5   900
plot(v, vinitial, cutoff = 1000, cex = 1.5)

fv <- fit.variogram(object = v,
                    model = vgm(psill = 0.5, model = "Sph",
                                range = 900, nugget = 0.1))
fv
##   model      psill    range
## 1   Nug 0.05066017   0.0000
## 2   Sph 0.59060556 897.0044
plot(v, fv, cex = 1.5)

1 Kriging

Create a gstat object for kriging and Predict zinc values on a grid ~1 means ordinary kriging (no trend)

k <- gstat(formula = log(zinc) ~ 1, data = meuse, model = fv)

kpred <- predict(k, meuse.grid)
## [using ordinary kriging]

**Plot predicted zinc values*

ggplot() + geom_sf(data = kpred, aes(color = var1.pred)) +
  geom_sf(data = meuse) +
  scale_color_viridis(name = "log(zinc)") + theme_bw()

Plot kriging variance (uncertainty)

ggplot() + geom_sf(data = kpred, aes(color = var1.var)) +
  geom_sf(data = meuse) +
  scale_color_viridis(name = "variance") + theme_bw()

2 Grid

st_bbox() gives the minimum and maximum x–y coordinates of the meuse dataset.

bbox <- st_bbox(meuse)
bbox
##   xmin   ymin   xmax   ymax 
## 178605 329714 181390 333611

Set grid cell size

cell_size <- 40

** Create sequences of x and y coordinates**

x <- seq(bbox$xmin, bbox$xmax, by=cell_size)
y <- seq(bbox$ymin, bbox$ymax, by=cell_size)

Create a full grid with every possible combinations of x–y

meuse_grid <- expand.grid(x=x, y=y)

Plot the grid points

plot(meuse_grid$x, meuse_grid$y, pch=19, cex=0.1)

Add a dummy variable: st_as_stars() requires at least one data column to convert to a stars raster.

meuse_grid$tmp <- 1
meuse_grid <- st_as_stars(meuse_grid, crs=st_crs(meuse))
st_crs(meuse_grid) <- st_crs(meuse)