SPIFA Urban Ipixuna: 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)

eifa_data <- readRDS(file.path(path_modelled, "eifa-ipixuna-urban.rds"))

samples <- readRDS(file.path(path_modelled, "spifa-ipixuna-urban-3gp.rds"))
iter <- nrow(samples[["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

set.seed(123)
thin <- 6

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

#>      user    system   elapsed 
#> 15847.114  6435.024 13298.966

# prediction as sf
fipred <- pred2df(predicted, fipred)

#> 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.

# visualize means
fipred |>
    mutate(`F1-mean` = rowMeans(dplyr::select(as_tibble(fipred), matches("^F1-")))) |>
    mutate(`F2-mean` = rowMeans(dplyr::select(as_tibble(fipred), matches("^F2-")))) |>
    mutate(`F3-mean` = rowMeans(dplyr::select(as_tibble(fipred), matches("^F3-")))) |>
    dplyr::select(matches("mean$")) |>
    pivot_longer(matches("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-3gp.rds"))

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>      user    system   elapsed 
#> 15875.014  6438.253 13351.754