2

The gstat package supports stratified ordinary kriging using sp objects as shown at https://github.com/r-spatial/gstat/blob/master/demo/examples.R (the code below). This year (https://r-spatial.github.io/gstat/news/index.html) gstat also got support for sf and stars objects. While sf and stars works well for many geostatistical methods, I am unable to make it work for stratified ordinary kriging.

# an example of gstat + sp
library(gstat)
library(sp)
data(meuse)
coordinates(meuse) = ~x+y
data(meuse.grid)
gridded(meuse.grid) = ~x+y

meuse$part.a = krige(part.a ~ 1, meuse.grid, meuse, nmax = 1)$var1.pred
#> [inverse distance weighted interpolation]
x1 = krige(
        log(zinc) ~ 1,
        meuse[meuse$part.a == 0, ],
        meuse.grid[meuse.grid$part.a == 0, ],
        model = vgm(.548, "Sph", 900, .0654),
        nmin = 20,
        nmax = 40,
        maxdist = 1000
)
#> [using ordinary kriging]
x2 = krige(
        log(zinc) ~ 1,
        meuse[meuse$part.a == 1, ],
        meuse.grid[meuse.grid$part.a == 1, ],
        model = vgm(.716, "Sph", 900),
        nmin = 20,
        nmax = 40,
        maxdist = 1000
)
#> [using ordinary kriging]
x = rbind(as.data.frame(x1), as.data.frame(x2))
gridded(x) = ~x+y
spplot(x["var1.pred"], main = "stratified kriging predictions")

The main problem is how to use only selected values from a stars object.

library(stars)
#> Loading required package: abind
#> Loading required package: sf
#> Linking to GEOS 3.7.1, GDAL 2.3.2, PROJ 5.2.0
library(sp)
data(meuse.grid)
gridded(meuse.grid) = ~x+y
meuse.grid_stars = st_as_stars(meuse.grid)
meuse.grid_stars[meuse.grid_stars$part.a == 1, ]
#> Warning in structure(eval(rlang::expr(x[[i]][!!!args])), levels = attr(x[[i]], : Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
#>   Consider 'structure(list(), *)' instead.

Is it possible to do stratified ordinary kriging using gstat and stars?

1
  • I don't know if it may help, but here a useful (and pretty recent approach) usign sf and stars. Did you find a way to solve that error?
    – Nico
    Commented Dec 7, 2020 at 10:19

1 Answer 1

0

I found similar question here. And the answer was published in github.

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.