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.

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.

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

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

Installing SyncroSim packages using 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

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

if (file.exists("stsimLibrary.ssim")){
  deleteLibrary("stsimLibrary.ssim", force = TRUE)
}
# Create a new library
myLibrary <- ssimLibrary(name = "stsimLibrary.ssim",
                         session = mySession,
                         packages = "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)
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.

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)

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)

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

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

Setting run parameters with 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)

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

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


syncrosim/rsyncrosim documentation built on Oct. 18, 2024, 1:29 a.m.