regenerated xy scatter for viticulture unweighted

Author: owenssam1

README.md

@@ -49,16 +49,18 @@ install.packages("blockCV")
 
 ### Sitemap
 
-This GitHub project is organized into general sections of our modeling pipeline. Each article can be found in [vignettes](https://github.com/ieco-lab/scari/tree/master/vignettes). Please see this sitemap for navigating both the GitHub repo and our site:
+This GitHub project is organized into two groups: vignettes which generate reports for SLF risk to viticulture based on our analysis, and our full modeling pipeline used to create these reports.
+
+Reports can be generated using vignettes 150-152, which contain example usage of our function [create_risk_report()](vignettes/150_create_risk_report.Rmd) (150) to create reports for global countries and states/provinces (151), and for the USA specifically (152).
+
+Please see this sitemap for a guide to our full modeling pipeline:
 
 * vignette 010: Initialize `scari` and usage of `renv` package for dependencies
-* 020-030: 1. Retrieve and tidy input data for MaxEnt
-* 040-090: 2. SDM modeling pipeline: train global and 3 regional-scale models
-* 100-110: 3. Ensemble Regional-scale SDMs
-* 120-130: 4. Quantify SLF risk
-* 140-142: 5. Quantify Model fit
-* 150-152: Example usage of our function [create_risk_report()](vignettes/150_create_risk_report.Rmd) (150) to create reports for global countries and states/provinces (151), and for the USA specifically (152)
-* 160: Generation of figures for the companion paper
+* 020-030:      1. Retrieve and tidy input data for MaxEnt
+* 040-090:      2. SDM modeling pipeline: train global and 3 regional-scale models
+* 100-110:      3. Ensemble Regional-scale SDMs
+* 120-130, 160: 4. Quantify SLF risk
+* 140-142:      5. Quantify Model fit
 
 ## How to Use this Project
 
@@ -70,6 +72,8 @@ Before diving into this project and our modeling workflow, an end user should:
 4. Run the first vignette, [010_initialize_renv](vignettes/010_initialize_pkg.Rmd), which initializes `renv` and lists our package's dependencies.
 5. See "Get Started" for help in using our package to produce localized reports on SLF risk to viticulture or to recreate our analysis for another invasive species of interest
 
+Once these steps are completed, the end user can get started either generating SLF reports, or following and editing the full modeling pipeline. 
+
 ### Notes about using this package's code
 
 I use some of the following conventions to ensure that the package's .html files render correctly, the code is not overly cumbersome to run, and that data aren't re-downloaded unnecessarily:

vignettes/110_ensemble_regional_models.Rmd

@@ -2712,7 +2712,7 @@ The weighted mean ensemble is more rigorous because it considers climatic extrap
 ```{r unweighted mean pred rasters}
 
 # historical
-pred_hist_mean <- terra::app(
+pred_hist_mean_unweighted <- terra::app(
   x = pred_hist_stack,
   fun = "mean",
   filename = file.path(mypath, "models", "slf_regional_ensemble_v2", "ensemble_regional_unweighted_mean_globe_1981-2010.asc"),
@@ -2722,15 +2722,15 @@ pred_hist_mean <- terra::app(
 
 # CMIP6
 # ssp126
-pred_126_mean <- terra::app(
+pred_126_mean_unweighted <- terra::app(
   x = pred_126_stack,
   fun = "mean",
   filename = file.path(mypath, "models", "slf_regional_ensemble_v2", "ensemble_regional_unweighted_mean_globe_2041-2070_GFDL_ssp126.asc"),
   filetype = "AAIGrid",
   overwrite = FALSE
 )
 # ssp126
-pred_370_mean <- terra::app(
+pred_370_mean_unweighted <- terra::app(
   x = pred_370_stack,
   fun = "mean",
   filename = file.path(mypath, "models", "slf_regional_ensemble_v2", "ensemble_regional_unweighted_mean_globe_2041-2070_GFDL_ssp370.asc"),
@@ -2739,7 +2739,7 @@ pred_370_mean <- terra::app(
 )
 
 # ssp585
-pred_585_mean <- terra::app(
+pred_585_mean_unweighted <- terra::app(
   x = pred_585_stack,
   fun = "mean",
   filename = file.path(mypath, "models", "slf_regional_ensemble_v2", "ensemble_regional_unweighted_mean_globe_2041-2070_GFDL_ssp585.asc"),
@@ -2748,6 +2748,24 @@ pred_585_mean <- terra::app(
 )
 
 
+```
+
+Now I will take an unweighted mean of the above three rasters.
+
+```{r take unweighted mean of weighted pred rasters- CMIP6}
+
+# stack rasters
+pred_CMIP6_mean_stack_unweighted <- c(pred_126_mean_unweighted, pred_370_mean_unweighted, pred_585_mean_unweighted)
+
+# apply mean between layers
+pred_CMIP6_mean_final_unweighted <- terra::app(
+  x = pred_CMIP6_mean_stack_unweighted,
+  fun = "mean",
+  filename = file.path(mypath, "models", "slf_regional_ensemble_v2", "ensemble_regional_unweighted_mean_globe_2041-2070_GFDL_ssp_averaged.asc"),
+  filetype = "AAIGrid",
+  overwrite = FALSE
+)
+
 ```
 
 I will also take the mean of the MTSS and MTP thresholds

vignettes/130_create_suitability_xy_plots_viticultural.Rmd

@@ -3,7 +3,7 @@ title: "Plot shift in SLF risk to important viticultural regions under climate c
 output: rmarkdown::html_document
 author: 
     - "Samuel M. Owens^[Temple University, sam.owens@temple.edu, [ORCID](https://orcid.org/0009-0001-2338-7928)]"
-date: "2025-08-04"
+date: "2026-01-27"
 ---
 
 # Overview
@@ -14,7 +14,7 @@ In this vignette, I will focus on making more specific predictions for SLF estab
 
 Gallien et al, 2012 originally defined this method for discerning the "stage of invasion" for invasive populations according to different scales of SDM. I will apply this method and re-interpret it as a method for assessing the risk of establishment. These two models will not completely agree on the suitability for a particular location, so agreement between the modeled scales can give us more confidence in that prediction. Where the modeled scales disagree, we might make different interpretations of the biological mechanism at play. 
 
-```{r example quadrant plot, echo = FALSE}
+```{r example quadrant plot, echo = FALSE, message = FALSE}
 
 library(tidyverse)
 
@@ -1381,3 +1381,174 @@ We now have risk quadrant plots and tables that we can use to assess the level o
 
 2. Smith, T. 2021, August 11. Evaluating Invasion Stage with SDMs - plantarum.ca. <https://plantarum.ca/2021/08/11/invasion-stage/>.
 
+<!--
+
+# Appendix
+
+## plot unweighted mean versions of suitability plots
+
+This workflow follows steps 1-3, except that it plots based on the unweighted mean of the calculated regional ensemble (rather than an ensemble created based on AUC and ExDet weights).
+
+```{r load in summary files}
+# summary file to extract thresholds from
+
+summary_regional_ensemble_unweighted <- read_csv(file = file.path(mypath, "slf_regional_ensemble_v2", "ensemble_threshold_values_unweighted.csv"))
+
+
+```
+
+Finally, I will load in the global and regional ensemble suitability maps. These will be used both for extracting the new xy suitability values and for plotting.
+
+```{r load in suitability rasters}
+
+# regional_ensemble
+regional_ensemble_1995_unweighted <- terra::rast(
+  x = file.path(mypath, "slf_regional_ensemble_v2", "ensemble_regional_unweighted_mean_globe_1981-2010.asc")
+  )
+
+regional_ensemble_2055_unweighted <- terra::rast(
+  x = file.path(mypath, "slf_regional_ensemble_v2", "ensemble_regional_unweighted_mean_globe_2041-2070_GFDL_ssp_averaged.asc")
+  )
+
+```
+
+### Retrieve/Import xy suitability
+
+These scatter plots will be based on the suitability for the IVR points in both the global and regional_ensemble models. I have already calculated the xy suitability for the global model based on these points, using the function `slfSpread::predict_xy_suitability()`. This function will not work for the regional_ensemble because it calls for a model object, which we did not use to predict the ensemble suitability. So, I will use `terra::extract()` to perform this action.
+
+I will retrieve the suitability values for the regional_ensemble. Instead of returning the coordinates from the map, I will join the coordinates from the original IVR_locations dataset so that the coordinates are exact for joining with other datasets.
+
+```{r retrieve xy suitability for regional_ensemble}
+
+# 1995
+xy_regional_ensemble_1995_unweighted <- terra::extract(
+  x = regional_ensemble_1995_unweighted,
+  y = dplyr::select(IVR_locations, x, y), # points
+  method = "simple",
+  xy = FALSE, # dont return coordinates
+  ID = TRUE,
+  raw = FALSE # return as df
+)
+
+
+# 2055
+xy_regional_ensemble_2055_unweighted <- terra::extract(
+  x = regional_ensemble_2055_unweighted,
+  y = dplyr::select(IVR_locations, x, y), # points
+  method = "simple",
+  xy = FALSE, # dont return coordinates 
+  ID = TRUE,
+  raw = FALSE # return as df
+) 
+
+```
+
+```{r join xy coordinates}
+
+# joining object
+IVR_coordinates <- dplyr::select(IVR_locations, ID, x, y)
+
+# perform join
+xy_regional_ensemble_1995_unweighted <- dplyr::left_join(xy_regional_ensemble_1995_unweighted, IVR_coordinates, by = "ID") %>%
+  dplyr::select(-ID)
+
+xy_regional_ensemble_2055_unweighted <- dplyr::left_join(xy_regional_ensemble_2055_unweighted, IVR_coordinates, by = "ID") %>%
+  dplyr::select(-ID)
+
+```
+
+Now, I will tidy and save the datasets
+
+```{r tidy datasets}
+
+# regional_ensemble datasets
+xy_regional_ensemble_1995_unweighted <- xy_regional_ensemble_1995_unweighted %>%
+  # rename the column for future joining
+  dplyr::rename("xy_regional_ensemble_1995_unweighted" = "mean") %>%
+  dplyr::filter(!is.na(xy_regional_ensemble_1995_unweighted)) %>%
+  dplyr::relocate(x, y) %>%
+  dplyr::mutate(
+      join_col_x = round(x, 5),
+      join_col_y = round(y, 4) # rounding to the 1000s (1km) place to prevent overly sensitive exclusions for UTM data
+    ) %>%
+  dplyr::select(-c(x, y))
+
+xy_regional_ensemble_2055_unweighted <- xy_regional_ensemble_2055_unweighted %>%
+  # rename the column for future joining
+  dplyr::rename("xy_regional_ensemble_2055_unweighted" = "mean") %>%
+  dplyr::filter(!is.na(xy_regional_ensemble_2055_unweighted)) %>%
+  dplyr::relocate(x, y) %>%
+  dplyr::mutate(
+      join_col_x = round(x, 5),
+      join_col_y = round(y, 4) # rounding to the 1000s (1km) place to prevent overly sensitive exclusions for UTM data
+    ) %>%
+  dplyr::select(-c(x, y))
+
+# we lost 1 record because of the NA removal, but that is OK
+
+```
+
+### plot untransformed suitability
+
+```{r join datasets}
+
+# join datasets for plotting
+xy_joined_unweighted <- dplyr::full_join(xy_global_1995, xy_regional_ensemble_1995_unweighted, by = c("join_col_x", "join_col_y")) %>%
+  # join CC datasets
+  dplyr::full_join(., xy_global_2055, by = c("join_col_x", "join_col_y")) %>%
+  dplyr::full_join(., xy_regional_ensemble_2055_unweighted, by = c("join_col_x", "join_col_y")) %>%
+  # order
+  dplyr::relocate(join_col_x, join_col_y, xy_global_1995, xy_global_2055) %>%
+  dplyr::select(-c(ID.x, ID.y))
+
+```
+
+```{r plot SLF suitability values, fig.asp = 1}
+
+# figure annotation title
+# "suitability for Lycorma delicatula establishment in globally important viticultural areas, projected for climate change"
+
+# plot
+xy_joined_unweighted_plot <- ggplot(data = xy_joined_unweighted) +
+   # threshold lines
+   # MTSS thresholds
+   geom_vline(xintercept = as.numeric(summary_global[42, ncol(summary_global)]), linetype = "dashed", linewidth = 0.7) + # global
+   geom_hline(yintercept = as.numeric(summary_regional_ensemble_unweighted[2, 2]), linetype = "dashed", linewidth = 0.7) + # regional_ensemble- there are two MTSS thresholds for this model, but the difference is so small that you will never see it on the plot
+   # historical data
+   geom_point(
+     aes(x = xy_global_1995, y = xy_regional_ensemble_1995_unweighted, shape = "Present"), 
+     size = 2, stroke = 0.7, color = "black", fill = "orchid1"
+     ) +
+   # GFDL ssp370 data
+   geom_point(
+     aes(x = xy_global_2055, y = xy_regional_ensemble_2055_unweighted, shape = "Future | GFDL-ESM4\nmean of ssp126/370/585"), 
+     size = 2, stroke = 0.7, color = "black", fill = "purple3"
+     ) +
+   # axes scaling
+   scale_x_continuous(name =  "'global' model cloglog suitability", limits = c(0, 1), breaks = breaks) + 
+   scale_y_continuous(name =  "'regional_ensemble' model cloglog suitability", limits = c(0, 1), breaks = breaks) +
+   # aesthetics
+   scale_shape_manual(name = "Time period", values = c(21, 21)) +
+   guides(shape = guide_legend(nrow = 1, override.aes = list(size = 2.5), reverse = TRUE)) +
+   theme_bw() +
+   theme(legend.position = "bottom", panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
+   coord_fixed(ratio = 1) 
+
+```
+
+```{r save scatterplot}
+
+ggsave(
+  xy_joined_unweighted_plot, 
+  filename = file.path(
+    here::here(), "vignette-outputs", "figures", "IVR_xy_suitability_IVR_regions_global_regional_ensemble_unweighted.jpg"
+    ),
+  height = 8, 
+  width = 8,
+  device = "jpeg",
+  dpi = "retina"
+  )
+
+```
+
+-->

vignettes/articles/articles_index.Rmd

@@ -45,7 +45,7 @@ Here is the index of articles which outlines the steps in our analysis. Start at
 1. [Create suitability maps per model scale](120_create_suitability_maps.html)
 2. [Create suitability xy plots for viticultural areas](130_create_suitability_xy_plots_viticultural.html)
 3. [Create suitability xy plots for SLF occurrences](131_create_suitability_xy_plots_SLF.html)
-4. [Generate and format figures of results](160_generate_format_figures.html)
+4. [Create viticultural risk tables](160_generate_format_figures.html)
 
 ## Step 6: Validate model fit