Merge pull request #18 from IDEELResearch/feature/map_slice_year_raster

Feature/map slice year raster

DESCRIPTION

@@ -21,9 +21,9 @@ Imports:
     terra,
     dplyr,
     sf,
-    purrr,
     rlang,
-    scales
+    scales,
+    purrr
 Suggests:
     knitr,
     rmarkdown,

NAMESPACE

@@ -20,6 +20,7 @@ importFrom(ggplot2,theme)
 importFrom(ggplot2,theme_void)
 importFrom(grDevices,dev.off)
 importFrom(grDevices,pdf)
+importFrom(here,here)
 importFrom(magrittr,"%>%")
 importFrom(purrr,imap_dfr)
 importFrom(rlang,.data)

R/create_raster_object.R

@@ -53,7 +53,7 @@
 
 create_raster_object <- function(shape_obj, raster_field, res = 0.05) {
   # Create an empty raster grid covering the extent of the shape object with the specified resolution
-  raster_grid <- rast(ext = ext(shape_obj), resolution = res)
+  raster_grid <- rast(ext = ext(shape_obj) + 0.0001, resolution = res)
 
   # Rasterize the shape object using the 'iso' column (country codes) as values for the raster cells
   rasterized <- rasterize(shape_obj, raster_grid, field = raster_field, background = NA)

R/plot_spatiotemporal_map_raster.R

@@ -1,42 +1,59 @@
-#' Plot spatiotemporal maps from raster data for specified variables
+#' Plot Spatiotemporal Maps from Raster Data for Specified Variables
 #'
-#' This function plots spatiotemporal maps from a raster object for each variable specified in the vector `variable_names`.
-#' It converts the raster data to a data frame, filters out missing values, and creates a plot using **ggplot2** for each variable.
-#' The plots are saved as PNG files in the 'analysis/plots' directory with filenames corresponding to each variable.
-#' Note that the raster object should contain layers corresponding to the variables specified in `variable_names`.
+#' This function plots spatiotemporal maps from a raster dataframe for each variable
+#' specified in the vector `variable_names`. It converts the raster data to a
+#' data frame, filters out missing values, and creates a plot using **ggplot2**
+#' for each variable. If `plot_by_year` is `TRUE`, the function will generate
+#' faceted plots by year, requiring a valid `year_array`. Otherwise, the plots
+#' are generated without faceting. The resulting plots are saved as PNG files
+#' in the 'analysis/plots' directory with filenames corresponding to each variable.
 #'
-#' @param raster A raster object (e.g., a SpatRaster from the **terra** package) containing the data to be plotted.
-#' @param variable_names A vector of variable names corresponding to the layers in the raster object to plot. The function will loop over these names.
-#' @return The function does not return a value; it saves plots as PNG files in the 'analysis/plots' directory.
-#' @importFrom ggplot2 ggplot aes geom_tile scale_fill_viridis_c coord_fixed theme_void theme element_rect ggsave
+#' @param raster_df A data frame containing the raster data to be plotted. The data frame should contain longitude (`long`), latitude (`lat`), time, and variables for plotting.
+#' @param variable_names A vector of variable names corresponding to columns in the `raster_df` to plot. The function loops over these names to generate plots.
+#' @param plot_by_year A logical flag indicating whether to create faceted plots by year (default is `TRUE`). If `TRUE`, the `raster_df` must contain a `time` variable for faceting.
+#' @return The function does not return a value but saves the plots as PNG files in the 'analysis/plots' directory.
+#' @importFrom ggplot2 ggplot aes geom_tile scale_fill_viridis_c coord_fixed theme_void theme element_rect ggsave facet_wrap
 #' @importFrom scales percent_format
 #' @importFrom dplyr filter
 #' @importFrom rlang sym .data
 #' @importFrom stringr str_glue
+#' @importFrom here here
 #' @examples
 #' # Example of how to use the function
-#' # library(terra)
-#' # Load your raster data (assuming 'stacked_rasters' is a SpatRaster with
-#' # layers 'prev_Q50', 'prev_Q2.5', 'prev_Q97.5')
-#' # grff_summary <- readRDS("data_derived/GRFF_summary_example.rds")
-#' # sf_grff_summary <- st_as_sf(grff_summary, coords = c("lon", "lat"), crs=4326)
+#' # Assuming 'raster_df' is a data frame with columns 'long', 'lat', 'time', and variables 'prev_Q50', 'prev_Q2.5', etc.
 #' # variable_names <- c("prev_Q50", "prev_Q2.5", "prev_Q97.5")
-#' # plot_spatiotemporal_map_raster(stacked_rasters, variable_names)
+#' # plot_spatiotemporal_map_raster(raster_df, variable_names, plot_by_year = TRUE)
+#'
+plot_spatiotemporal_map_raster <- function(raster_df, variable_names, plot_by_year = TRUE) {
+
+  if (plot_by_year) {
+    for (v in variable_names) {
+      p <- raster_df %>%
+        filter(!is.na(.data[[v]])) %>%
+        ggplot(aes(x = .data$long, y = .data$lat)) +
+        geom_tile(aes(fill = !!sym(v)), color = NA) +
+        scale_fill_viridis_c(option = "viridis", labels = scales::percent_format(accuracy = 1)) +
+        coord_fixed() +
+        facet_wrap(~time) +
+        theme_void(base_size = 14) +
+        theme(plot.background = element_rect(fill = "white", color = "white"))
 
-plot_spatiotemporal_map_raster <- function(raster, variable_names){
-  raster_df <- as.data.frame(raster, xy = TRUE, na.rm = TRUE)
-  for (v in variable_names) {
-    p <- raster_df %>%
-      filter(!is.na(.data[[v]])) %>%
-      ggplot(aes(x = .data$x, y = .data$y)) +
-      geom_tile(aes(fill = !!sym(v)), color = NA) +
-      scale_fill_viridis_c(option = "viridis", labels = scales::percent_format(accuracy = 1)) +
-      coord_fixed() +
-      #facet_wrap(~year) +
-      theme_void(base_size = 14) +
-      theme(plot.background = element_rect(fill = "white", color = "white"))
+      filename <- str_glue("analysis/plots/map_slice_raster_grff_{v}.png")
+      ggsave(filename = here(filename), plot = p, width = 8, height = 6, dpi = 300)
+    }
+  } else {
+    for (v in variable_names) {
+      p <- raster_df %>%
+        filter(!is.na(.data[[v]])) %>%
+        ggplot(aes(x = .data$long, y = .data$lat)) +
+        geom_tile(aes(fill = !!sym(v)), color = NA) +
+        scale_fill_viridis_c(option = "viridis", labels = scales::percent_format(accuracy = 1)) +
+        coord_fixed() +
+        theme_void(base_size = 14) +
+        theme(plot.background = element_rect(fill = "white", color = "white"))
 
-    filename <- str_glue("analysis/plots/map_raster_grff_{v}.png")
-    ggsave(filename, p, width = 8, height = 6, dpi = 300)
+      filename <- str_glue("analysis/plots/map_raster_grff_{v}.png")
+      ggsave(filename = here(filename), plot = p, width = 8, height = 6, dpi = 300)
+    }
   }
 }

analysis/04_plotting_grff_sliced_rasterg.Rmd

@@ -0,0 +1,84 @@
+---
+title: "Plotting Posterior Median of GRFF"
+author: "Cecile Meier-Scherling"
+date: "`r Sys.Date()`"
+output: 
+  html_document: 
+    code_folding: hide 
+    fig_width: 8
+    fig_height: 6
+    theme: cosmo
+---
+
+```{r setup, include=FALSE}
+# Set global chunk options to show code and hide warnings and messages
+knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
+
+# Load necessary libraries
+# devtools::install_github("malaria-atlas-project/malariaAtlas")
+library(ggplot2)  # For creating plots
+library(raster)   # For working with raster data
+library(terra)    # For spatial data manipulation
+library(dplyr)    # For data manipulation
+library(sf)       # For handling spatial features (simple features)
+library(viridis)  # For color scales
+library(purrr)    # For functional programming tools like map and imap
+```
+
+Map of GRFF output sliced by year
+```{r}
+# Load the GRFF summary data from an RDS file
+grff_summary_year <- readRDS("data_derived/GRFF_example1.rds")
+
+# Load shape file for Uganda
+sf_uganda <- readRDS("data_derived/sf_admin0_africa.rds") %>%
+  filter(iso == "UGA")
+
+# Combine post_summary data frames for all time points
+# and add a 'time' column using imap_dfr to retain the time information
+post_summary_combined <- imap_dfr(grff_summary_year, ~ .x$post_summary %>%
+  mutate(time = .x$time))
+
+# Convert the combined dataframe to a spatial (sf) object, specifying longitude and latitude
+sf_grff_summary_year <- st_as_sf(post_summary_combined, coords = c("lon", "lat"), crs = 4326)
+
+# Create an empty dataframe to store rasterized data for each time point
+raster_time_slice_df <- data.frame()
+
+# Loop through each unique time point in the data
+for (time_point in unique(sf_grff_summary_year$time)) {
+  
+  # Filter the sf object to get data for the current time point
+  sf_time_slice <- sf_grff_summary_year %>%
+    filter(time == time_point)
+  
+  # Rasterize the sf object for different variables (prev_Q2.5, prev_Q50, etc.)
+  raster_Q2.5 <- create_raster_object(sf_time_slice, "prev_Q2.5")  # Raster for 2.5% quantile
+  raster_Q50 <- create_raster_object(sf_time_slice, "prev_Q50")    # Raster for median (50% quantile)
+  raster_Q97.5 <- create_raster_object(sf_time_slice, "prev_Q97.5")  # Raster for 97.5% quantile
+  raster_sd <- create_raster_object(sf_time_slice, "prev_sd")      # Raster for standard deviation
+  
+  # Combine (stack) the rasters for this time point
+  raster_time_slice <- c(raster_Q2.5, raster_Q50, raster_Q97.5, raster_sd)
+  
+  # Mask the raster using the Uganda shapefile to limit the data to the region of interest
+  masked_raster_grff_summary_year <- mask(raster_time_slice, vect(sf_uganda))
+  
+  # Convert the masked raster to a dataframe and add the 'time' column
+  raster_time_slice_df_current <- as.data.frame(masked_raster_grff_summary_year, xy = TRUE) %>%
+    rename("long" = x, "lat" = y) %>%
+    mutate(time = time_point)
+  
+  # Append the current time slice dataframe to the main dataframe
+  raster_time_slice_df <- rbind(raster_time_slice_df, raster_time_slice_df_current)
+}
+
+# Define the variable names to plot from the rasterized data
+variable_names <- c("prev_Q2.5", "prev_Q50", "prev_Q97.5", "prev_sd")
+
+# Filter the raster data to include time slices every 6 months (roughly 365-day intervals)
+raster_time_slice_df_6m <- raster_time_slice_df %>% filter(time %in% seq(1, 6991, 365))
+
+# Plot the spatiotemporal maps for the filtered raster data, faceting by year
+plot_spatiotemporal_map_raster(raster_time_slice_df_6m, variable_names, plot_by_year = TRUE)
+```