This vignette will cover how to develop and run spatially-explicit, stochastic state-and-transition simulation models of landscape change using the rsyncrosim
package within the SyncroSim software framework. For an overview of SyncroSim and rsyncrosim
, as well as a basic usage tutorial for rsyncrosim
, see the Introduction to rsyncrosim
vignette.
stsim
To demonstrate how to create and run spatially-explicit, stochastic state-and-transition simulation models (STSMs) of landscape change, we will be using the stsim SyncroSim package. STSMs are well-suited for integrating uncertainty into model projections, and have been applied to a variety of landscapes and management questions. In this stsim
example below, we will model changes in forest cover types under two scenarios: no harvest within a landscape, and harvest of 20 acres per year. To do this, we will set a number of model parameters describing how many timesteps and iterations the model will run for, transition types and their probabilities, transition targets (particularly, for the harvest transition type), etc. Next, we will provide stsim
with a set of initial conditions that describe the starting landscape at time 0, and specify parameters for the model's output. After running both scenarios, we will view the output data in tabular form and as a raster layer. Spatial and non-spatial examples of this exercise are provided.
For more details on stsim
, consult the ST-Sim online documentation.
Before using rsyncrosim
you will first need to download and install the SyncroSim software. Versions of SyncroSim exist for both Windows and Linux.
Note: this tutorial was developed using rsyncrosim
version 2.0. To use rsyncrosim
version 2.0 or greater, SyncroSim version 3.0 or greater is required.
You will need to install the rsyncrosim
R package, either using CRAN or from the rsyncrosim
GitHub repository. Versions of rsyncrosim
are available for both Windows and Linux. You may need to install the terra
package from CRAN as well.
In a new R script, load the necessary packages. This includes the rsyncrosim
and terra
R packages.
# Load R packages library(rsyncrosim) # package for working with SyncroSim library(terra) # package for working with spatial data
session()
Finish setting up the R environment for the rsyncrosim
workflow by creating a SyncroSim Session object. Use the session()
function to connect R to your installed copy of the SyncroSim software.
mySession <- session("path/to/install_folder") # Create a Session based SyncroSim install folder mySession <- session() # Using default install folder (Windows only) mySession # Displays the Session object
# Results of this code shown for above mySession <- session() # Using default install folder (Windows only) mySession # Displays the Session object
Use the version()
function to ensure you are using the latest version of SyncroSim.
version(mySession)
installPackage()
Install stsim
using the rsyncrosim
function installPackage()
. This function takes a package name as input and then queries the SyncroSim package server for the specified package.
installedPackages <- packages() if (is.element( "stsim", installedPackages$name)) uninstallPackage( "stsim")
# Install stsim installPackage("stsim")
stsim
should now be included in the package list returned by the packages()
function in rsyncrosim
:
# Get list of installed packages packages()
installedPackages <- packages() stsim_pkg <- installedPackages[installedPackages$name == "stsim", ] row.names(stsim_pkg) <- NULL stsim_pkg
When creating a new modeling workflow from scratch, we need to create objects of the following scopes:
For more information on these scopes, see the
Introduction to rsyncrosim
vignette.
if (file.exists("stsimLibrary.ssim")){ deleteLibrary("stsimLibrary.ssim", force = TRUE) }
# Create a new library myLibrary <- ssimLibrary(name = "stsimLibrary.ssim", session = mySession, packages = "stsim") # Open the default project myProject <- rsyncrosim::project(ssimObject = myLibrary, project = "Definitions") # Create a new scenario (associated with the default project) myScenario <- scenario(ssimObject = myProject, scenario = "My spatial scenario")
datasheet()
View the datasheets associated with your new project and scenario using the datasheet()
function from rsyncrosim
.
# View all datasheets associated with a library, project, or scenario datasheet_list <- datasheet(myScenario) head(datasheet_list) tail(datasheet_list)
From this list of datasheets, we can check which datasheets specific to the stsim
package we would like to modify. This will include a number of Initial Conditions datasheets, Output datasheets, and a Run Control datasheet. For more information about all stsim
datasheets, see the online documentation.
datasheet()
and addRow()
Now, we will add some values to these datasheets so we can run our models.
Inputs for stsim
come in two forms: scenario inputs (which can vary by scenario) and project inputs (which must be fixed across all scenarios). In general most inputs are scenario inputs; project inputs are typically reserved for values that must be shared by all scenarios (e.g. constants, shared lookup values). We refer to project input datasheets as project-scoped datasheets and scenario input datasheets as scenario-scoped datasheets. We will start by retrieving some project-scoped input datasheets to add and edit their values.
Terminology: The project-scoped datasheet called Terminology
specifies terms used across all scenarios in the same project.
# Load the Terminology datasheet to a new R data frame terminology <- datasheet(myProject, name = "stsim_Terminology") # Check the columns of the Terminology data frame str(terminology)
We can change the terminology of the StateLabelX
and AmountUnits
columns in this datasheet, and then save those changes back to the SyncroSim library file.
# Edit the values of the StateLabelX and AmountUnits columns terminology$AmountUnits <- "hectares" terminology$StateLabelX <- "Forest Type" # Saves edits as a SyncroSim datasheet saveDatasheet(ssimObject = myProject, data = terminology, name = "stsim_Terminology")
Similarly, we can edit other project-scoped datasheets for 'stsim'.
Stratum: Primary Strata in the model
# Load a copy of the Stratum datasheet. # To load an empty copy of this datasheet, specify the argument empty = TRUE stratum <- datasheet(myProject, "stsim_Stratum", empty = TRUE) # Use the addRow() to add a value to the stratum data frame stratum <- addRow(stratum, "Entire Forest") # Save edits as a SyncroSim datasheet saveDatasheet(myProject, stratum, "stsim_Stratum", force = TRUE)
StateLabelX: First dimension of labels for State Classes
It is also possible to add values to a datasheet without loading it into R. Below we create a vector of forestTypes
and add this to the stsim_StateLabelX
datasheet as a data frame using saveDatasheet()
.
# Create a vector containing the State Class labels forestTypes <- c("Coniferous", "Deciduous", "Mixed") # Add values as a data frame to a SyncroSim datasheet saveDatasheet(myProject, data.frame(Name = forestTypes), "stsim_StateLabelX", force = TRUE)
StateLabelY: Second dimension of labels for State Classes
# Add values as a data frame directly to an stsim datasheet saveDatasheet(myProject, data.frame(Name = c("All")), "stsim_StateLabelY", force = TRUE)
State Classes: Combine StateLabelX
and StateLabelY
, and assign each class a unique name and Id
# Create a new R data frame containing the names of the State Classes # and their corresponding data stateClasses <- data.frame(Name = forestTypes) stateClasses$StateLabelXId <- stateClasses$Name stateClasses$StateLabelYId <- "All" stateClasses$Id <- c(1, 2, 3) # Save stateClasses R data frame to a SyncroSim datasheet saveDatasheet(myProject, stateClasses, "stsim_StateClass", force = TRUE)
Transition Types: Assign a unique name and Id to each type of transition in our model.
# Create an R data frame containing transition type data transitionTypes <- data.frame(Name = c("Fire", "Harvest", "Succession"), Id = c(1, 2, 3)) # Save transitionTypes R data frame to a SyncroSim datasheet saveDatasheet(myProject, transitionTypes, "stsim_TransitionType", force = TRUE)
Transition Groups: Create Transition Groups identical to the Types
# Create an R data frame containing a column of transition type names transitionGroups <- data.frame(Name = c("Fire", "Harvest", "Succession")) # Save transitionGroups R data frame to a SyncroSim datasheet saveDatasheet(myProject, transitionGroups, "stsim_TransitionGroup", force = T)
Transition Types by Groups: Assign each Type to its Group
# Create an R data frame that contains Transition Type Group names transitionTypesGroups <- data.frame(TransitionTypeId = transitionTypes$Name, TransitionGroupId = transitionGroups$Name) # Save transitionTypesGroups R data frame to a SyncroSim datasheet saveDatasheet(myProject, transitionTypesGroups, "stsim_TransitionTypeGroup", force = T)
Ages: Define the basic parameters to control the age reporting in the model
# Define values for age reporting ageFrequency <- 1 ageMax <- 101 ageGroups <- c(20, 40, 60, 80, 100) # Add values as R data frames to the appropriate SyncroSim datasheet saveDatasheet(myProject, data.frame(Frequency = ageFrequency, MaximumAge = ageMax), "stsim_AgeType", force = TRUE) saveDatasheet(myProject, data.frame(MaximumAge = ageGroups), "stsim_AgeGroup", force = TRUE)
Now that we have defined all our project-scoped input datasheets, we can move on to specifying scenario-specific model inputs. We begin by using the scenario
function to create a new scenario in our project.
# Create a new SyncroSim scenario myScenario <- scenario(myProject, "No Harvest")
Once again we can use the datasheet()
function (with summary=TRUE
) to display all the scenario-scoped datasheets.
# Subset the full datasheet list to show only scenario-scoped datasheets scenario_datasheet_list <- subset(datasheet(myScenario, summary = TRUE), scope == "scenario") head(scenario_datasheet_list) tail(scenario_datasheet_list)
We can now use the datasheet()
function to retrieve, one at a time, each of our scenario-scoped datasheets from our library.
Run Control: Define the length of the run and whether or not it is a spatial run (requires spatial inputs to be set, see below). Here we make the run spatial.
# Create an R data frame specifying to run the simulation for # 7 realizations and 10 timesteps runControl <- data.frame(MaximumIteration = 7, MinimumTimestep = 0, MaximumTimestep = 10, IsSpatial = TRUE) # Save transitionTypesGroups R data frame to a SyncroSim datasheet saveDatasheet(myScenario, runControl, "stsim_RunControl",append = FALSE)
Deterministic Transitions: Define transitions that take place in the absence of probabilistic transitions. Here we also set the age boundaries for each State Class.
# Load an empty Deterministic Transitions datasheet to a new R data frame dTransitions <- datasheet(myScenario, "stsim_DeterministicTransition", optional = T, empty = T, lookupsAsFactors = F) # Add all Deterministic Transitions to the R data frame dTransitions <- addRow(dTransitions, data.frame( StateClassIdSource = "Coniferous", StateClassIdDest = "Coniferous", AgeMin = 21, Location = "C1")) dTransitions <- addRow(dTransitions, data.frame( StateClassIdSource = "Deciduous", StateClassIdDest = "Deciduous", Location = "A1")) dTransitions <- addRow(dTransitions, data.frame( StateClassIdSource = "Mixed", StateClassIdDest = "Mixed", AgeMin = 11, Location = "B1")) # Save dTransitions R data frame to a SyncroSim datasheet saveDatasheet(myScenario, dTransitions, "stsim_DeterministicTransition")
Probabilistic Transitions: Define transitions between State Classes and assigns a probability to each.
# Load an empty Probabilistic Transitions datasheet to a new R data frame pTransitions <- datasheet(myScenario, "stsim_Transition", optional = T, empty = T, lookupsAsFactors = F) # Add all Probabilistic Transitions to the R data frame pTransitions <- addRow(pTransitions, data.frame( StateClassIdSource = "Coniferous", StateClassIdDest = "Deciduous", TransitionTypeId = "Fire", Probability = 0.01)) pTransitions <- addRow(pTransitions, data.frame( StateClassIdSource = "Coniferous", StateClassIdDest = "Deciduous", TransitionTypeId = "Harvest", Probability = 1, AgeMin = 40)) pTransitions <- addRow(pTransitions, data.frame( StateClassIdSource = "Deciduous", StateClassIdDest = "Deciduous", TransitionTypeId = "Fire", Probability = 0.002)) pTransitions <- addRow(pTransitions, data.frame( StateClassIdSource = "Deciduous", StateClassIdDest = "Mixed", TransitionTypeId = "Succession", Probability = 0.1, AgeMin = 10)) pTransitions <- addRow(pTransitions, data.frame( StateClassIdSource = "Mixed", StateClassIdDest = "Deciduous", TransitionTypeId = "Fire", Probability = 0.005)) pTransitions <- addRow(pTransitions, data.frame( StateClassIdSource = "Mixed", StateClassIdDest = "Coniferous", TransitionTypeId = "Succession", Probability = 0.1, AgeMin = 20)) # Save pTransitions R data frame to a SyncroSim datasheet saveDatasheet(myScenario, pTransitions, "stsim_Transition")
Initial Conditions: Set the starting conditions of the model at time 0. There are two options for setting initial conditions: either spatial or non-spatial. In this example we will use spatial initial conditions; however, we also demonstrate how to set initial conditions non-spatially below.
rsyncrosim
package. # Load sample .tif files stratumTif <- "C:/gitprojects/rsyncrosim/vignettes/initial-stratum.tif" sclassTif <- "C:/gitprojects/rsyncrosim/vignettes/initial-sclass.tif" ageTif <- "C:/gitprojects/rsyncrosim/vignettes/initial-age.tif"
# Load sample .tif files stratumTif <- "initial-stratum.tif" sclassTif <- "initial-sclass.tif" ageTif <- "initial-age.tif"
# Create raster layers from the .tif files rStratum <- rast(stratumTif) rSclass <- rast(sclassTif) rAge <- rast(ageTif) # Plot raster layers plot(rStratum) plot(rSclass) plot(rAge)
We can add these rasters as model inputs using the stsim_InitialConditionsSpatial
datasheet.
# Create an data.frame of the input raster layers ICSpatial <- data.frame(StratumFileName = stratumTif, StateClassFileName = sclassTif, AgeFileName = ageTif) # Save initialConditionsSpatial R list to a SyncroSim datasheet saveDatasheet(myScenario, ICSpatial, "stsim_InitialConditionsSpatial")
stsim_InitialConditionsNonSpatial
and stsim_InitialConditionsNonSpatialDistribution
datasheets:# Create non-spatial initial conditions data and add it to an R data frame ICNonSpatial <- data.frame(TotalAmount = 100, NumCells = 100, CalcFromDist = F) # Save the ICNonSpatial R data frame to a SyncroSim datasheet saveDatasheet(myScenario, ICNonSpatial, "stsim_InitialConditionsNonSpatial")
# Create non-spatial initial conditions distribution data and add it to an R data frame ICNonSpatialDistribution <- data.frame(StratumId = "Entire Forest", StateClassId = "Coniferous", RelativeAmount = 1) # Save the ICNonSpatial R data frame to a SyncroSim datasheet saveDatasheet(myScenario, ICNonSpatialDistribution, "stsim_InitialConditionsNonSpatialDistribution")
Transition Targets: Define targets, in units of area, to be reached by the allocation procedure within SyncroSim.
# Set the transition target for harvest to 0 saveDatasheet(myScenario, data.frame(TransitionGroupId = "Harvest", Amount = 0), "stsim_TransitionTarget")
Output Options: Regulate the model outputs and determine the frequency at which syncrosim saves the model outputs.
# Create output options for spatial model and add it to an R data frame outputOptionsSpatial <- data.frame( RasterOutputSC = T, RasterOutputSCTimesteps = 1, RasterOutputTR = T, RasterOutputTRTimesteps = 1, RasterOutputAge = T, RasterOutputAgeTimesteps = 1 ) # Save the outputOptionsSpatial R data frame to a SyncroSim datasheet saveDatasheet(myScenario, outputOptionsSpatial, "stsim_OutputOptionsSpatial") # Create output options for non-spatial model and add it to an R data frame outputOptionsNonSpatial <- data.frame( SummaryOutputSC = T, SummaryOutputSCTimesteps = 1, SummaryOutputTR = T, SummaryOutputTRTimesteps = 1 ) # Save the outputOptionsNonSpatial R data frame to a SyncroSim datasheet saveDatasheet(myScenario, outputOptionsNonSpatial, "stsim_OutputOptions")
Pipeline datasheet: Specify which transformer stage the scenarios will run and in which order.
To learn more about pipelines, see the rsyncrosim
: introduction to pipelines vignette and the SyncroSim Enhancing a Package: Linking Models tutorial.
To implement pipelines in our package, we need to specify the order in which to run the transformers (i.e. models) in our pipeline by editing the Pipeline
datasheet. The Pipeline
datasheet is part of the built-in SyncroSim core, so we access it using the "core_" prefix with the datasheet()
function.
From viewing the structure of the Pipeline
datasheet we know that the StageNameId
is a factor with a single level: ST-Sim
.
# Load Pipeline datasheet to an R data frame myPipelineDataframe <- datasheet(myScenario, name = "core_Pipeline") # Check the columns of the Pipeline data frame str(myPipelineDataframe) # Create Pipeline data and add it to the Pipeline data frame myPipelineRow <- data.frame(StageNameId = c("ST-Sim"), RunOrder = c(1)) myPipelineDataframe <- addRow(myPipelineDataframe, myPipelineRow) # Check values myPipelineDataframe # Save Pipeline R data frame to a SyncroSim datasheet saveDatasheet(ssimObject = myScenario, data = myPipelineDataframe, name = "core_Pipeline")
We are done parameterizing our simple "No Harvest" scenario. Let's now define a new scenario that implements forest harvesting. Below, we create a second "Harvest" scenario that is a copy of the first scenario, but with a harvest level of 20 acres/year.
# Create a copy of the no harvest scenario (i.e myScenario) and name it myScenarioHarvest myScenarioHarvest <- scenario(myProject, scenario = "Harvest", sourceScenario = myScenario) # Set the transition target for harvest to 20 acres/year saveDatasheet(myScenarioHarvest, data.frame(TransitionGroupId = "Harvest", Amount = 20), "stsim_TransitionTarget")
We can display the harvest levels for both scenarios.
# View the transition targets for the Harvest and No Harvest scenarios datasheet(myProject, scenario = c("Harvest", "No Harvest"), name = "stsim_TransitionTarget")
run()
We will now run both scenarios using the run()
function in rsyncrosim
.
If we have a large model and we want to parallelize the run using multiprocessing, we can modify the library-scoped "core_Multiprocessing" datasheet. Since we are using five iterations in our model, we will set the number of jobs to five so each multiprocessing core will run a single iteration.
# Load list of available library-scoped datasheets datasheet(myLibrary) # Load the library-scoped multiprocessing datasheet multiprocess <- datasheet(myLibrary, name = "core_Multiprocessing") # Check required inputs str(multiprocess) # Enable multiprocessing multiprocess$EnableMultiprocessing <- TRUE # Set maximum number of jobs to 5 multiprocess$MaximumJobs <- 5 # Save multiprocessing configuration saveDatasheet(ssimObject = myLibrary, data = multiprocess, name = "core_Multiprocessing")
Now, when we run our scenario, it will use the desired multiprocessing configuration. Running a scenario generates a corresponding new child scenario, called a results scenario, which contains the results of the run along with a snapshot of all the model inputs.
# Run both scenarios myResultScenario <- run(myProject, scenario = c("Harvest", "No Harvest"), summary = TRUE)
The next step is to view the output datasheets added to the result scenario when it was run. To look at the results we first need to retrieve the unique scenarioId
for each child result scenario.
# Retrieve scenario IDs resultIDNoHarvest <- subset(myResultScenario, ParentId == scenarioId(myScenario))$ScenarioId resultIDHarvest <- subset(myResultScenario, ParentId == scenarioId(myScenarioHarvest))$ScenarioId
We can now retrieve tabular output regarding the projected State Class over time (for both scenarios combined) from the stsim_OutputStratumState
datasheet.
# Retrieve output projected State Class for both Scenarios in tabular form outputStratumState <- datasheet( myProject, scenario = c(resultIDNoHarvest, resultIDHarvest), name = "stsim_OutputStratumState")
Finally, we can get the State Class raster output using the datasheetRaster()
function (here for the Harvest scenario only).
# Retrieve the output State Class raster for the Harvest scenario at timestep 5 myRastersTimestep5 <- datasheetSpatRaster(ssimObject = myProject, scenario = resultIDHarvest, datasheet = "stsim_OutputSpatialState", timestep = 5) myRastersTimestep5 # Plot raster for the first realization of timestep 5 plot(myRastersTimestep5[[1]])
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.