# Kriging example with sf object

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
1: 'newdata' had 100 rows but variables found have 4 rows
2: 'newdata' had 100 rows but variables found have 4 rows
``````
• Does that `krige` function call work or produce an error message? – Spacedman Jun 30 '18 at 16:52
• Currently krige function returns error, probably because the data is incompatible with what it expects (i pasted the error message above). – Jacek Kotowski Jun 30 '18 at 19:13
• Have you understood the example in the help for `krige`? Have you read an introduction to kriging and geostatistics? Probably better you do that than we duplicate what kriging is here. – Spacedman Jul 1 '18 at 8:58
• @Spacedman I am starting now to understand the process. Pasted below. – Jacek Kotowski Jul 3 '18 at 13:22
• @JacekKotowski maybe this can help gis.stackexchange.com/questions/239301/… – Guzmán Jul 3 '18 at 14:32

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")
``````