[R-sig-Geo] NA values in predictions from gstat::krige()

[Robin Lovelace]

Hi all,

When using the gstat package I get NA values for the majority of cells
predicted and I'm not sure what I'm doing wrong. The data is quite noisy
and I plan to do some aggregation on the input data before trying this
again with the interpolate() function from raster/terra, but thought it
worth asking as the answer may be of use to others.

# reprex
library(sf)
library(gstat)
library(stars)
u = "
https://github.com/saferactive/saferactive/releases/download/0.1.1/rnet_lnd_1pcnt.Rds
"
rnet_lnd = readRDS(url(u))
rnet_lnd
grd = st_bbox(rnet_lnd) %>%
  st_as_stars(dx = 500, dy = 500) %>%
  st_set_crs(27700) %>%
  st_crop(rnet_lnd)
grd
v = variogram(bicycle~1, rnet_lnd, cutoff = 10000)
plot(v)
vm = fit.variogram(v, vgm(psill = "Sph", model = "Exp", range = 10000,
nugget = 1))
plot(vm, cutoff = 10000)
rnet_krige = gstat::krige(bicycle~1, rnet_lnd, grd, vm, nmax = 100)
plot(rnet_lnd$geometry)
plot(rnet_krige, add = TRUE)

Outcome on my computer: predictions only in raster grid cells with
observations but I expected from the Spatial Data Science book a continuous
surface with predictions for all of the grid cells, as shown here:
https://keen-swartz-3146c4.netlify.app/interpolation.html

I've also asked this question on stack-overflow where the output from the
above commands can be found: https://stackoverflow.com/questions/64541951/

Thanks for reading and apologies if I'm missing something obvious.

Robin

	[[alternative HTML version deleted]]