`rsyncrosim`: introduction to `ST-Sim`

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.

SyncroSim Package: `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.

Setup

Install SyncroSim

Before using rsyncrosim you will first need to download and install the SyncroSim software. Versions of SyncroSim exist for both Windows and Linux.

Installing and loading R packages

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 raster and rgdal packages from CRAN as well.

In a new R script, load the necessary packages. This includes the rsyncrosim and raster R packages.

# Load R packages
library(rsyncrosim)  # package for working with SyncroSim

## Warning: package 'rsyncrosim' was built under R version 4.3.3

library(terra)       # package for working with spatial data

Connecting R to SyncroSim using `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

## class               : Session
## filepath [character]: C:/Program Files/SyncroSim
## silent [logical]    : TRUE
## printCmd [logical]  : FALSE
## condaFilepath [NULL]:

Use the version() function to ensure you are using the latest version of SyncroSim.

version(mySession)

## [1] "2.5.11"

Installing SyncroSim packages using `addPackage()`

Install stsim using the rsyncrosim function addPackage(). This function takes a package name as input and then queries the SyncroSim package server for the specified package.

# Install stsim
addPackage("stsim")

## Package <stsim> installed

stsim should now be included in the package list returned by the package() function in rsyncrosim:

# Get list of installed packages
package()

##    name                                      description version
## 1 stsim The ST-Sim state-and-transition simulation model   3.4.5
##                                       location status
## 1 C:\\Users\\sarah\\SyncroSim\\Packages\\stsim     OK

Create a modeling workflow

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.

Set up Library, Project, and Scenario

# Create a new Library
myLibrary <- ssimLibrary(name = "stsimLibrary.ssim",
                         session = mySession,
                         package = "stsim",
                         overwrite = TRUE)

# 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")

View model inputs using `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)

##     scope                 name              displayName
## 1 library          core_Backup                   Backup
## 2 library     core_CondaConfig      Conda Configuration
## 3 library        core_JlConfig      Julia Configuration
## 4 library      core_LNGPackage Last Known Good Packages
## 5 library core_Multiprocessing          Multiprocessing
## 6 library         core_Options                  Options

tail(datasheet_list)

##        scope                                        name
## 107 scenario             stsim_TransitionSlopeMultiplier
## 108 scenario stsim_TransitionSpatialInitiationMultiplier
## 109 scenario           stsim_TransitionSpatialMultiplier
## 110 scenario          stsim_TransitionSpreadDistribution
## 111 scenario                      stsim_TransitionTarget
## 112 scenario        stsim_TransitionTargetPrioritization
##                                   displayName
## 107              Transition Slope Multipliers
## 108 Transition Spatial Initiation Multipliers
## 109            Transition Spatial Multipliers
## 110            Transition Spread Distribution
## 111                        Transition Targets
## 112          Transition Target Prioritization

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.

Configure model inputs using `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.

Project-scoped Datasheets

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)

## 'data.frame':    1 obs. of  8 variables:
##  $ AmountLabel          : chr "Area"
##  $ AmountUnits          : Factor w/ 5 levels "acres","hectares",..: 1
##  $ StateLabelX          : chr "Class"
##  $ StateLabelY          : chr "Subclass"
##  $ PrimaryStratumLabel  : chr "Primary Stratum"
##  $ SecondaryStratumLabel: chr "Secondary Stratum"
##  $ TertiaryStratumLabel : chr "Tertiary Stratum"
##  $ TimestepUnits        : chr "Timestep"

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(myProject, terminology, "stsim_Terminology")

## Datasheet <stsim_Terminology> saved

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 in the datasheet() function.
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)

## Datasheet <stsim_Stratum> saved

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)

## Datasheet <stsim_StateLabelX> saved

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)

## Datasheet <stsim_StateLabelY> saved

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)

## Datasheet <stsim_StateClass> saved

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)

## Datasheet <stsim_TransitionType> saved

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, "TransitionGroup", force = T)

## Datasheet <stsim_TransitionGroup> saved

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, 
              "TransitionTypeGroup", 
              force = T)

## Datasheet <stsim_TransitionTypeGroup> saved

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),
              "AgeType", 
              force = TRUE)

## Datasheet <stsim_AgeType> saved

saveDatasheet(myProject, 
              data.frame(MaximumAge = ageGroups), 
              "stsim_AgeGroup", 
              force = TRUE)

## Datasheet <stsim_AgeGroup> saved

Scenario-scoped Datasheets

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)

##       scope                            name                displayName
## 46 scenario            core_AutoGenTagValue Auto Generation Tag Values
## 47 scenario                   core_Pipeline                   Pipeline
## 48 scenario     corestime_DistributionValue              Distributions
## 49 scenario              corestime_External                   External
## 50 scenario corestime_ExternalVariableValue         External Variables
## 51 scenario       corestime_Multiprocessing    Spatial Multiprocessing

tail(scenario_datasheet_list)

##        scope                                        name
## 107 scenario             stsim_TransitionSlopeMultiplier
## 108 scenario stsim_TransitionSpatialInitiationMultiplier
## 109 scenario           stsim_TransitionSpatialMultiplier
## 110 scenario          stsim_TransitionSpreadDistribution
## 111 scenario                      stsim_TransitionTarget
## 112 scenario        stsim_TransitionTargetPrioritization
##                                   displayName
## 107              Transition Slope Multipliers
## 108 Transition Spatial Initiation Multipliers
## 109            Transition Spatial Multipliers
## 110            Transition Spread Distribution
## 111                        Transition Targets
## 112          Transition Target Prioritization

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)

## Datasheet <stsim_RunControl> saved

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)

# 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")

## Datasheet <stsim_DeterministicTransition> saved

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)

# 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")

## Datasheet <stsim_Transition> saved

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.

Option 1 - Spatial. Let’s first take a look at our rasters. We will retrieve the sample raster files (in GeoTIFF format) provided with the rsyncrosim package.

# 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 R list of the input raster layers
ICSpatial <- list(StratumFileName = stratumTif, 
                  StateClassFileName = sclassTif, 
                  AgeFileName = ageTif)

# Save initialConditionsSpatial R list to a SyncroSim Datasheet
saveDatasheet(myScenario, ICSpatial, "stsim_InitialConditionsSpatial")

## Datasheet <stsim_InitialConditionsSpatial> saved

Option 2 - Non-spatial. The second option is to set the proportions of each class, making this a non-spatial parameterization. To do so, we use the 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")

## Datasheet <stsim_InitialConditionsNonSpatial> saved

# 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")

## Datasheet <stsim_InitialConditionsNonSpatialDistribution> saved

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")

## Datasheet <stsim_TransitionTarget> saved

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")

## Datasheet <stsim_OutputOptionsSpatial> saved

# 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")

## Datasheet <stsim_OutputOptions> saved

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")

## Datasheet <stsim_TransitionTarget> saved

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")

##   ScenarioID ProjectID ScenarioName ParentID ParentName TransitionGroupID
## 1          2         1   No Harvest       NA       <NA>           Harvest
## 2          3         1      Harvest       NA       <NA>           Harvest
##   Amount
## 1      0
## 2     20

Run Scenarios

Setting run parameters with `run()`

We will now run both Scenarios using the run() function in rsyncrosim. 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. If we have a large modeling workflow and we want to parallelize the run using multiprocessing, we can set the jobs argument to be a value greater than one.

# Run both Scenarios
myResultScenario <- run(myProject, scenario = c("Harvest", "No Harvest"), 
                     jobs = 7, summary = TRUE)

## [1] "Running scenario [3] Harvest"
## [1] "Running scenario [2] No Harvest"

View results

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

## class       : SpatRaster 
## dimensions  : 10, 10, 7  (nrow, ncol, nlyr)
## resolution  : 1, 1  (x, y)
## extent      : 0, 10, 0, 10  (xmin, xmax, ymin, ymax)
## coord. ref. : UTM Zone 10, Northern Hemisphere 
## sources     : sc.it1.ts5.tif  
##               sc.it2.ts5.tif  
##               sc.it3.ts5.tif  
##               ... and 4 more source(s)
## names       : sc.it1.ts5, sc.it2.ts5, sc.it3.ts5, sc.it4.ts5, sc.it5.ts5, sc.it6.ts5, ...

# Plot raster for the first realization of timestep 5
plot(myRastersTimestep5[[1]])

rsyncrosim: introduction to ST-Sim

SyncroSim Package: stsim