9

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 
5
  • Does that krige function call work or produce an error message?
    – Spacedman
    Jun 30, 2018 at 16:52
  • Currently krige function returns error, probably because the data is incompatible with what it expects (i pasted the error message above). Jun 30, 2018 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, 2018 at 8:58
  • @Spacedman I am starting now to understand the process. Pasted below. Jul 3, 2018 at 13:22
  • 1
    @JacekKotowski maybe this can help gis.stackexchange.com/questions/239301/…
    – Guz
    Jul 3, 2018 at 14:32

1 Answer 1

7

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 and acknowledge you have read our privacy policy.

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