Food insecurity maps: season models
Maps of the prediction of the latent abilities and the spatial component.
Load required libraries and data
rm(list = ls())
library(day2day)
library(spifa)
library(xtable)
library(tibble)
library(sf)#> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUElibrary(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, unionlibrary(tidyr)
library(purrr)
library(ggplot2)
library(patchwork)
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_fig <- file.path(path_data, "figures")
path_src <- file.path(path_main, "src")
source(file.path(path_src, "itemnames.R"))
source(file.path(path_src, "ggtheme-publication.R"))
source(file.path(path_src, "maps.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 84ipi_osm <- readRDS(file.path(path_processed, "ipixuna-osm.rds"))
pred_spifa <- readRDS(file.path(path_modelled, "spifa-prediction-ipixuna-urban-season-3gp.rds"))
pred_spifa_sp <- readRDS(file.path(path_modelled, "spifa-predictionsp-ipixuna-urban-season-3gp.rds"))Get summary
# custom function
fun_sum <- function(data, regex, fun){
dplyr::select(data, matches(regex)) |>
st_drop_geometry() |>
apply(1, fun)
}
# compute summaries
pred_summary <- pred_spifa["id"]
pred_summary[paste0("F", 1:3, "-exceedance-dry")] <- lapply(c("^F1-.+dry$", "^F2-.+dry$", "^F3-.+dry$"),
function(z) fun_sum(pred_spifa, z, function(x) sum(x > 0) / length(x)))
pred_summary[paste0("F", 1:3, "-exceedance-wet")] <- lapply(c("^F1-.+wet$", "^F2-.+wet$", "^F3-.+wet$"),
function(z) fun_sum(pred_spifa, z, function(x) sum(x > 0) / length(x)))
pred_summary[paste0("F", 1:3, "-median-dry")] <- lapply(c("^F1-.+dry$", "^F2-.+dry$", "^F3-.+dry$"),
function(z) fun_sum(pred_spifa, z, median))
pred_summary[paste0("F", 1:3, "-median-wet")] <- lapply(c("^F1-.+wet$", "^F2-.+wet$", "^F3-.+wet$"),
function(z) fun_sum(pred_spifa, z, median))
# compute summaries
pred_summary_sp <- pred_spifa_sp["id"]
pred_summary_sp[paste0("F", 1:3, "-exceedance-dry")] <- lapply(c("^F1-.+dry$", "^F2-.+dry$", "^F3-.+dry$"),
function(z) fun_sum(pred_spifa_sp, z, function(x) sum(x > 0) / length(x)))
pred_summary_sp[paste0("F", 1:3, "-exceedance-wet")] <- lapply(c("^F1-.+wet$", "^F2-.+wet$", "^F3-.+wet$"),
function(z) fun_sum(pred_spifa_sp, z, function(x) sum(x > 0) / length(x)))
pred_summary_sp[paste0("F", 1:3, "-median-dry")] <- lapply(c("^F1-.+dry$", "^F2-.+dry$", "^F3-.+dry$"),
function(z) fun_sum(pred_spifa_sp, z, median))
pred_summary_sp[paste0("F", 1:3, "-median-wet")] <- lapply(c("^F1-.+wet$", "^F2-.+wet$", "^F3-.+wet$"),
function(z) fun_sum(pred_spifa_sp, z, median))Detect high insecure areas
mylabs <- c(
"F1-exceedance" = "atop('(F1) Severe food insecurity', 'including children')",
"F2-exceedance" = "atop('(F2) Anxiety and relying on', 'social relations')",
"F3-exceedance" = "atop('(F3) Adults eating less and', 'going hungry')"
)
myseason <- c("dry" = "Dry~season", "wet" = "Rainy~season")
pred_summary_sp |>
dplyr::select(id, matches("exceedance")) |>
tidyr::pivot_longer(matches("exceedance"),
names_to = c("factor_exceed", "season"), names_pattern = "(F[1-3]-exceedance)-(.+)") |>
mutate(factor_exceed = mylabs[factor_exceed], season = myseason[season]) |>
mutate(value = ifelse(value >= 0.75, 1, NA)) |>
fi_ipi_exceed(osm = ipi_osm, title = "title", license = TRUE) +
facet_grid(season ~ factor_exceed, labeller = ggplot2::label_parsed) +
theme_map_Publication2(base_size = 9) +
theme(legend.position = "none")#> Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
#> ℹ Please use the `linewidth` argument instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
Text to add
text_df <- tibble(
# text = c(toupper(letters[1:10])),
text = c("C2", "D2", "B", "A", "C", "F2", "G2", "D", "I2", "J2", "E", "F", "G"),
x1 = c(-71.695, -71.69, -71.6975, -71.6837, -71.6915, -71.686, -71.698, -71.694,
-71.6875, -71.6875, -71.6849, -71.699, -71.6915),
y1 = c(-7.045, -7.045, -7.051, -7.0585, -7.0488, -7.048, -7.057, -7.045, -7.0452,
-7.0558, -7.0572, -7.0585, -7.0483),
x2 = c(-71.6975, -71.6845, -71.6993, -71.6817, -71.699, -71.683, -71.699, -71.6975,
-71.684, -71.682, -71.683, -71.696, -71.699),
y2 = c(-7.045, -7.044, -7.048, -7.055, -7.049, -7.049, -7.06, -7.045, -7.045,
-7.055, -7.0535, -7.0565, -7.0483),
factor = c(rep("F1-median", 4), rep("F2-median", 3), rep("F3-median", 3), rep("F1-median", 2),
"F3-median"),
factor_exceed = c(rep(mylabs[[1]], 4), rep(mylabs[[2]], 3),
rep(mylabs[[3]], 3), rep(mylabs[[1]], 2), mylabs[[3]]),
season = myseason[c(rep(NA, 2), "wet", "dry", rep(NA, 6), rep("wet", 3))],
) |>
dplyr::filter(text %in% toupper(letters[1:7])) |>
st_as_sf(coords = c("x1", "y1"), crs = st_crs(4326)) |>
st_transform(crs = st_crs(3857))
text_df[c("x1", "y1")] <- st_coordinates(text_df)
text_df <- text_df |>
st_set_geometry(NULL) |>
st_as_sf(coords = c("x2", "y2"), crs = st_crs(4326)) |>
st_transform(crs = st_crs(3857))
text_df[c("x2", "y2")] <- st_coordinates(text_df)
text_df <- text_df |> st_set_geometry(NULL)
text_df <- subset(text_df, !is.na(season))
text_df#> # A tibble: 5 × 8
#> text factor factor_exceed season x1 y1 x2 y2
#> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 B F1-median atop('(F1) Severe food insecurity', 'inclu… Rainy… -7.98e6 -7.87e5 -7.98e6 -7.87e5
#> 2 A F1-median atop('(F1) Severe food insecurity', 'inclu… Dry~s… -7.98e6 -7.88e5 -7.98e6 -7.87e5
#> 3 E F1-median atop('(F1) Severe food insecurity', 'inclu… Rainy… -7.98e6 -7.88e5 -7.98e6 -7.87e5
#> 4 F F1-median atop('(F1) Severe food insecurity', 'inclu… Rainy… -7.98e6 -7.88e5 -7.98e6 -7.88e5
#> 5 G F3-median atop('(F3) Adults eating less and', 'going… Rainy… -7.98e6 -7.87e5 -7.98e6 -7.87e5text_dry <- subset(text_df, season == myseason["dry"])
text_wet <- subset(text_df, season == myseason["wet"])Visualize prediction of latent abilities
gg1 <- fi_ipi(mutate(pred_summary, `F1-median` = `F1-median-dry`), `F1-median`, ipi_osm,
title = mylabs[[1]], axis.x = FALSE, text_df = text_dry) +
fi_ipi(mutate(pred_summary, `F2-median` = `F2-median-dry`), `F2-median`, ipi_osm,
title = mylabs[[2]], axis.y = FALSE, axis.x = FALSE, text_df = text_dry) +
fi_ipi(mutate(pred_summary, `F3-median` = `F3-median-dry`), `F3-median`, ipi_osm,
title = mylabs[[3]], axis.y = FALSE, axis.x = FALSE, text_df = text_dry,
license = FALSE, title2 = "Dry~season") &
theme_map_Publication(base_size = 9) &
theme(
strip.text.x = element_text(margin=margin(l=0, r=0, t=0, b=2)),
strip.text.y = element_text(margin=margin(l=2, r=0, t=0, b=0))
)#> # A tibble: 1 × 8
#> text factor factor_exceed season x1 y1 x2 y2
#> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 A F1-median atop('(F1) Severe food insecurity', 'inclu… Dry~s… -7.98e6 -7.88e5 -7.98e6 -7.87e5#> Warning: `aes_()` was deprecated in ggplot2 3.0.0.
#> ℹ Please use tidy evaluation idioms with `aes()`
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.#> # A tibble: 0 × 8
#> # ℹ 8 variables: text <chr>, factor <chr>, factor_exceed <chr>, season <chr>, x1 <dbl>, y1 <dbl>,
#> # x2 <dbl>, y2 <dbl>
#> # A tibble: 0 × 8
#> # ℹ 8 variables: text <chr>, factor <chr>, factor_exceed <chr>, season <chr>, x1 <dbl>, y1 <dbl>,
#> # x2 <dbl>, y2 <dbl>gg2 <- fi_ipi(mutate(pred_summary, `F1-median` = `F1-median-wet`), `F1-median`, ipi_osm,
title = mylabs[[1]], text_df = text_wet) +
fi_ipi(mutate(pred_summary, `F2-median` = `F2-median-wet`), `F2-median`, ipi_osm,
title = mylabs[[2]], axis.y = FALSE, text_df = text_wet) +
fi_ipi(mutate(pred_summary, `F3-median` = `F3-median-wet`), `F3-median`, ipi_osm,
title = mylabs[[3]], axis.y = FALSE, text_df = text_wet, license = TRUE, title2 = "Rainy~season") &
theme_map_Publication(base_size = 9) &
theme(strip.background = element_blank(), strip.text.x = element_blank(),
strip.text.y = element_text(margin=margin(l=2, r=0, t=0, b=0)))#> # A tibble: 3 × 8
#> text factor factor_exceed season x1 y1 x2 y2
#> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 B F1-median atop('(F1) Severe food insecurity', 'inclu… Rainy… -7.98e6 -7.87e5 -7.98e6 -7.87e5
#> 2 E F1-median atop('(F1) Severe food insecurity', 'inclu… Rainy… -7.98e6 -7.88e5 -7.98e6 -7.87e5
#> 3 F F1-median atop('(F1) Severe food insecurity', 'inclu… Rainy… -7.98e6 -7.88e5 -7.98e6 -7.88e5
#> # A tibble: 0 × 8
#> # ℹ 8 variables: text <chr>, factor <chr>, factor_exceed <chr>, season <chr>, x1 <dbl>, y1 <dbl>,
#> # x2 <dbl>, y2 <dbl>
#> # A tibble: 1 × 8
#> text factor factor_exceed season x1 y1 x2 y2
#> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 G F3-median atop('(F3) Adults eating less and', 'going… Rainy… -7.98e6 -7.87e5 -7.98e6 -7.87e5(gg1 / gg2) & theme(legend.text = element_text(size = rel(0.7)))
ggsave(file.path(path_fig, "osmmap-season-theta-median.pdf"), width = 7, height = 3.8 * 1.7)
myseason <- c("dry" = "Dry~season", "wet" = "Rainy~season")
pred_summary |>
dplyr::select(id, matches("exceedance")) |>
tidyr::pivot_longer(matches("exceedance"),
names_to = c("factor_exceed", "season"), names_pattern = "(F[1-3]-exceedance)-(.+)") |>
mutate(factor_exceed = mylabs[factor_exceed], season = myseason[season]) |>
within({value[value < 0.5] <- NA}) |>
fi_ipi_exceed(osm = ipi_osm, title = "title", text_df = text_df, license = TRUE) +
facet_grid(season ~ factor_exceed, labeller = ggplot2::label_parsed) +
theme_map_Publication2(base_size = 9)
ggsave(file.path(path_fig, "osmmap-season-theta-exceedance.pdf"), width = 7, height = 3.7 * 1.7)Visualize prediction of the spatial component of latent abilities
gg1 <- fi_ipi(mutate(pred_summary_sp, `F1-median` = `F1-median-dry`), `F1-median`, ipi_osm,
title = mylabs[[1]], axis.x = FALSE, text_df = text_dry) +
fi_ipi(mutate(pred_summary_sp, `F2-median` = `F2-median-dry`), `F2-median`, ipi_osm,
title = mylabs[[2]], axis.y = FALSE, axis.x = FALSE, text_df = text_dry) +
fi_ipi(mutate(pred_summary_sp, `F3-median` = `F3-median-dry`), `F3-median`, ipi_osm,
title = mylabs[[3]], axis.y = FALSE, axis.x = FALSE, text_df = text_dry,
license = FALSE, title2 = "Dry~season") &
theme_map_Publication(base_size = 9) &
theme(
strip.text.x = element_text(margin=margin(l=0, r=0, t=0, b=2)),
strip.text.y = element_text(margin=margin(l=2, r=0, t=0, b=0))
)#> # A tibble: 1 × 8
#> text factor factor_exceed season x1 y1 x2 y2
#> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 A F1-median atop('(F1) Severe food insecurity', 'inclu… Dry~s… -7.98e6 -7.88e5 -7.98e6 -7.87e5
#> # A tibble: 0 × 8
#> # ℹ 8 variables: text <chr>, factor <chr>, factor_exceed <chr>, season <chr>, x1 <dbl>, y1 <dbl>,
#> # x2 <dbl>, y2 <dbl>
#> # A tibble: 0 × 8
#> # ℹ 8 variables: text <chr>, factor <chr>, factor_exceed <chr>, season <chr>, x1 <dbl>, y1 <dbl>,
#> # x2 <dbl>, y2 <dbl>gg2 <- fi_ipi(mutate(pred_summary_sp, `F1-median` = `F1-median-wet`), `F1-median`, ipi_osm,
title = mylabs[[1]], text_df = text_wet) +
fi_ipi(mutate(pred_summary_sp, `F2-median` = `F2-median-wet`), `F2-median`, ipi_osm,
title = mylabs[[2]], axis.y = FALSE, text_df = text_wet) +
fi_ipi(mutate(pred_summary_sp, `F3-median` = `F3-median-wet`), `F3-median`, ipi_osm,
title = mylabs[[3]], axis.y = FALSE, text_df = text_wet, license = TRUE, title2 = "Rainy~season") &
theme_map_Publication(base_size = 9) &
theme(strip.background = element_blank(), strip.text.x = element_blank(),
strip.text.y = element_text(margin=margin(l=2, r=0, t=0, b=0)))#> # A tibble: 3 × 8
#> text factor factor_exceed season x1 y1 x2 y2
#> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 B F1-median atop('(F1) Severe food insecurity', 'inclu… Rainy… -7.98e6 -7.87e5 -7.98e6 -7.87e5
#> 2 E F1-median atop('(F1) Severe food insecurity', 'inclu… Rainy… -7.98e6 -7.88e5 -7.98e6 -7.87e5
#> 3 F F1-median atop('(F1) Severe food insecurity', 'inclu… Rainy… -7.98e6 -7.88e5 -7.98e6 -7.88e5
#> # A tibble: 0 × 8
#> # ℹ 8 variables: text <chr>, factor <chr>, factor_exceed <chr>, season <chr>, x1 <dbl>, y1 <dbl>,
#> # x2 <dbl>, y2 <dbl>
#> # A tibble: 1 × 8
#> text factor factor_exceed season x1 y1 x2 y2
#> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 G F3-median atop('(F3) Adults eating less and', 'going… Rainy… -7.98e6 -7.87e5 -7.98e6 -7.87e5(gg1 / gg2) & theme(legend.text = element_text(size = rel(0.7)))
ggsave(file.path(path_fig, "osmmap-season-thetaspatial-median.pdf"), width = 7, height = 3.8 * 1.7)
myseason <- c("dry" = "Dry~season", "wet" = "Rainy~season")
pred_summary_sp |>
dplyr::select(id, matches("exceedance")) |>
tidyr::pivot_longer(matches("exceedance"),
names_to = c("factor_exceed", "season"), names_pattern = "(F[1-3]-exceedance)-(.+)") |>
mutate(factor_exceed = mylabs[factor_exceed], season = myseason[season]) |>
within({value[value < 0.5] <- NA}) |>
fi_ipi_exceed(osm = ipi_osm, title = "title", text_df = text_df, license = TRUE) +
facet_grid(season ~ factor_exceed, labeller = ggplot2::label_parsed) +
theme_map_Publication2(base_size = 9)
ggsave(file.path(path_fig, "osmmap-season-thetaspatial-exceedance.pdf"), width = 7, height = 3.7 * 1.7)Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 31.897 2.569 31.923