6. Project Models to Multiple
Scenarios
Description
Once selected models have been fit and explored, projections to
single or multiple scenarios can be performed. The
project_selected() function is designed for projections to
multiple scenarios (i.e., multiple sets of data, each representing
distinct environmental scenarios). This vignette contains examples of
how to use many of the options available for model projections.
Getting ready
At this point it is assumed that kuenm2 is installed (if
not, see the Main guide). Load
kuenm2 and any other required packages, and define a
working directory (if needed).
Note: functions from other packages (i.e., not from base R or
kuenm2) used in this guide will be displayed as
package::function().
# Load packages
library(kuenm2)
library(terra)
# Current directory
getwd()
# Define new directory
#setwd("YOUR/DIRECTORY") # uncomment and modify if setting a new directory
# Saving original plotting parameters
original_par <- par(no.readonly = TRUE)
Fitted models
To project selected models, a fitted_models object is
required. For detailed information on model fitting, check the vignette
Fit and Explore Selected Models.
The fitted_models object generated in that vignette is
included as an example dataset within the package. Let’s load it.
# Import fitted_model_maxnet
data("fitted_model_maxnet", package = "kuenm2")
# Print fitted models
fitted_model_maxnet
#> fitted_models object summary
#> ============================
#> Species: Myrcia hatschbachii
#> Algortihm: maxnet
#> Number of fitted models: 2
#> Models fitted with 4 replicates
Variables for projections
Predicting models for a single scenario requires a single
SpatRaster object containing the variables (as detailed in
Predict Models to a Single
Scenario). In contrast, projecting models to multiple scenarios
requires a folder that stores variables for each scenario organized in a
certain way.
To ensure the following automated process can correctly track
variables, the data must follow a strict hierarchical directory
structure. At the top level, a Base Directory serves as the
primary container for all project data. Inside this base folder, the
first level of organization consists of subfolders for distinct Time
Periods, such as future years (e.g., “2070”, “2100”) or
paleoclimate eras (e.g., “Mid-holocene”, “LGM”). Within each period
folder, if applicable, you should include subfolders at the
second level for each Emission Scenario (e.g., “ssp126”,
“ssp585”). Finally, within each emission scenario or time period folder,
the third level should include a separate folder for each General
Circulation Model (GCM) (e.g., “BCC-CSM2-MR”, “MIROC6”) to house
the actual raster variables. This structured organization enables the
function to automatically access and process the data according to
period, emission scenario, and GCM.
Variables from WorldClim
The package provides a function to import future climate variables
downloaded from WorldClim (version 2.1). This function renames the files
and organizes them into folders categorized by period/year, emission
scenario (Shared Socioeconomic Pathways; SSPs), and General Circulation
Model (GCM). This simplifies the preparation of climate data, ensuring
all required variables are properly structured for modeling
projections.
To use this function, download the future raster
variables from WorldClim 2.1 and save them all within the same
folder. DO NOT rename the files or the variables, as
the function relies on the patterns provided in the original files to
work properly.
The package also provides an example of raw variables downloaded from
WorldClim 2.1. This example includes bioclimatic predictions for the
periods “2041-2060” and “2081-2100”, for two SSPs (125 and 585) and two
GCMs (ACCESS-CM2 and MIROC6), at 10 arc-minutes resolution.
# See raster files with future variables provided as example
# The data is located in the "inst/extdata" folder.
in_dir <- system.file("extdata", package = "kuenm2")
list.files(in_dir)
#> [1] "CHELSA_LGM_CCSM4.tif"
#> [2] "CHELSA_LGM_CNRM-CM5.tif"
#> [3] "CHELSA_LGM_FGOALS-g2.tif"
#> [4] "CHELSA_LGM_IPSL-CM5A-LR.tif"
#> [5] "CHELSA_LGM_MIROC-ESM.tif"
#> [6] "CHELSA_LGM_MPI-ESM-P.tif"
#> [7] "CHELSA_LGM_MRI-CGCM3.tif"
#> [8] "Current_CHELSA.tif"
#> [9] "Current_variables.tif"
#> [10] "bias_file.tif"
#> [11] "m.gpkg"
#> [12] "wc2.1_10m_bioc_ACCESS-CM2_ssp126_2041-2060.tif"
#> [13] "wc2.1_10m_bioc_ACCESS-CM2_ssp126_2081-2100.tif"
#> [14] "wc2.1_10m_bioc_ACCESS-CM2_ssp585_2041-2060.tif"
#> [15] "wc2.1_10m_bioc_ACCESS-CM2_ssp585_2081-2100.tif"
#> [16] "wc2.1_10m_bioc_MIROC6_ssp126_2041-2060.tif"
#> [17] "wc2.1_10m_bioc_MIROC6_ssp126_2081-2100.tif"
#> [18] "wc2.1_10m_bioc_MIROC6_ssp585_2041-2060.tif"
#> [19] "wc2.1_10m_bioc_MIROC6_ssp585_2081-2100.tif"
#> [20] "world.gpkg"
Note that all variables are in the same folder and retain the
original names provided by WorldClim. You can download these variables
directly from WorldClim
or by using the geodata R package (see example code
below):
# Install geodata if necessary
if (!require("geodata")) {
install.packages("geodata")
}
# Load geodata
library(geodata)
# Create folder to save the raster files
# Here, in a temporary directory
geodata_dir <- file.path(tempdir(), "Future_worldclim")
dir.create(geodata_dir)
# Define GCMs, SSPs and time periods
gcms <- c("ACCESS-CM2", "MIROC6")
ssps <- c("126", "585")
periods <- c("2041-2060", "2061-2080")
# Create a grid of combination of periods, ssps and gcms
g <- expand.grid("period" = periods, "ssps" = ssps, "gcms" = gcms)
g # Each line is a specific scenario for future
# Loop to download variables for each scenario (It can take a while)
lapply(1:nrow(g), function(i) {
cmip6_world(model = g$gcms[i],
ssp = g$ssps[i],
time = g$period[i],
var = "bioc",
res = 10, path = geodata_dir)
})
Let’s check the variables inside the “geodata_dir” folder:
list.files(geodata_dir, recursive = TRUE)
#> [1] "climate/wc2.1_10m/wc2.1_10m_bioc_ACCESS-CM2_ssp126_2041-2060.tif"
#> [2] "climate/wc2.1_10m/wc2.1_10m_bioc_ACCESS-CM2_ssp126_2061-2080.tif"
#> [3] "climate/wc2.1_10m/wc2.1_10m_bioc_ACCESS-CM2_ssp585_2041-2060.tif"
#> [4] "climate/wc2.1_10m/wc2.1_10m_bioc_ACCESS-CM2_ssp585_2061-2080.tif"
#> [5] "climate/wc2.1_10m/wc2.1_10m_bioc_MIROC6_ssp126_2041-2060.tif"
#> [6] "climate/wc2.1_10m/wc2.1_10m_bioc_MIROC6_ssp126_2061-2080.tif"
#> [7] "climate/wc2.1_10m/wc2.1_10m_bioc_MIROC6_ssp585_2041-2060.tif"
#> [8] "climate/wc2.1_10m/wc2.1_10m_bioc_MIROC6_ssp585_2061-2080.tif"
#>
#> #Set climate as input directory
#> in_dir <- file.path(geodata_dir, "climate")
Now, we can organize and structure the files using the
organize_future_worldclim() function.
Static variables
When predicting to other times, some variables could be static (i.e.,
they remain unchanged in scenarios of projections). The
static_variables argument allows users to append static
variables alongside the Bioclimatic ones. First, let’s bring
soilType, which will remain static in future scenarios (we
will use it in a later step).
# Import raster layers (same used to calibrate and fit final models)
var <- rast(system.file("extdata", "Current_variables.tif", package = "kuenm2"))
# Get soilType
soiltype <- var$SoilType
Organizing files
Now, let’s organize the WorldClim files with the
organize_future_worldclim() function:
# Create folder to save structured files
out_dir_future <- file.path(tempdir(), "Future_raw") # a temporary directory
# Organize
organize_future_worldclim(input_dir = in_dir, # Path to variables from WorldClim
output_dir = out_dir_future,
name_format = "bio_", # Name format
static_variables = var$SoilType) # Static variables
#> | | | 0% | |========= | 12% | |================== | 25% | |========================== | 38% | |=================================== | 50% | |============================================ | 62% | |==================================================== | 75% | |============================================================= | 88% | |======================================================================| 100%
#>
#> Variables successfully organized in directory:
#> C:\Users\wever\AppData\Local\Temp\RtmpGMPg46/Future_raw
# Check files organized
dir(out_dir_future, recursive = TRUE)
#> [1] "2041-2060/ssp126/ACCESS-CM2/Variables.tif"
#> [2] "2041-2060/ssp126/MIROC6/Variables.tif"
#> [3] "2041-2060/ssp585/ACCESS-CM2/Variables.tif"
#> [4] "2041-2060/ssp585/MIROC6/Variables.tif"
#> [5] "2081-2100/ssp126/ACCESS-CM2/Variables.tif"
#> [6] "2081-2100/ssp126/MIROC6/Variables.tif"
#> [7] "2081-2100/ssp585/ACCESS-CM2/Variables.tif"
#> [8] "2081-2100/ssp585/MIROC6/Variables.tif"
We can check the files structured hierarchically in nested folders
using the dir_tree() function from the fs
package:
# Install package if necessary
if (!require("fs")) {
install.packages("fs")
}
dir_tree(out_dir_future)
#> Temp\RtmpkhmGWN/Future_raw
#> ├── 2041-2060
#> │ ├── ssp126
#> │ │ ├── ACCESS-CM2
#> │ │ │ └── Variables.tif
#> │ │ └── MIROC6
#> │ │ └── Variables.tif
#> │ └── ssp585
#> │ ├── ACCESS-CM2
#> │ │ └── Variables.tif
#> │ └── MIROC6
#> │ └── Variables.tif
#> └── 2081-2100
#> ├── ssp126
#> │ ├── ACCESS-CM2
#> │ │ └── Variables.tif
#> │ └── MIROC6
#> │ └── Variables.tif
#> └── ssp585
#> ├── ACCESS-CM2
#> │ └── Variables.tif
#> └── MIROC6
#> └── Variables.tif
After organizing variables, the next step is to create an object that
prepares things for projections.
Preparing data for projections
To prepare data for model projections across multiple scenarios,
storing the paths to the raster layers representing each scenario, we
use the function prepare_projection().
In contrast to predict_selected(), which requires a
SpatRaster object, when projecting to multiple scenarios,
we need the paths to the folders where the raster files are stored. This
includes the variables for the present time, which were used to
calibrate and fit the models. Currently, we only have the future climate
files. The present-day variables must reside in the same base directory
as the processed future variables. Let’s copy the raster layers used for
model fitting to a folder we can easily direct the function that runs
the next step:
# Create a "Current_raw" folder in a temporary directory
# and copy the rawvariables there.
out_dir_current <- file.path(tempdir(), "Current_raw")
dir.create(out_dir_current, recursive = TRUE)
# Save current variables in temporary directory
terra::writeRaster(var, file.path(out_dir_current, "Variables.tif"))
# Check folder
list.files(out_dir_current)
#> [1] "Variables.tif"
Now, we can prepare the data for projections (see the functions
documentation with help(prepare_projection)). In addition
to storing the paths to the variables for each scenario, the function
also verifies if all variables used to fit the final models are
available across all scenarios. To perform this check, you need to
provide either the fitted_models object you intend to use
for projection or simply the variable names. We strongly suggest using
the fitted_models object to prevent any errors.
We also need to define the root directory containing the scenarios
for projection (present, past, and/or future), along with additional
information regarding time periods, SSPs, and GCMs.
# Prepare projections using fitted models to check variables
pr <- prepare_projection(models = fitted_model_maxnet,
present_dir = out_dir_current, # Directory with present-day variables
past_dir = NULL, # NULL because we won't project to the past
past_period = NULL, # NULL because we won't project to the past
past_gcm = NULL, # NULL because we won't project to the past
future_dir = out_dir_future, # Directory with future variables
future_period = c("2041-2060", "2081-2100"),
future_pscen = c("ssp126", "ssp585"),
future_gcm = c("ACCESS-CM2", "MIROC6"))
The projection_data object summarizes information about
all the scenarios we will project to, and shows the root directory where
the variables are stored:
pr
#> projection_data object summary
#> =============================
#> Variables prepared to project models for Present and Future
#> Future projections contain the following periods, scenarios and GCMs:
#> - Periods: 2041-2060 | 2081-2100
#> - Scenarios: ssp126 | ssp585
#> - GCMs: ACCESS-CM2 | MIROC6
#> All variables are located in the following root directory:
#> C:/Users/wever/AppData/Local/Temp/RtmpGMPg46
If we check the structure of the prepared_projection
object, we can see it’s a list containing:
- Paths to all variables representing distinct scenarios in
subfolders.
- The pattern used to identify the format of raster files within the
folders (by default,
*.tif).
- The names of the variables.
- A list of class
prcomp if a Principal Component
Analysis (PCA) was performed on the set of variables with
prepare_data().
# Open prepared_projection in a new window
View(pr)
Projecting to multiple scenarios
After preparing the data, we can use the
project_selected() function to predict selected models
across the scenarios specified in prepare_projection:
# Create a folder to save projection results
# Here, in a temporary directory
out_dir <- file.path(tempdir(), "Projection_results/maxnet")
dir.create(out_dir, recursive = TRUE)
# Project selected models to multiple scenarios
p <- project_selected(models = fitted_model_maxnet,
projection_data = pr,
out_dir = out_dir,
write_replicates = TRUE,
progress_bar = FALSE) # Do not print progress bar
The function returns a model_projections object. This
object is similar to the prepared_data object, storing
information about the projections performed and the folder where results
were saved.
print(p)
#> model_projections object summary
#> ================================
#> Models projected for Present and Future
#> Future projections contain the following periods, scenarios and GCMs:
#> - Periods: 2041-2060 | 2081-2100
#> - Scenarios:
#> - GCMs: ACCESS-CM2 | MIROC6
#> All raster files containing the projection results are located in the following root directory:
#> C:/Users/wever/AppData/Local/Temp/RtmpGMPg46/Projection_results/maxnet
Note that results were saved hierarchically in nested subfolders,
each representing a distinct scenario. In the base directory, the
function also saves a file named “Projection_paths.RDS”, which is the
model_projections object. This object can be imported into
R using the readRDS() function.
dir_tree(out_dir)
#> Temp\Projection_results/maxnet
#> ├── Future
#> │ ├── 2041-2060
#> │ │ ├── ssp126
#> │ │ │ ├── ACCESS-CM2
#> │ │ │ │ ├── General_consensus.tif
#> │ │ │ │ ├── Model_192_consensus.tif
#> │ │ │ │ ├── Model_192_replicates.tif
#> │ │ │ │ ├── Model_219_consensus.tif
#> │ │ │ │ └── Model_219_replicates.tif
#> │ │ │ └── MIROC6
#> │ │ │ ├── General_consensus.tif
#> │ │ │ ├── Model_192_consensus.tif
#> │ │ │ ├── Model_192_replicates.tif
#> │ │ │ ├── Model_219_consensus.tif
#> │ │ │ └── Model_219_replicates.tif
#> │ │ └── ssp585
#> │ │ ├── ACCESS-CM2
#> │ │ │ ├── General_consensus.tif
#> │ │ │ ├── Model_192_consensus.tif
#> │ │ │ ├── Model_192_replicates.tif
#> │ │ │ ├── Model_219_consensus.tif
#> │ │ │ └── Model_219_replicates.tif
#> │ │ └── MIROC6
#> │ │ ├── General_consensus.tif
#> │ │ ├── Model_192_consensus.tif
#> │ │ ├── Model_192_replicates.tif
#> │ │ ├── Model_219_consensus.tif
#> │ │ └── Model_219_replicates.tif
#> │ └── 2081-2100
#> │ ├── ssp126
#> │ │ ├── ACCESS-CM2
#> │ │ │ ├── General_consensus.tif
#> │ │ │ ├── Model_192_consensus.tif
#> │ │ │ ├── Model_192_replicates.tif
#> │ │ │ ├── Model_219_consensus.tif
#> │ │ │ └── Model_219_replicates.tif
#> │ │ └── MIROC6
#> │ │ ├── General_consensus.tif
#> │ │ ├── Model_192_consensus.tif
#> │ │ ├── Model_192_replicates.tif
#> │ │ ├── Model_219_consensus.tif
#> │ │ └── Model_219_replicates.tif
#> │ └── ssp585
#> │ ├── ACCESS-CM2
#> │ │ ├── General_consensus.tif
#> │ │ ├── Model_192_consensus.tif
#> │ │ ├── Model_192_replicates.tif
#> │ │ ├── Model_219_consensus.tif
#> │ │ └── Model_219_replicates.tif
#> │ └── MIROC6
#> │ ├── General_consensus.tif
#> │ ├── Model_192_consensus.tif
#> │ ├── Model_192_replicates.tif
#> │ ├── Model_219_consensus.tif
#> │ └── Model_219_replicates.tif
#> ├── Present
#> │ └── Present
#> │ ├── General_consensus.tif
#> │ ├── Model_192_consensus.tif
#> │ ├── Model_192_replicates.tif
#> │ ├── Model_219_consensus.tif
#> │ └── Model_219_replicates.tif
#> └── Projection_paths.RDS
By default, separated by each scenario, the function computes
consensus metrics (mean, median, range, and standard deviation) for each
model across its replicates (if replicates were performed), as well as a
general consensus across all models (if more than one model was
selected).
By default, the function does not write these individual replicates
unless write_replicates = TRUE is set. It is important to
write the replicates if you intend to compute the variability deriving
from them in later steps (check Explore Variability and
Uncertainty in Projections).
The function project_selected() accepts several other
parameters that control how predictions are done, such as consensus to
compute, extrapolation type (free extrapolation (E),
extrapolation with clamping (EC), and no extrapolation
(NE)), variables to restrict, and the format of prediction
values (raw, cumulative,
logistic, or the default cloglog). For more
details, consult the vignette for Predict models to single scenario.
Importing results from projections
The model_projections object stores only the paths to
the resulting raster layers. To import the results, we can use the
import_results() function. By default, it imports all
consensus metrics (“median”, “range”, “mean”, and “stdev”) and all
scenarios (time periods, SSPs, and GCMs) available in the
model_projections object. Let’s import the mean for all
scenarios:
# Import mean of each projected scenario
p_mean <- import_results(projection = p, consensus = "mean")
# Plot all scenarios
terra::plot(p_mean, cex.main = 0.8)
Alternatively, we can import results from specific scenarios. For
example, let’s import the results only for the “2041-2060” time period
under the SSP 126:
# Importing
p_2060_ssp126 <- import_results(projection = p, consensus = "mean",
present = FALSE, # Do not import present projections
future_period = "2041-2060",
future_pscen = "ssp126")
# Plot all scenarios
terra::plot(p_2060_ssp126, cex.main = 0.8)
With the model_projections object, we can compute
changes in suitable areas (projection_changes()), explore
variability patterns coming from replicates, parameterizations, and GCMs
(projection_variability()), and perform analysis of
extrapolation risks (projection_mop()). For more details,
check Explore Variability and
Uncertainty in Projections.
Detecting changes in projections
When projecting a model to different scenarios (past or future),
suitable areas can be classified into three categories relative to the
current baseline: gain, loss and
stability. The interpretation of these categories
depends on the temporal direction of the projection.
When projecting to future scenarios:
- Gain: Areas currently unsuitable become suitable in the
future.
- Loss: Areas currently suitable become unsuitable in the
future.
- Stability: Areas remain the same in the future, suitable or
unsuitable.
When projecting to past scenarios:
- Gain: Areas unsuitable in the past are suitable in the
present.
- Loss: Areas suitable in the past are unsuitable in the
present.
- Stability: Areas remain the same in the present, suitable
or unsuitable.
These outcomes may vary across different General Circulation Models
(GCMs) within each time scenario (e.g., various Shared Socioeconomic
Pathways (SSPs) for the same period).
The projection_changes() function summarizes the number
of GCMs predicting gain, loss, and stability for each time scenario.
By default, this function writes the summary results to disk (unless
write_results is set to FALSE), but it does
not save binary layers for individual GCMs. In the example below, we
demonstrate how to configure the function to return the raster layers
with changes and write the binary results to disk.
# Run analysis to detect changes in suitable areas
changes <- projection_changes(model_projections = p,
output_dir = out_dir,
write_bin_models = TRUE, # Write individual binary results
return_raster = TRUE)
Setting colors for maps
Before plotting the results, we can use the
colors_for_changes() function to assign custom colors to
areas of gain, loss, and stability. By default, the function uses ‘teal
green’ for gains, ‘orange-red’ for losses, ‘Oxford blue’ for areas that
remain suitable, and ‘grey’ for areas that remain unsuitable. However,
you can customize these colors as needed.
The intensity of each color is automatically adjusted based on the
number of GCMs: highest (as defined by max_alpha) when all
GCMs agree on a prediction, and decreases progressively (down to
min_alpha) as fewer GCMs support that outcome.
# Set colors for change maps
changes_col <- colors_for_changes(changes)
The function returns the same changes_projections
object, but with color tables embedded in its SpatRasters.
These colors are automatically applied when visualizing the data using
terra::plot().
Types os results
The projection_changes() function returns four types of
results: binary model prediction, results by GCM, results by change, and
a general summary considering all GCMs:
- Binary prediction for each GCM: These are
suitable/unsuitable maps for each individual GCM. By default, the
omission error threshold used when selecting the best models is used
(e.g., 10%). You can specify a different threshold using the
user_threshold argument.
terra::plot(changes_col$Binarized, cex.main = 0.8)
- Results by GCM: provides the changes identified in
comparisons (gain, loss, stability) for each GCM individually.
terra::plot(changes_col$Results_by_gcm, cex.main = 0.8)
- Results by change: a list where each
SpatRaster represents a specific type of change (e.g.,
gain, loss, stability) across all GCMs for a given scenario.
# Results by change for the scenario of 2041-2060 (ssp126)
terra::plot(changes_col$Results_by_change$`Future_2041-2060_ssp126`)
- Summary changes: provides a general summary
indicating how many GCMs project gain, loss, and stability for each
scenario.
terra::plot(changes_col$Summary_changes, plg = list(cex = 0.75)) # Decrease size of legend text
Importing Results
When return_raster = TRUE is set, the resulting
SpatRaster objects are returned within the
changes_projections object. By default, however,
return_raster = FALSE and the object only contains the
directory path where the results were saved. In that case, the results
can be imported using the import_results() function. You
can specify the type of results to import, along with the target period
and emission scenario.
A changes_projections object imported using
import_results() can also be used as input to
colors_for_changes() to customize the colors used for
plotting.
For example, below we import only the general summary for the
2041–2060 period under the SSP5-8.5 scenario:
# Import changes detected for 2041–2060 SSP5-8.5
general_changes <- import_results(projection = changes,
future_period = "2041-2060",
future_pscen = "ssp585",
change_types = "summary")
# Set colors
general_changes <- colors_for_changes(general_changes)
# Plot
terra::plot(general_changes$Summary, main = names(general_changes$Summary),
plg = list(cex = 0.75)) # Decrease size of legend text
Saving results
The changes_projections object is a list that contains
the resulting SpatRaster objects (if
return_raster = TRUE) and the directory path where results
were saved (if write_results = TRUE).
If results were not written to disk during the initial run, you can
save SpatRaster objects afterward using the
writeRaster() function. For example, to save the general
summary raster:
writeRaster(changes$Summary_changes,
file.path(out_dir, "Summary_changes.tif"))
If the results were saved to disk, the
changes_projections object is automatically stored in a
folder named Projection_changes inside the specified
output_dir. You can load it back into R using
readRDS():
changes <- readRDS(file.path(out_dir, "Projection_changes/changes_projections.rds"))
After saving, this object can be used to import specific results with
the import_results() function.
# Reset plotting parameters
par(original_par)