diff --git a/.DS_Store b/.DS_Store index 9f53299..aaad3ff 100644 Binary files a/.DS_Store and b/.DS_Store differ diff --git a/DESCRIPTION b/DESCRIPTION index c8028d8..d4a24ec 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -21,9 +21,9 @@ Imports: terra, dplyr, sf, - purrr, rlang, - scales + scales, + purrr Suggests: knitr, rmarkdown, diff --git a/NAMESPACE b/NAMESPACE index 19bb5a5..921e6a7 100644 --- a/NAMESPACE +++ b/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) diff --git a/R/create_raster_object.R b/R/create_raster_object.R index 5d5827e..34525b3 100644 --- a/R/create_raster_object.R +++ b/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) diff --git a/R/plot_spatiotemporal_map_raster.R b/R/plot_spatiotemporal_map_raster.R index be6b229..a444e14 100644 --- a/R/plot_spatiotemporal_map_raster.R +++ b/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) + } } } diff --git a/analysis/.DS_Store b/analysis/.DS_Store index 2eae1d3..caf5a4e 100644 Binary files a/analysis/.DS_Store and b/analysis/.DS_Store differ diff --git a/analysis/02_grff_posterior_plotting.Rmd b/analysis/02_grff_posterior_plotting.Rmd deleted file mode 100644 index 2917a31..0000000 --- a/analysis/02_grff_posterior_plotting.Rmd +++ /dev/null @@ -1,52 +0,0 @@ ---- -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} -knitr::opts_chunk$set(echo = TRUE, warning = F, message = F) - -# devtools::install_github("malaria-atlas-project/malariaAtlas") - -library(ggplot2) -library(raster) -library(terra) -library(dplyr) -library(sf) -library(viridis) -``` - -Loading grff data -```{r load data} -grff_draws <- readRDS("data_derived/GRFF_draws_example.rds") -grff_summary <- readRDS("data_derived/GRFF_summary_example.rds") -``` - -Convert dataframe to raster -```{r convert df to sf object} -# Convert dataframe to sf object -sf_grff_summary <- st_as_sf(grff_summary, coords = c("lon", "lat"), crs=4326) - -# Rasterize the sf object -raster_Q2.5 <- create_raster_object(sf_grff_summary, "prev_Q2.5", res = 0.1) -raster_Q50 <- create_raster_object(sf_grff_summary, "prev_Q50", res = 0.1) -raster_Q97.5 <- create_raster_object(sf_grff_summary, "prev_Q97.5", res = 0.1) -raster_sd <- create_raster_object(sf_grff_summary, "prev_sd", res = 0.1) - -# Stack the Rasters Together -raster_grff_summary <- c(raster_Q2.5, raster_Q50, raster_Q97.5, raster_sd) -``` - -Plot rasters -```{r map Q50} -#plot_spatiotemporal -plot_spatiotemporal_map_raster(raster_grff_summary, names(raster_grff_summary)) - -``` diff --git a/analysis/04_plotting_grff_sliced_rasterg.Rmd b/analysis/04_plotting_grff_sliced_rasterg.Rmd new file mode 100644 index 0000000..5382833 --- /dev/null +++ b/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) +``` diff --git a/analysis/data_derived/GRFF_draws_example.rds b/analysis/data_derived/GRFF_draws_example.rds new file mode 100644 index 0000000..4ff3faf Binary files /dev/null and b/analysis/data_derived/GRFF_draws_example.rds differ diff --git a/analysis/data_derived/GRFF_summary_example.rds b/analysis/data_derived/GRFF_summary_example.rds new file mode 100644 index 0000000..0d6c9e2 Binary files /dev/null and b/analysis/data_derived/GRFF_summary_example.rds differ diff --git a/analysis/data_derived/uganda_raster_admin0.tif.aux.xml b/analysis/data_derived/uganda_raster_admin0.tif.aux.xml new file mode 100644 index 0000000..b3803ff --- /dev/null +++ b/analysis/data_derived/uganda_raster_admin0.tif.aux.xml @@ -0,0 +1,263 @@ + + + iso + + UGA + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/analysis/plots/africa_map_sample_fill.png b/analysis/plots/africa_map_sample_fill.png new file mode 100644 index 0000000..c2c824d Binary files /dev/null and b/analysis/plots/africa_map_sample_fill.png differ diff --git a/analysis/plots/africa_map_sample_size_points.png b/analysis/plots/africa_map_sample_size_points.png new file mode 100644 index 0000000..5fa984e Binary files /dev/null and b/analysis/plots/africa_map_sample_size_points.png differ diff --git a/analysis/plots/map_slice_raster_grff_prev_Q2.5.png b/analysis/plots/map_slice_raster_grff_prev_Q2.5.png new file mode 100644 index 0000000..3e91749 Binary files /dev/null and b/analysis/plots/map_slice_raster_grff_prev_Q2.5.png differ diff --git a/analysis/plots/map_slice_raster_grff_prev_Q50.png b/analysis/plots/map_slice_raster_grff_prev_Q50.png new file mode 100644 index 0000000..22ba18c Binary files /dev/null and b/analysis/plots/map_slice_raster_grff_prev_Q50.png differ diff --git a/analysis/plots/map_slice_raster_grff_prev_Q97.5.png b/analysis/plots/map_slice_raster_grff_prev_Q97.5.png new file mode 100644 index 0000000..0123118 Binary files /dev/null and b/analysis/plots/map_slice_raster_grff_prev_Q97.5.png differ diff --git a/analysis/plots/map_slice_raster_grff_prev_sd.png b/analysis/plots/map_slice_raster_grff_prev_sd.png new file mode 100644 index 0000000..c9cd02e Binary files /dev/null and b/analysis/plots/map_slice_raster_grff_prev_sd.png differ diff --git a/man/plot_spatiotemporal_map_raster.Rd b/man/plot_spatiotemporal_map_raster.Rd new file mode 100644 index 0000000..f049bde --- /dev/null +++ b/man/plot_spatiotemporal_map_raster.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_spatiotemporal_map_raster.R +\name{plot_spatiotemporal_map_raster} +\alias{plot_spatiotemporal_map_raster} +\title{Plot Spatiotemporal Maps from Raster Data for Specified Variables} +\usage{ +plot_spatiotemporal_map_raster(raster_df, variable_names, plot_by_year = TRUE) +} +\arguments{ +\item{raster_df}{A data frame containing the raster data to be plotted. The data frame should contain longitude (\code{long}), latitude (\code{lat}), time, and variables for plotting.} + +\item{variable_names}{A vector of variable names corresponding to columns in the \code{raster_df} to plot. The function loops over these names to generate plots.} + +\item{plot_by_year}{A logical flag indicating whether to create faceted plots by year (default is \code{TRUE}). If \code{TRUE}, the \code{raster_df} must contain a \code{time} variable for faceting.} +} +\value{ +The function does not return a value but saves the plots as PNG files in the 'analysis/plots' directory. +} +\description{ +This function plots spatiotemporal maps from a raster dataframe for each variable +specified in the vector \code{variable_names}. It converts the raster data to a +data frame, filters out missing values, and creates a plot using \strong{ggplot2} +for each variable. If \code{plot_by_year} is \code{TRUE}, the function will generate +faceted plots by year, requiring a valid \code{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. +} +\examples{ +# Example of how to use the function +# 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(raster_df, variable_names, plot_by_year = TRUE) + +}