Merge pull request #16 from IDEELResearch/feature/map_raster

updated create_raster_object.R function, added maps of grff fake data…

DESCRIPTION

@@ -14,8 +14,8 @@ Depends:
     R (>= 2.10)
 Imports:
     ggplot2,
-    here,
     magrittr,
+    here,
     ragg,
     stringr,
     terra,

NAMESPACE

@@ -7,10 +7,12 @@ importFrom(dplyr,mutate)
 importFrom(dplyr,relocate)
 importFrom(dplyr,rename)
 importFrom(ggplot2,aes)
+importFrom(ggplot2,coord_fixed)
 importFrom(ggplot2,coord_sf)
 importFrom(ggplot2,element_rect)
 importFrom(ggplot2,facet_wrap)
 importFrom(ggplot2,geom_sf)
+importFrom(ggplot2,geom_tile)
 importFrom(ggplot2,ggplot)
 importFrom(ggplot2,ggsave)
 importFrom(ggplot2,scale_fill_viridis_c)

R/create_raster_object.R

@@ -1,54 +1,65 @@
-#' Function to convert a shapefile object to a raster object
+#' Create a Raster Object from a Spatial Object
 #'
-#' This function converts a sf shape file into a raster object.
+#' This function converts a spatial object into a raster object by rasterizing it based on a specified attribute field.
 #'
-#' @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'.
+#' @param shape_obj An `sf` or `SpatVector` object representing the spatial data to be rasterized.
+#' @param raster_field A character string specifying the name of the attribute field in `shape_obj` to use for raster cell values.
+#' @param res Numeric value specifying the resolution of the raster grid (default is `0.05`).
+#' @return A `SpatRaster` object created by rasterizing the input spatial data. Each cell in the raster corresponds to an area in the spatial object, with cell values based on the specified attribute field.
 #' @importFrom terra rast rasterize ext
 #' @examples
+#' \dontrun{
+#' # Load required libraries
+#' library(terra)
+#' library(sf)
+#'
 #' # Example of usage:
-#' # Read in Africa shape file
-#' # rds_file = "analysis/data_derived/sf_admin0_africa.rds"
-#' # africa_shp <- readRDS(file = rds_file)
+#' # Read in Africa shapefile
+#' africa_shp <- st_read("path_to_africa_shapefile.shp")
 #'
 #' # 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")
+#' uganda_shp <- africa_shp[africa_shp$iso == "UGA", ]
+#' uganda_raster <- create_raster_object(shape_obj = uganda_shp, raster_field = "iso")
+#'
+#' # Plot Uganda raster with Ugandan border
+#' 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"
-#'   # )
+#' writeRaster(
+#'   uganda_raster,
+#'   filename = "uganda_raster_admin0.tif",
+#'   filetype = "GTiff",
+#'   overwrite = TRUE
+#' )
+#'
+#' # Plot raster via ggplot2
+#' library(ggplot2)
+#' 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, fill = iso)) +
+#'   coord_equal() +
+#'   theme_minimal() +
+#'   labs(
+#'     title = "Uganda Raster",
+#'     x = "Longitude",
+#'     y = "Latitude"
+#'   )
 #'
 #' # Create African raster
-#' # africa_raster <- create_raster_object(shape_obj=africa_shp)
-#' # plot(africa_raster)
+#' africa_raster <- create_raster_object(shape_obj = africa_shp, raster_field = "iso")
+#' plot(africa_raster)
+#' }
 
-create_raster_object <- function(shape_obj, res = 0.05) {
+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)
 
   # Rasterize the shape object using the 'iso' column (country codes) as values for the raster cells
-  rasterized <- rasterize(shape_obj, raster_grid, field = "iso", background = NA)
+  rasterized <- rasterize(shape_obj, raster_grid, field = raster_field, background = NA)
+
+  #change name of raster of the raster layer
+  names(rasterized) <- raster_field
 
   # Return the rasterized object
   return(rasterized)

R/plot_spatiotemporal_map_raster.R

@@ -0,0 +1,42 @@
+#' 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`.
+#'
+#' @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
+#' @importFrom scales percent_format
+#' @importFrom dplyr filter
+#' @importFrom rlang sym .data
+#' @importFrom stringr str_glue
+#' @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)
+#' # variable_names <- c("prev_Q50", "prev_Q2.5", "prev_Q97.5")
+#' # plot_spatiotemporal_map_raster(stacked_rasters, variable_names)
+
+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_raster_grff_{v}.png")
+    ggsave(filename, p, width = 8, height = 6, dpi = 300)
+  }
+}

analysis/02_grff_posterior_plotting.Rmd

@@ -0,0 +1,52 @@
+---
+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))
+
+```

analysis/XX_map_of_resistance_data.R

@@ -0,0 +1,89 @@
+library(sf)       # for shapefile handling
+library(ggplot2)  # for plotting
+library(dplyr)    # for data manipulation
+
+# Expanded dataframe with 15 African countries and varying admin 1 regions
+prevalence_data <- data.frame(
+  country = c("Uganda", "Uganda", "Uganda",
+              "Kenya", "Kenya", "Kenya",
+              "Tanzania", "Tanzania", "Tanzania",
+              "Nigeria", "Nigeria", "Nigeria", "Nigeria",
+              "Ghana", "Ghana",
+              "South Africa", "South Africa", "South Africa",
+              "Zambia", "Zambia",
+              "Zimbabwe", "Zimbabwe", "Zimbabwe",
+              "Malawi",
+              "Ethiopia", "Ethiopia", "Ethiopia", "Ethiopia",
+              "Senegal", "Senegal",
+              "Côte d'Ivoire", "Côte d'Ivoire", "Côte d'Ivoire", "Côte d'Ivoire",
+              "Angola", "Angola", "Angola", "Angola", "Angola",
+              "Rwanda", "Rwanda", "Rwanda",
+              "Burkina Faso", "Burkina Faso", "Burkina Faso", "Burkina Faso",
+              "Togo", "Togo",
+              "Cameroon", "Cameroon", "Cameroon", "Cameroon", "Cameroon",
+              "Mali", "Mali", "Mali",
+              "Gambia",
+              "Somalia", "Somalia", "Somalia",
+              "Namibia", "Namibia", "Namibia", "Namibia"),
+  name_1 = c("Central", "Eastern", "Northern",  # Uganda
+             "Nairobi", "Mombasa", "Kisumu",  # Kenya
+             "Dodoma", "Arusha", "Mwanza",  # Tanzania
+             "Lagos", "Kano", "Abuja", "Port Harcourt",  # Nigeria
+             "Greater Accra", "Ashanti",  # Ghana
+             "Gauteng", "Western Cape", "KwaZulu-Natal",  # South Africa
+             "Lusaka", "Copperbelt",  # Zambia
+             "Harare", "Bulawayo", "Manicaland",  # Zimbabwe
+             "Blantyre",  # Malawi
+             "Addis Ababa", "Oromia", "Amhara", "Tigray",  # Ethiopia
+             "Dakar", "Thies",  # Senegal
+             "Abidjan", "Bouaké", "Yamoussoukro", "Korhogo",  # Côte d'Ivoire
+             "Luanda", "Huambo", "Benguela", "Malanje", "Kwanza Sul",  # Angola
+             "Kigali", "Western", "Northern",  # Rwanda
+             "Centre", "Boucle du Mouhoun", "Hauts-Bassins", "Sud-Ouest",  # Burkina Faso
+             "Maritime", "Plateaux",  # Togo
+             "Littoral", "Ouest", "Nord-Ouest", "Sud-Ouest", "Adamaoua",  # Cameroon
+             "Bamako", "Kayes", "Sikasso",  # Mali
+             "Banjul",  # Gambia
+             "Mogadishu", "South West", "Galmudug",  # Somalia
+             "Khomas", "Erongo", "Otjozondjupa", "Hardap"),  # Namibia,  # Namibia
+  year = sample(2018:2022, 64, replace = TRUE),  # Random years between 2018 and 2022
+  sample_size = sample(0:180, 64, replace = TRUE),
+  prevalence = sample(seq(0,1,0.001), 64, replace = TRUE)
+)
+
+# Read in Africa shape file
+rds_file_admin0 = "analysis/data_derived/sf_admin0_africa.rds"
+rds_file_admin1 = "analysis/data_derived/sf_admin1_africa.rds"
+africa_shp_admin0 <- readRDS(file = rds_file_admin0)
+africa_shp_admin1 <- readRDS(file = rds_file_admin1)
+
+# Calculate the centroids for each MULTIPOLYGON
+sf_use_s2(FALSE)
+africa_admin1_longlat <- africa_shp_admin1 %>%
+  st_geometry() %>% #obtain geometry
+  st_centroid() %>% #obtain centroid of geometry
+  st_coordinates() %>% #obtain coordinates of the centroids (long, lat)
+  cbind(africa_shp_admin1) %>%
+  rename("lon" = "X",
+         "lat" = "Y")
+
+# Add prevalence data to admin 1 shape file
+africa_admin1_longlat_prev <- africa_admin1_longlat %>%
+  left_join(prevalence_data, by = "name_1") %>%
+  filter(!is.na(sample_size))
+
+# Plot the Africa map coloring districts by sample_size
+africa_map_sample_fill <- ggplot() +
+  facet_wrap(~year) +
+  geom_sf(data = africa_admin1_longlat_prev, aes(fill = sample_size), color = "darkgrey", lwd = 0.05) +
+  geom_sf(data = africa_shp_admin0, fill = NA, color = "black", show.legend = FALSE, lwd = 0.05) +
+  theme_void(base_size = 14) +
+  labs(fill = "Sample Size (N)") +
+  scale_fill_viridis_c() +
+  theme(legend.position = "bottom",
+        plot.background = element_rect(fill = "white", color="white"))
+
+ggsave(filename="analysis/plots/africa_map_sample_fill.png", africa_map_sample_fill)
+
+
+