Chapter 16 Spatial Interpolation and geostatistics

16.1 Load meuse data

# this reloads meuse as data.frame, so
library(sp)
demo(meuse, ask = FALSE, echo = FALSE)
library(sf)
meuse_sf = st_as_sf(meuse)

16.2 fit variogram

library(gstat)
v = variogram(log(zinc)~1, meuse_sf)
(v.fit = fit.variogram(v, vgm(1, "Sph", 900, 1)))
#>   model  psill range
#> 1   Nug 0.0507     0
#> 2   Sph 0.5906   897

16.3 kriging to point locations:

(k_sf = krige(log(zinc)~1, meuse_sf, meuse_sf, v.fit))
#> [using ordinary kriging]
#> Simple feature collection with 155 features and 2 fields
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 179000 ymin: 330000 xmax: 181000 ymax: 334000
#> epsg (SRID):    NA
#> proj4string:    +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.2369,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812 +units=m +no_defs
#> First 10 features:
#>    var1.pred  var1.var              geometry
#> 1       6.93  1.96e-33 POINT (181072 333611)
#> 2       7.04  1.96e-33 POINT (181025 333558)
#> 3       6.46  0.00e+00 POINT (181165 333537)
#> 4       5.55  1.96e-33 POINT (181298 333484)
#> 5       5.59 -1.11e-16 POINT (181307 333330)
#> 6       5.64  1.96e-33 POINT (181390 333260)
#> 7       5.85  1.11e-16 POINT (181165 333370)
#> 8       6.01  1.11e-16 POINT (181027 333363)
#> 9       5.85 -2.22e-16 POINT (181060 333231)
#> 10      5.21 -1.11e-16 POINT (181232 333168)
plot(k_sf[1], pch = 16)

16.4 kriging to stars grid:

library(stars)
meuse_stars = st_as_stars(meuse.grid)
k_st = krige(log(zinc)~1, meuse_sf, meuse_stars, v.fit)
#> [using ordinary kriging]
plot(k_st, breaks = "equal", col = sf.colors())

16.5 Spatio-temporal

see gstat/tests/windst.R see gstat/tests/stars.R