4

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 
  • 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
  • 1
    @JacekKotowski maybe this can help gis.stackexchange.com/questions/239301/… – Guzmán Jul 3 '18 at 14:32
3

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.