D.4 Predicción espacial (KU)

Para generar la rejilla de predicción consideramos un buffer de radio 40 en torno a las posiciones espaciales:

buffer <- aquifer_sf %>% st_geometry() %>% st_buffer(40)

En lugar de emplear una rejilla sf:

# grid <- buffer %>% st_make_grid(n = c(50, 50), what = "centers") %>% st_intersection(buffer)

por comodidad es preferible emplear una rejilla stars:

library(stars)
grid <- buffer %>%  st_as_stars(nx = 50, ny = 50)

Si suponemos un modelo (no constante) para la tendencia, es necesario añadir los valores de las variables explicativas a la rejilla de predicción:

coord <- st_coordinates(grid)
grid$lon <- coord$x
grid$lat <- coord$y

Además, en este caso recortamos la rejilla para filtrar predicciones alejadas de las observaciones:

grid <- grid %>% st_crop(buffer)

Obtenemos las predicciones mediante kriging universal (Sección 4.3 y Sección 4.4):

pred <- krige(formula = head ~ lon + lat, locations = aquifer_sf, model = fit,
              newdata = grid)

## [using universal kriging]

ERROR en krige: cambia las coordenadas del objeto stars

summary(st_coordinates(grid))

##        x                 y         
##  Min.   :-181.86   Min.   :-28.03  
##  1st Qu.:-100.73   1st Qu.: 33.25  
##  Median : -16.22   Median : 97.09  
##  Mean   : -16.22   Mean   : 97.09  
##  3rd Qu.:  68.29   3rd Qu.:160.93  
##  Max.   : 149.42   Max.   :222.21

summary(st_coordinates(pred))

##        x                 y         
##  Min.   :-181.86   Min.   :-28.03  
##  1st Qu.:-100.73   1st Qu.: 33.25  
##  Median : -16.22   Median : 97.09  
##  Mean   : -16.22   Mean   : 97.09  
##  3rd Qu.:  68.29   3rd Qu.:160.93  
##  Max.   : 149.42   Max.   :222.21

Posible solución: añadir el resultado a grid:

grid$var1.pred <- pred$var1.pred
grid$var1.var <- pred$var1.var

Finalmente representamos las predicciones y las varianzas kriging:

plot(grid["var1.pred"], breaks = "equal", col = sf.colors(64), key.pos = 4,
     main = "Predicciones kriging")

plot(grid["var1.var"], breaks = "equal", col = sf.colors(64), key.pos = 4,
     main = "Varianzas kriging")

También podríamos emplear el paquete ggplot2:

library(ggplot2)
library(gridExtra)
p1 <- ggplot() + geom_stars(data = grid, aes(fill = var1.pred, x = x, y = y)) +
    scale_fill_viridis_c() + geom_sf(data = aquifer_sf) +
    coord_sf(lims_method = "geometry_bbox")
p2 <- ggplot() + geom_stars(data = grid, aes(fill = var1.var, x = x, y = y)) +
    scale_fill_viridis_c() + geom_sf(data = aquifer_sf) +
    coord_sf(lims_method = "geometry_bbox")
grid.arrange(p1, p2, ncol = 2)