Update, I had missed something obvious: the unnecessary st_crop. Many thanks to Jannes, who pointed this out: https://github.com/Robinlovelace/geocompr/issues/585#issuecomment-716780818 And apologies for 'spamming' the list, hope my next comment/question will be more enlightening! Robin
On Mon, Oct 26, 2020 at 5:31 PM Robin Lovelace <rob00x at gmail.com> wrote:
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