code to make multipanel figure of predicted species presence

R/map_predicted_distributions.R

@@ -0,0 +1,89 @@
+#' Make a multi-panel grid of predicted species presence
+#'
+#'
+#' @title map_predicted_distributions
+#' @return 
+#' @author Cecilia Sanchez
+
+map_predicted_distributions <- function(){
+  # import a world countries map
+  countries <- geodata::world(resolution = 3, path = tempdir())
+
+  # countries of interest
+  WA_names <- c("Afghanistan", "Armenia", "Azerbaijan", "Bahrain", "Georgia",
+                "Iran", "Iraq", "Israel", "Jordan", "Kuwait", "Lebanon", "Oman",
+                "Pakistan", "Palestine", "Qatar", "Saudi Arabia", "Syria", 
+                "Turkey", "United Arab Emirates", "Yemen")
+
+  # get ISO3 country codes for WA countries
+  WA_codes <- geodata::country_codes() %>% 
+    dplyr::filter(NAME %in% WA_names) %>% 
+    pull(ISO3)
+
+  # subset the full countries map to just WA countries
+  WA_countries <- countries[countries$GID_0 %in% WA_codes]
+  # make an sf object for later plotting
+  WA_countries_sf <- st_as_sf(WA_countries)
+
+  # get list of all generated predictions rasters
+  bat_ras_names <- list.files(path = "outputs/", pattern = "full", 
+                              full.names = TRUE)
+
+  maps_list <- list()
+
+  for(i in 1:length(bat_ras_names)){
+
+    # pull one raster
+    bat_ras <- terra::rast(bat_ras_names[i])
+    # get the species name for adding plot title
+    species_name <- strsplit(bat_ras_names[i], '[_.]')[[1]][2]
+
+    # crop/mask to WA countries
+    bat_ras_cropped <- terra::crop(bat_ras, WA_countries, mask = T)
+    # align CRS
+    crs(bat_ras_cropped) <- crs(WA_countries)
+
+    # plot using tidyterra functionality
+    maps_list[[i]] <- ggplot() +
+      geom_spatvector(data = WA_countries, fill = "gray95") +
+      geom_spatraster(data = bat_ras_cropped) +
+      scale_fill_viridis_c(option = "D", na.value = "transparent", 
+                           name = "Predicted probability \nof species presence",
+                           breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) +
+      geom_spatvector(data = WA_countries, fill = "transparent", lwd = 0.5) +
+      labs(x = "Longitude", y = "Latitude") +
+      ggtitle(species_name) +
+      theme_bw() +
+      theme(axis.text = element_text(color = "black"))
+  }
+
+  # create initial grid with no legends
+  plot_grid(maps_list[[1]] + theme(legend.position = "none"), 
+            maps_list[[2]] + theme(legend.position = "none"), 
+            maps_list[[3]] + theme(legend.position = "none"), 
+            maps_list[[4]] + theme(legend.position = "none"), 
+            maps_list[[5]] + theme(legend.position = "none"), 
+            maps_list[[6]] + theme(legend.position = "none"), 
+            maps_list[[7]] + theme(legend.position = "none"), 
+            maps_list[[8]] + theme(legend.position = "none"), 
+            maps_list[[9]] + theme(legend.position = "none"), 
+            maps_list[[10]] + theme(legend.position = "none"), 
+            maps_list[[11]] + theme(legend.position = "none"), 
+            maps_list[[12]] + theme(legend.position = "none"), 
+            ncol = 3) -> my_grid
+
+  # make one legend to go at the bottom
+  legend_below <- get_legend(maps_list[[1]] + 
+                               guides(color = guide_legend(nrow = 1)) +
+                               theme(legend.position = "bottom"))
+
+  # add initial grid plus the single legend
+  plot_grid(my_grid, legend_below, ncol = 1, rel_heights = c(1, 0.1))
+
+
+  rm(list = "bat_ras", "bat_ras_cropped", "countries", "WA_countries", 
+     "WA_countries_sf", "my_grid", "maps_list")
+
+  ggsave("outputs/species_predictions.png", width = 8.5, height = 11, 
+         units = "in", dpi = 600)
+}

_targets.R

@@ -119,9 +119,8 @@ analysis_targets <- tar_plan(
 
 ## Outputs
 outputs_targets <- tar_plan(
-  ## This is a placeholder for any targets that produces outputs such as
-  ## tables of model outputs, plots, etc. Delete or keep empty if you will not
-  ## produce any of these types of outputs
+
+  multipanel = map_predicted_distributions()
 )
 
 
@@ -136,28 +135,6 @@ report_targets <- tar_plan(
   # )
 )
 
-## Deploy targets
-deploy_targets <- tar_plan(
-  ## This is a placeholder for any targets that are meant to deploy reports or
-  ## any outputs externally e.g., website, Google Cloud Storage, Amazon Web
-  ## Services buckets, etc. Delete or keep empty if you will not perform any
-  ## deployments. The aws_s3_upload function requires AWS credentials to be loaded
-  ## but will print a warning and do nothing if not
-
-  # html_files = containerTemplateUtils::get_file_paths(tar_obj = example_report,
-  #                                                     pattern = "\\.html$"),
-  # uploaded_report = containerTemplateUtils::aws_s3_upload(html_files,
-  #                                                       bucket = Sys.getenv("AWS_BUCKET"),
-  #                                                       error = FALSE,
-  #                                                       file_type = "html"),
-  # email_updates= 
-  #   containerTemplateUtils::send_email_update(
-  #     to = strsplit(Sys.getenv("EMAIL_RECIPIENTS"),";")[[1]],
-  #     from = Sys.getenv("EMAIL_SENDER"),
-  #     project_name = "My Project",
-  #     attach = TRUE
-  #   )
-)
 
 # List targets -----------------------------------------------------------------
 
@@ -166,6 +143,5 @@ list(
   data_processing_targets,
   analysis_targets,
   outputs_targets,
-  report_targets,
-  deploy_targets
+  report_targets
 )

packages.R

@@ -39,3 +39,6 @@ library(sf)
 library(ENMeval)
 library(rJava)
 library(geodata)
+library(tidyterra)
+library(viridis)
+library(cowplot)

renv.lock

@@ -363,6 +363,19 @@
         "rlang"
       ]
     },
+    "cowplot": {
+      "Package": "cowplot",
+      "Version": "1.1.1",
+      "Source": "Repository",
+      "Repository": "CRAN",
+      "Hash": "b418e8423699d11c7f2087c2bfd07da2",
+      "Requirements": [
+        "ggplot2",
+        "gtable",
+        "rlang",
+        "scales"
+      ]
+    },
     "cpp11": {
       "Package": "cpp11",
       "Version": "0.4.5",
@@ -710,6 +723,16 @@
       "Hash": "15e9634c0fcd294799e9b2e929ed1b86",
       "Requirements": []
     },
+    "geodata": {
+      "Package": "geodata",
+      "Version": "0.5-8",
+      "Source": "Repository",
+      "Repository": "CRAN",
+      "Hash": "93d238706f3fe7c326a7d0bbc2fecdb8",
+      "Requirements": [
+        "terra"
+      ]
+    },
     "geosphere": {
       "Package": "geosphere",
       "Version": "1.5-18",
@@ -2118,6 +2141,26 @@
         "withr"
       ]
     },
+    "tidyterra": {
+      "Package": "tidyterra",
+      "Version": "0.4.0",
+      "Source": "Repository",
+      "Repository": "CRAN",
+      "Hash": "c50e45d8fd32d10e40ff6d9f03e22249",
+      "Requirements": [
+        "cli",
+        "data.table",
+        "dplyr",
+        "ggplot2",
+        "magrittr",
+        "rlang",
+        "scales",
+        "sf",
+        "terra",
+        "tibble",
+        "tidyr"
+      ]
+    },
     "tidyverse": {
       "Package": "tidyverse",
       "Version": "1.3.2",