NA values in predictions from gstat::krige()

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

NA values in predictions from gstat::krige()

robinlovelace
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]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: NA values in predictions from gstat::krige()

robinlovelace
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 <[hidden email]> 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
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo