SPIFA Urban Ipixuna Season: spatial prediction

Load required libraries and data

rm(list = ls())
library(day2day)
library(dplyr)

#> 
#> Attaching package: 'dplyr'

#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag

#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

library(mirt)

#> Loading required package: stats4

#> Loading required package: lattice

library(spifa)
library(tidyr)
library(sf)

#> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE

library(purrr)
library(ggplot2)

path_main <- git_path()
path_data <- file.path(path_main, "data")
path_raw <- file.path(path_data, "raw")
path_processed <- file.path(path_data, "processed")
path_modelled <- file.path(path_data, "modelled")
path_src <- file.path(path_main, "src")

devtools::load_all("/home/rstudio/documents/repositories/spifa")

#> ℹ Loading spifa

source(file.path(path_src, "prediction.R"))

fidata <- file.path(path_processed, "fi-items-ipixuna-urban.gpkg") |>
    st_read(as_tibble = TRUE)

#> Reading layer `fi-items-ipixuna-urban' from data source 
#>   `/home/rstudio/documents/projects/food-insecurity-mapping/data/processed/fi-items-ipixuna-urban.gpkg' 
#>   using driver `GPKG'
#> Simple feature collection with 200 features and 36 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -71.70038 ymin: -7.06058 xmax: -71.68109 ymax: -7.03724
#> Geodetic CRS:  WGS 84

fidata_3857 <- sf::st_transform(fidata, crs = 3857)

samples_dry <- readRDS(file.path(path_modelled, "spifa-ipixuna-urban-dry-3gp.rds"))
samples_wet <- readRDS(file.path(path_modelled, "spifa-ipixuna-urban-wet-3gp.rds"))
iter <- nrow(samples_dry[["theta"]])

Create grid for prediction

# create area of interest for prediction
buffer <- sf::st_geometry(fidata_3857) %>%
  sf::st_union() %>%
  sf::st_convex_hull() %>%
  sf::st_buffer(50) %>%
  sf::st_sf()

# craselect only points inside the area of prediction
fipred <- st_sf(geom = st_make_grid(buffer, cellsize = c(30,30))) |>
    sf::st_join(buffer, join = st_intersects, left = FALSE)

# visualize
ggplot(fipred) +
    geom_sf(data = buffer) +
    geom_sf(fill = NA, color = 4) +
    geom_sf(data = fidata_3857) +
    theme_bw()

Perform spatial prediction dry

set.seed(123)
thin <- 6

# predict
system.time(
predicted_wet <- as_tibble(samples_wet, burnin = 1 * iter / 5, thin = thin) %>%
        spifa:::predict.spifa(newcoords = st_coordinates(st_centroid(fipred)))
)

#>     user   system  elapsed 
#> 12937.39  5851.08 12732.36

system.time(
predicted_dry <- as_tibble(samples_dry, burnin = 1 * iter / 5, thin = thin) %>%
        spifa:::predict.spifa(newcoords = st_coordinates(st_centroid(fipred)))
)

#>      user    system   elapsed 
#> 12628.583  5634.133 12393.434

# prediction as sf
fipred_wet  <- pred2df(predicted_wet, fipred) |>
    rename_with(function(x) paste0(x, "-wet"), matches("^F[1-3]-"))

#> Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if `.name_repair` is omitted as of tibble 2.0.0.
#> ℹ Using compatibility `.name_repair`.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

fipred_dry  <- pred2df(predicted_dry, fipred) |>
    rename_with(function(x) paste0(x, "-dry"), matches("^F[1-3]-")) |>
    st_drop_geometry()
fipred <- left_join(fipred_wet, fipred_dry, by = join_by(id))

# visualize means
fipred |>
    mutate(`F1-dry-mean` = rowMeans(dplyr::select(as_tibble(fipred), matches("^F1-.+dry$")))) |>
    mutate(`F2-dry-mean` = rowMeans(dplyr::select(as_tibble(fipred), matches("^F2-.+dry$")))) |>
    mutate(`F3-dry-mean` = rowMeans(dplyr::select(as_tibble(fipred), matches("^F3-.+dry$")))) |>
    mutate(`F1-wet-mean` = rowMeans(dplyr::select(as_tibble(fipred), matches("^F1-.+wet$")))) |>
    mutate(`F2-wet-mean` = rowMeans(dplyr::select(as_tibble(fipred), matches("^F2-.+wet$")))) |>
    mutate(`F3-wet-mean` = rowMeans(dplyr::select(as_tibble(fipred), matches("^F3-.+wet$")))) |>
    dplyr::select(matches("mean$")) |>
    pivot_longer(matches("mean$")) |>
    mutate(name = factor(name, levels = c("F1-dry-mean", "F2-dry-mean", "F3-dry-mean",
        "F1-wet-mean", "F2-wet-mean", "F3-wet-mean"))) |>
    ggplot() +
        geom_sf(aes(fill = value), color = NA) +
        scale_fill_distiller(palette = "RdBu") +
        facet_wrap(~ name)

Save predicted values

saveRDS(fipred, file = file.path(path_modelled, "spifa-prediction-ipixuna-urban-season-3gp.rds"))

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>     user   system  elapsed 
#> 25607.91 11489.94 25201.74