Merge pull request #14 from IDEELResearch/feature/covars

Feature/covars

.gitignore

@@ -13,3 +13,4 @@ analysis/data-dhs
 analysis/data-derived/sims
 wc2-5
 *deprecated*
+.DS_Store

DESCRIPTION

@@ -21,13 +21,16 @@ Imports:
     terra,
     dplyr,
     sf,
-    raster
+    purrr,
+    rlang,
+    scales
 Suggests:
     knitr,
     rmarkdown,
     testthat
 Remotes:
-    mrc-ide/STAVE
+    mrc-ide/STAVE,
+    malaria-atlas-project/malariaAtlas
 URL: https://IDEELResearch.github.io/stitch/
 Encoding: UTF-8
 LazyData: true

NAMESPACE

@@ -1,6 +1,33 @@
 # Generated by roxygen2: do not edit by hand
 
 export("%>%")
+importFrom(dplyr,filter)
+importFrom(dplyr,left_join)
+importFrom(dplyr,mutate)
+importFrom(dplyr,relocate)
+importFrom(dplyr,rename)
+importFrom(ggplot2,aes)
+importFrom(ggplot2,coord_sf)
+importFrom(ggplot2,element_rect)
+importFrom(ggplot2,facet_wrap)
+importFrom(ggplot2,geom_sf)
+importFrom(ggplot2,ggplot)
+importFrom(ggplot2,ggsave)
+importFrom(ggplot2,scale_fill_viridis_c)
+importFrom(ggplot2,theme)
+importFrom(ggplot2,theme_void)
 importFrom(grDevices,dev.off)
 importFrom(grDevices,pdf)
 importFrom(magrittr,"%>%")
+importFrom(purrr,imap_dfr)
+importFrom(rlang,.data)
+importFrom(rlang,sym)
+importFrom(scales,percent_format)
+importFrom(sf,sf_use_s2)
+importFrom(sf,st_centroid)
+importFrom(sf,st_coordinates)
+importFrom(sf,st_geometry)
+importFrom(stringr,str_glue)
+importFrom(terra,ext)
+importFrom(terra,rast)
+importFrom(terra,rasterize)

R/append_lat_lon.R

@@ -0,0 +1,34 @@
+#' Get 'centroid' coordinates in lat and lon
+#'
+#' This function extracts the centroid coordinates from the geometries of a shape file object and appends them as columns
+#'
+#' @param shp The sf object
+#' @return A shape file with lat lon 'centroid' coordinates appended as columns before the geometries
+#' @importFrom sf sf_use_s2 st_geometry st_centroid st_coordinates
+#' @importFrom dplyr rename relocate
+#' @importFrom rlang .data
+#'
+#' @examples
+#' # shp_adm0 <- readRDS("analysis/data_derived/sf_admin0_africa.rds")
+#' # shp_adm0 <- shp_adm0 %>% append_lat_lon()
+
+append_lat_lon <- function(shp){
+
+  # switching off spherical geometry
+  sf_use_s2(F)
+
+  shp <- shp %>%
+    # get geometry col
+    st_geometry() %>%
+    # get centroid of polygon
+    st_centroid() %>%
+    # convert to lat/lon coords
+    st_coordinates() %>%
+    # bind with map_data
+    cbind(shp) %>%
+    rename("lon" = "X",
+           "lat" = "Y") %>%
+    relocate(c(.data$lat:.data$lon), .before = .data$geometry)
+
+  return(shp)
+}

R/create_raster_object.R

@@ -1,50 +1,47 @@
-# Function to convert a shapefile object to a raster object
-#
-# Args:
-#   shape_obj: An sf or SpatVector object that represents the geographic shapefile data
-#              (e.g., African countries or Uganda).
-#   res: Numeric value that represents the resolution of the raster. A smaller number gives higher
-#        resolution and more detail, but increases computational load (default is 0.05).
-#
-# Returns:
-#   A raster object created by rasterizing the input shapefile. Each cell in the raster corresponds
-#   to an area in the shapefile, and the values in the raster are based on the attribute 'iso'.
-#
-# Example usage:
-# # Read in Africa shape file
-# rds_file = "analysis/data_derived/sf_admin0_africa.rds"
-# africa_shp <- readRDS(file = rds_file)
-#
-# # Create a raster for Uganda
-# uganda_shp <- africa_shp %>% dplyr::filter(iso=="UGA")
-# uganda_raster <- create_raster_object(shape_obj=uganda_shp)
-# # Plot Uganda raster with Ugandan boarder
-# plot(uganda_raster)
-# plot(st_geometry(uganda_shp), add=TRUE, border="black")
-#
-# # Save the raster to a file
-# writeRaster(uganda_raster,
-#             filename="analysis/data_derived/uganda_raster_admin0.tif",
-#             filetype="GTiff",
-#             overwrite=TRUE)
-#
-# # Plot raster via ggplot
-# uganda_raster_df <- as.data.frame(uganda_raster, xy=TRUE, na.rm=TRUE)
-# ggplot() +
-#   geom_raster(data = uganda_raster_df, aes(x = x, y = y)) +
-#   #geom_sf(data=uganda_shp, color = "black", size = 0.5) +
-#   scale_fill_viridis_c() +  # Optional: Nice color scale
-#   coord_equal() +
-#   theme_minimal() +
-#   labs(
-#     title = "Africa Admin 0 Raster",
-#     x = "Longitude",
-#     y = "Latitude"
-#   )
-#
-# # Create African raster
-# africa_raster <- create_raster_object(shape_obj=africa_shp)
-# plot(africa_raster)
+#' Function to convert a shapefile object to a raster object
+#'
+#' This function converts a sf shape file into a raster object.
+#'
+#' @param shape_obj An sf or SpatVector object that represents the geographic shapefile data (e.g., African countries or Uganda).
+#' @param res Numeric value that represents the resolution of the raster (default is 0.05).
+#' @return A raster object created by rasterizing the input shapefile. Each cell in the raster corresponds to an area in the shapefile, and the values in the raster are based on the attribute 'iso'.
+#' @importFrom terra rast rasterize ext
+#' @examples
+#' # Example of usage:
+#' # Read in Africa shape file
+#' # rds_file = "analysis/data_derived/sf_admin0_africa.rds"
+#' # africa_shp <- readRDS(file = rds_file)
+#'
+#' # Create a raster for Uganda
+#' # uganda_shp <- africa_shp %>% dplyr::filter(iso=="UGA")
+#' # uganda_raster <- create_raster_object(shape_obj=uganda_shp)
+#' # Plot Uganda raster with Ugandan boarder
+#' # plot(uganda_raster)
+#' # plot(st_geometry(uganda_shp), add=TRUE, border="black")
+#'
+#' # Save the raster to a file
+#' # writeRaster(uganda_raster,
+#'             # filename="analysis/data_derived/uganda_raster_admin0.tif",
+#'             # filetype="GTiff",
+#'             # overwrite=TRUE)
+#'
+#' # Plot raster via ggplot
+#' # uganda_raster_df <- as.data.frame(uganda_raster, xy=TRUE, na.rm=TRUE)
+#' # ggplot() +
+#'   # geom_raster(data = uganda_raster_df, aes(x = x, y = y)) +
+#'   # #geom_sf(data=uganda_shp, color = "black", size = 0.5) +
+#'   # scale_fill_viridis_c() +  # Optional: Nice color scale
+#'   # coord_equal() +
+#'   # theme_minimal() +
+#'   # labs(
+#'     # title = "Africa Admin 0 Raster",
+#'     # x = "Longitude",
+#'     # y = "Latitude"
+#'   # )
+#'
+#' # Create African raster
+#' # africa_raster <- create_raster_object(shape_obj=africa_shp)
+#' # plot(africa_raster)
 
 create_raster_object <- function(shape_obj, res = 0.05) {
   # Create an empty raster grid covering the extent of the shape object with the specified resolution

R/get_map_covar_data.R

@@ -0,0 +1,31 @@
+#' Title: Get map and covariate data in long-format
+#'
+#' This unlists the covariate data frames and returns a sf object with the covariate data in long format.
+#' It expects an object with a list of `covars` covariate data frames (eg by year)
+#' and the `map` shp file from `scrape`. Bespoke function to wrangle `final_covariates.rds` from `scrape`.
+#'
+#' @param covars_obj Description of the parameter 'x'.
+#' @return A sf object with covariate data in long format from all covar list
+#' @importFrom purrr imap_dfr
+#' @importFrom dplyr mutate left_join
+#'
+#' @examples
+#' # Example of how to use the function
+#' # covars <- readRDS("analysis/data_raw/final_covariates.rds")
+#' # map_data <- get_map_covar_data(covars)
+
+get_map_covar_data <- function(covars_obj) {
+  # Unlist and combine ------------------------------------------------------
+  map <- covars_obj$map
+  covars <- covars_obj$covars
+
+  covars_df <- imap_dfr(covars, ~ {
+    year <- as.integer(sub("Y", "", .y))
+    .x %>% mutate(year = year)
+  })
+
+  # Join map geometries with covar data -------------------------------------
+  map_data <- map %>% left_join(covars_df)
+
+  return(map_data)
+}

R/plot_spatiotemporal_map.R

@@ -0,0 +1,45 @@
+#' Title: Plotting spatiotemporal maps (time facets) using long sf objects
+#'
+#' Plot spatiotemporal maps with time facets based on a vector of variable names to loop over.
+#' This expects a long format sf object (eg covariate data in long format returned from `get_map_covar_data()`)
+#' and a border shp sf object (this could be the same as the long sf if you want the same borders), alternatively,
+#' you could supply eg admin0 sf to plot only the country borders. Because rendering takes a lot of time, the function
+#' outputs to 'analysis/plots' directory. Will likely edit this to be more felxible (eg user-supplied titles, filenames, color palettes)
+#'
+#' @param long_sf A long-format sf object
+#' @param border_sf Shape files for the borders (eg adm0)
+#' @param variable_names A vector of variable names that you wish to plot, it will loop over the names if >1
+#' @return Description of the value returned by the function.
+#' @importFrom ggplot2 ggplot geom_sf scale_fill_viridis_c coord_sf facet_wrap theme_void ggsave theme aes element_rect
+#' @importFrom stringr str_glue
+#' @importFrom dplyr filter
+#' @importFrom rlang sym
+#' @importFrom scales percent_format
+#' @importFrom rlang .data
+#'
+#' @examples
+#' # Example of how to use the function
+#' # covars <- readRDS("analysis/data_raw/final_covariates.rds")
+#' # border_sf <- readRDS("analysis/data_derived/sf_admin0_africa.rds")
+#' # map_data <- get_map_covar_data(covars)
+#' # variable_names <- c("pfpr210_mean", "ft", "AL")
+#' # plot_spatiotemporal_map(map_data, border_sf, variable_names)
+
+plot_spatiotemporal_map <- function(long_sf, border_sf, variable_names){
+
+  for (v in variable_names) {
+    p <- long_sf %>%
+      filter(!is.na(.data[[v]])) %>%
+      ggplot() +
+      geom_sf(aes(fill = !!sym(v)), color = NA) +
+      geom_sf(data = border_sf, color = "black", fill = NA) +
+      scale_fill_viridis_c(option = "viridis", labels = scales::percent_format(accuracy = 1)) +
+      coord_sf() +
+      facet_wrap(~year) +
+      theme_void(base_size = 14) +
+      theme(plot.background = element_rect(fill = "white", color = "white"))
+
+    filename <- str_glue("analysis/plots/map_{v}.png")
+    ggsave(filename, p, width = 8, height = 6, dpi = 300)
+  }
+}

analysis/01_generating_base_maps.Rmd → analysis/01_generating_base_maps.R

@@ -1,65 +1,45 @@
----
-title: "Generating base maps"
-author: "Shazia Ruybal-Pesántez"
-date: "`r Sys.Date()`"
-output: 
-  html_document: 
-    code_folding: hide 
-    fig_width: 8
-    fig_height: 6
-    theme: cosmo
----
-
-```{r setup, include=FALSE}
-knitr::opts_chunk$set(echo = TRUE, warning = F, message = F)
+# title: "Generating base maps"
+# author: "Shazia Ruybal-Pesántez"
 
+# Load packages -----------------------------------------------------------
 # devtools::install_github("malaria-atlas-project/malariaAtlas")
 
 library(malariaAtlas)
-library(tidyverse)
-library(here)
-```
+library(dplyr)
+library(sf)
 
-Use `malariaAtlas` R package to get shape files for Admin 0 
-```{r adm0 shape files}
-adm0_shp_available <- malariaAtlas::listShp(admin_level = "admin0", printed = F) # 249 countries
+# Get shape files ---------------------------------------------------------
 
+# Use `malariaAtlas` R package to get shape files for Admin 0
+adm0_shp_available <- malariaAtlas::listShp(admin_level = "admin0", printed = F) # 249 countries
 shp_adm0 <- malariaAtlas::getShp(ISO = adm0_shp_available$iso, admin_level = "admin0")
-```
-
 
-Use `malariaAtlas` R package to get shape files for Admin 0 
-```{r adm1 shape files}
 adm1_shp_available <- malariaAtlas::listShp(admin_level = "admin1", printed = F) # 2022 adm1
 
 shp_adm1 <- malariaAtlas::getShp(ISO = adm0_shp_available$iso, admin_level = "admin1")
-```
 
-Convert to `sf` objects
-```{r}
+
+# Convert to sf object ----------------------------------------------------
+
 sf_adm0 <- shp_adm0 %>% sf::st_as_sf()
 sf_adm1 <- shp_adm1 %>% sf::st_as_sf()
-```
 
-Subset sSA
-```{r}
+# Subset SSA --------------------------------------------------------------
+
 africa_iso3 <- c(
-  "DZA", "AGO", "BEN", "BWA", "BFA", "BDI", "CPV", "CMR", "CAF", "TCD", 
-  "COM", "COD", "COG", "DJI", "EGY", "GNQ", "ERI", "SWZ", "ETH", "GAB", 
-  "GMB", "GHA", "GIN", "GNB", "CIV", "KEN", "LSO", "LBR", "LBY", "MDG", 
-  "MWI", "MLI", "MRT", "MUS", "MAR", "MOZ", "NAM", "NER", "NGA", "RWA", 
-  "STP", "SEN", "SYC", "SLE", "SOM", "ZAF", "SSD", "SDN", "TGO", "TUN", 
+  "DZA", "AGO", "BEN", "BWA", "BFA", "BDI", "CPV", "CMR", "CAF", "TCD",
+  "COM", "COD", "COG", "DJI", "EGY", "GNQ", "ERI", "SWZ", "ETH", "GAB",
+  "GMB", "GHA", "GIN", "GNB", "CIV", "KEN", "LSO", "LBR", "LBY", "MDG",
+  "MWI", "MLI", "MRT", "MUS", "MAR", "MOZ", "NAM", "NER", "NGA", "RWA",
+  "STP", "SEN", "SYC", "SLE", "SOM", "ZAF", "SSD", "SDN", "TGO", "TUN",
   "UGA", "TZA", "ZMB", "ZWE", "ESH"
 )
 
 sf_adm0_africa <- sf_adm0 %>% filter(iso %in% africa_iso3)
 sf_adm1_africa <- sf_adm1 %>% filter(iso %in% africa_iso3)
-```
 
-Save base maps
-```{r}
-saveRDS(sf_adm0, "data_derived/sf_admin0.rds")
-saveRDS(sf_adm1, "data_derived/sf_admin1.rds")
-saveRDS(sf_adm0_africa, "data_derived/sf_admin0_africa.rds")
-saveRDS(sf_adm1_africa, "data_derived/sf_admin1_africa.rds")
-```
+# Save base maps ----------------------------------------------------------
+saveRDS(sf_adm0, "analysis/data_derived/sf_admin0.rds")
+saveRDS(sf_adm1, "analysis/data_derived/sf_admin1.rds")
+saveRDS(sf_adm0_africa, "analysis/data_derived/sf_admin0_africa.rds")
+saveRDS(sf_adm1_africa, "analysis/data_derived/sf_admin1_africa.rds")

analysis/02_covariate_maps.R

@@ -0,0 +1,22 @@
+# Load pkgs ---------------------------------------------------------------
+
+library(sf)
+library(raster)
+library(tidyverse)
+
+# Read data ---------------------------------------------------------------
+
+covars <- readRDS("analysis/data_raw/final_covariates.rds")
+sf_adm0 <- readRDS("analysis/data_derived/sf_admin0_africa.rds") # adm0 map data
+
+
+# Unlist and combine dfs in long format with map --------------------------
+
+# Use the custom (but very specific) get_map_covar_data() fxn
+map_data <- get_map_covar_data(covars)
+
+# Plot all covariates over time facets ------------------------------------
+
+covars_to_plot <- map_data %>% st_drop_geometry() %>% select(-c(iso3c:year)) %>% names()
+
+plot_spatiotemporal_map(map_data, sf_adm0, covars_to_plot)