[GIS] Kriging example with sf object

gstatkrigingrsfsp

I am new to geo spacial analysis. I would like to learn kriging from a simple example. Say, we have four points with some z value.

library(sf)
library(gstat)
library(tidyverse)


dt<-tibble::tribble(
  ~id,  ~lon,  ~lat,    ~z,
  "A",   500,   500,  12,
  "B",  1000,   500,  13,
  "C",   500,  1000,  15,
  "D",  1000,  1000,  17
  )

Here I create an sf object and visualize it

DT_sf  <-  st_as_sf(dt, coords = c("lon", "lat"), 
                 crs = 4326, agr = "constant")

ggplot() +
  geom_sf(data=DT_sf, aes(color=z), size=10)

How can I interpolate the z values, some average measurements onto the whole area?

So far I have explored I need to prepare a grid: like this (?) and a matrix of distances

grd <- st_sf(geom=st_make_grid(DT_sf), crs=4326)

dist<-spDists(as.matrix(dt[2:3]), longlat = TRUE)

coef = lm(log(z)~sqrt(dist), dt)$coef

I tried this and it is probably nonsense (i do not know what I am doing).

k = krige(log(z)~dist, as_Spatial(DT_sf), as_Spatial(grd), vgm(.6, "Sph", 900), beta = coef)

What package and command would result in extrapolation of z values onto the whole surface? What parameters would be needed?

The krige function now returns error:

Error in gstat.formula.predict(d$formula, newdata, na.action = na.action,  : 
  NROW(locs) != NROW(X): this should not occur
In addition: Warning messages:
1: 'newdata' had 100 rows but variables found have 4 rows 
2: 'newdata' had 100 rows but variables found have 4 rows

Best Answer

After the study of https://rpubs.com/nabilabd/118172 and https://rpubs.com/nabilabd/134781 I assembled the following solution.

library(sf)
library(sp)
library(gstat)
library(tidyverse)

Here are points I would like to play with.

dt<-tibble::tribble(
  ~id,  ~lon,  ~lat,    ~z,
  "a",  1.1,   1.3,  12.5,
  "b",  2.4,   1.7,  13.0,
  "c",  3.2,   1.4,  12.0,
  "d",  4.2,   1.4,  16,
  "e",  1.2,   2.3,  13.0,
  "f",  2.3,   2.7,  15.5,
  "g",  3.7,   2.5,  19.0,
  "h",  4.5,   2.2,  17.5,
  "i",  1.1,   3.2,  16.5,
  "j",  2.2,   3.4,  18.5,
  "k",  3.7,   3.3,  18.2,
  "l",  4.7,   3.3,  11.5,
  "m",  1.1,   4.1,  17.5,
  "n",  2.2,   4.2,  18.5,
  "o",  3.7,   4.3,  19.2,
  "p",  4.7,   4.2,  8  
)

DT_sf  <-  st_as_sf(dt, coords = c("lon", "lat"), 
                    crs = 4326, agr = "constant")

ggplot() +
  geom_sf(data=DT_sf, aes(color=z), size=10)


DT_sp  <- as_Spatial(DT_sf)

Now I need to create a grid of points I would like to predict the values for.

lon <- seq(1.0, 5.0, length.out = 100)
lat <- seq(1.0, 5.0, length.out = 100)
grd <- expand.grid(lon = lon, lat = lat)

grd_sf  <-  st_as_sf(grd, coords = c("lon", "lat"), 
                     crs = 4326, agr = "constant")

grd_sp <- as_Spatial(grd_sf)

Now I need to create a variogram. Like in lm values can be dependent on some feature like distance from the river (meuse dataset) and if there is no such variable, use 1.

dt.vgm <- variogram(z~1, DT_sp)

class(dt.vgm)

dt.fit <-
  fit.variogram(dt.vgm, model = vgm(1,"Lin",900,1)) # fit model

# vgm() list of models

plot(dt.vgm, dt.fit)

Now I can perform kriging and plot the result

lzn.kriged <- krige((z) ~ 1, DT_sp, grd_sp, model=dt.fit)

lzn.kriged %>% as.data.frame %>% rename(lon=coords.x1, lat=coords.x2) %>% 
  ggplot(aes(x=lon, y=lat)) + geom_tile(aes(fill=var1.pred)) + coord_equal() +
  scale_fill_gradient2(low="green", mid = "yellow",  high="red",midpoint = 15) +
  theme_bw()+
  geom_point(data=dt, aes(color=z), size=10)+
  geom_text(data=dt, aes(label=z), color="white")
Related Question