ST-Sim in rsyncrosim
vignettes/rsyncrosim_stsim_vignette.Rmd
rsyncrosim_stsim_vignette.Rmdrsyncrosim?rsyncrosim is a package designed to facilitate the development of modeling workflows in R for the SyncroSim software framework.
SyncroSim is a software platform that helps you turn your data into forecasts. At the core of SyncroSim is an engine that automatically structures your existing data, regardless of its original format. SyncroSim then feeds this structured data through a pipeline of calculations (i.e. a suite of models) in order to transform the data into forecasts. As you develop alternative “what-if” forecasting scenarios, SyncroSim keeps track of all of your model inputs and outputs for you. Finally, SyncroSim provides a rich interface for you and others to interact with your data and models, in order to explore and track the consequences of alternative “what-if” forecasting scenarios.
For more details consult the SyncroSim online documentation.
rsyncrosim
Before using rsyncrosim you will first need to download and install the SyncroSim software. Versions of SyncroSim exist for both Windows and Linux.
Next 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.
The following tutorial introduces rsyncrosim using the stsim SyncroSim Package. st-sim is a SyncroSim Package that lets you develop and run spatially-explicit, stochastic state-and-transition simulation models of landscape change. For more details on stsim consult the ST-Sim online documentation.
Let’s first load the necessary packages. This includes both rsyncrosim and the raster R package.
For this tutorial, we will first need to retrieve the sample raster files (in GeoTIFF format) provided with the rsyncrosim package.
stratumTif <- system.file("extdata", "initial-stratum.tif",
package = "rsyncrosim")
sclassTif <- system.file("extdata", "initial-sclass.tif",
package = "rsyncrosim")
ageTif <- system.file("extdata", "initial-age.tif",
package = "rsyncrosim")The first step in a typical rsyncrosim workflow is to create a SyncroSim Session object in R that provides the connection to your installed copy of the SyncroSim software. A new Session is created using the session function, for which the first argument is a path to the folder on your computer into which the SyncroSim software has been installed.
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 objectThe next step is to create a SyncroSim Library. A Library is a file (with .ssim extension) that stores all of your model inputs and outputs. The format of each SyncroSim Library is unique to the SyncroSim Package with which it is associated. We use the ssimLibrary function to create a new SsimLibrary object in R that is connected (through your Session) to a SyncroSim Library file. Note that an existing Library file can be later opened again by setting overwrite=FALSE for the ssimLibrary function.
addPackage("stsim", mySession) # Get stsim package from SyncroSim online server (only the first time)
myLibrary <- ssimLibrary(name = "demoLibrary.ssim", session = mySession, overwrite = TRUE)Each SyncroSim Library contains one or more SyncroSim Projects, each represented by a Project object in R. Projects typically store model inputs that are common to all your scenarios. In most situations you will need only a single Project for your Library; by default each new Library starts with a single Project named “Definitions” (with a unique projectId= 1). The project function is used to both create and retrieve Projects.
# Loads the existing default Project
myProject = project(myLibrary, project = "Definitions") # Using name for Project
myProject = project(myLibrary, project = 1) # Using projectId for ProjectFinally, each SyncroSim Project contains one or more Scenarios, each represented by a Scenario object in R. Scenarios store the specific model inputs and outputs associated with each model run in SyncroSim. Each Scenario can be identified by its unique scenarioId. The scenario function is used to both create and retrieve Scenarios.
# Creates a Scenario object (associated with the default Project)
myScenario = scenario(myProject, scenario = "My first scenario")
myScenarioFinally, each SyncroSim Library contains multiple SyncroSim Datasheets. A SyncroSim Datasheet is simply a table of data stored in the Library. Datasheets each have a scope: either library, project or scenario.
For example use the datasheet function (with the argument summary=TRUE) to list all the available SyncroSim Datasheets for your Library. These typically represent general settings for your SyncroSim Library.
datasheet(myLibrary, summary = TRUE)Note that each datasheet has a name. This is important - you will use this later to retrieve the data associated with a specific Datasheet.
Similarly you can view the list of Datasheets associated with your Project and your Scenario
In order to run a model in SyncroSim, we first need to setup all of the required model inputs. As outlined above, inputs 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).
Let’s begin by configuring our Project inputs. All inputs in SyncroSim are stored in Datasheets. We say that inputs associated with a Project are project-scoped Datasheets. For example, in stsim there is a project-scoped Datasheet called Terminology specifying terms used across all Scenarios. We will use the datasheet function to retrieve the current values for this Datasheet, as currently stored in the SyncroSim Library file.
# Returns a dataframe of the Terminology Datasheet
sheetData <- datasheet(myProject, name="stsim_Terminology")
str(sheetData)We can now change the terminology of the StateLabelX and AmountUnits columns in our Datasheet, and then save those changes back to the SyncroSim Library file.
# Edit the dataframe values
sheetData$AmountUnits <- "Hectares"
sheetData$StateLabelX <- "Forest Type"
# Saves the edits back to the Library file
saveDatasheet(myProject, sheetData, "stsim_Terminology") Similarly we can edit other Project-scoped Datasheets for stsim.
Stratum: Primary Strata in the model
# Retrieves an empty copy of the Datasheet
sheetData <- datasheet(myProject, "stsim_Stratum", empty = TRUE)
# Helper function in rsyncrosim to add rows to dataframes
sheetData <- addRow(sheetData, "Entire Forest")
# Save edits to file
saveDatasheet(myProject, sheetData, "stsim_Stratum", force = TRUE)StateLabelX: First dimension of labels for State Classes
forestTypes <- c("Coniferous", "Deciduous", "Mixed")
saveDatasheet(myProject, data.frame(Name = forestTypes), "stsim_StateLabelX", force = TRUE)StateLabelY: Second dimension of labels for State Classes
saveDatasheet(myProject, data.frame(Name = c("All")), "stsim_StateLabelY", force = TRUE)State Classes: Combines StateLabelX and StateLabelY; assigns each one a unique name and ID
stateClasses <- data.frame(Name = forestTypes)
stateClasses$StateLabelXID <- stateClasses$Name
stateClasses$StateLabelYID <- "All"
stateClasses$ID <- c(1, 2, 3)
saveDatasheet(myProject, stateClasses, "stsim_StateClass", force = TRUE)Transition Types: Defines the types of transition in our model; assigns each one a unique name and ID
transitionTypes <- data.frame(Name = c("Fire", "Harvest", "Succession"), ID = c(1, 2, 3))
saveDatasheet(myProject, transitionTypes, "stsim_TransitionType", force = TRUE)Transition Groups: Makes the Groups identical to the Types
transitionGroups <- data.frame(Name = c("Fire", "Harvest", "Succession"))
saveDatasheet(myProject, transitionGroups, "TransitionGroup", force = T)Transition Types by Groups: Assigns each Type to its Group
transitionTypesGroups <- data.frame(TransitionTypeID = transitionTypes$Name,
TransitionGroupID = transitionGroups$Name)
saveDatasheet(myProject, transitionTypesGroups, "TransitionTypeGroup", force = T)Ages defines the basic parameters to control the age reporting in the model
ageFrequency <- 1
ageMax <- 101
ageGroups <- c(20, 40, 60, 80, 100)
saveDatasheet(myProject, data.frame(Frequency = ageFrequency, MaximumAge = ageMax),
"AgeType", force = TRUE)
saveDatasheet(myProject, data.frame(MaximumAge = ageGroups), "stsim_AgeGroup", force = TRUE)Now that we have defined all our Project level definitions 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.
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 list to show only Scenario Datasheets
subset(datasheet(myScenario, summary = TRUE), scope == "scenario")We can now use the datasheet function to retrieve, one at a time, each of our Scenario-scoped Datasheets from our Library.
Run Control defines 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:
# Run simulation for 7 realizations and 10 timesteps
sheetName <- "stsim_RunControl"
sheetData <- data.frame(MaximumIteration = 7,
MinimumTimestep = 0,
MaximumTimestep = 10,
isSpatial = TRUE)
saveDatasheet(myScenario, sheetData, sheetName)Deterministic Transitions first define transitions that take place in the absence of probabilistic transitions. Here we also set the age boundaries for each State Class:
# Add all the deterministic transitions
sheetName <- "stsim_DeterministicTransition"
sheetData <- datasheet(myScenario, sheetName, optional = T, empty = T)
sheetData <- addRow(sheetData, data.frame(StateClassIDSource = "Coniferous",
StateClassIDDest = "Coniferous",
AgeMin = 21,
Location = "C1"))
sheetData <- addRow(sheetData, data.frame(StateClassIDSource = "Deciduous",
StateClassIDDest = "Deciduous",
Location = "A1"))
sheetData <- addRow(sheetData, data.frame(StateClassIDSource = "Mixed",
StateClassIDDest = "Mixed",
AgeMin = 11,
Location = "B1"))
saveDatasheet(myScenario, sheetData, sheetName)Probabilistic Transitions define the transitions between State Classes and assigns a probability to each:
# Add all the probabilistic transition pathways
sheetName <- "stsim_Transition"
sheetData <- datasheet(myScenario, sheetName, optional = T, empty = T)
sheetData <- addRow(sheetData, data.frame(StateClassIDSource = "Coniferous",
StateClassIDDest = "Deciduous",
TransitionTypeID = "Fire",
Probability = 0.01))
sheetData <- addRow(sheetData, data.frame(StateClassIDSource = "Coniferous",
StateClassIDDest = "Deciduous",
TransitionTypeID = "Harvest",
Probability = 1,
AgeMin = 40))
sheetData <- addRow(sheetData, data.frame(StateClassIDSource = "Deciduous",
StateClassIDDest = "Deciduous",
TransitionTypeID = "Fire",
Probability = 0.002))
sheetData <- addRow(sheetData, data.frame(StateClassIDSource = "Deciduous",
StateClassIDDest = "Mixed",
TransitionTypeID = "Succession",
Probability = 0.1,
AgeMin = 10))
sheetData <- addRow(sheetData, data.frame(StateClassIDSource = "Mixed",
StateClassIDDest = "Deciduous",
TransitionTypeID = "Fire",
Probability = 0.005))
sheetData <- addRow(sheetData, data.frame(StateClassIDSource = "Mixed",
StateClassIDDest = "Coniferous",
TransitionTypeID = "Succession",
Probability = 0.1,
AgeMin = 20))
saveDatasheet(myScenario, sheetData, sheetName)Initial Conditions sets 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 demonstrate below how also to set initial conditions non-spatially.
rStratum <- raster(stratumTif)
rSclass <- raster(sclassTif)
rAge <- raster(ageTif)
plot(rStratum)
plot(rSclass)
plot(rAge)We can add these rasters as model inputs using the stsim_InitialConditionsSpatial datasheet.
sheetName <- "stsim_InitialConditionsSpatial"
sheetData <- list(StratumFileName = stratumTif,
StateClassFileName = sclassTif,
AgeFileName = ageTif)
saveDatasheet(myScenario, sheetData, sheetName)Let’s check if the rasters were entered correctly. We can extract rasters with datasheetRaster function.
rStratumTest <- datasheetRaster(myScenario, sheetName, "StratumFileName")
rSclassTest <- datasheetRaster(myScenario, sheetName, "StateClassFileName")
rAgeTest <- datasheetRaster(myScenario, sheetName, "AgeFileName")
plot(rStratumTest)
plot(rSclassTest)
plot(rAgeTest)Initial Conditions: 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:
sheetName <- "stsim_InitialConditionsNonSpatial"
sheetData <- data.frame(TotalAmount = 100,
NumCells = 100,
CalcFromDist = F)
saveDatasheet(myScenario, sheetData, sheetName)
datasheet(myScenario, sheetName)
sheetName <- "stsim_InitialConditionsNonSpatialDistribution"
sheetData <- data.frame(StratumID = "Entire Forest",
StateClassID = "Coniferous",
RelativeAmount = 1)
saveDatasheet(myScenario, sheetData, sheetName)
datasheet(myScenario, sheetName)Transition Targets defines targets, in units of area, to be reached by the allocation procedure within SyncroSim:
# Transition targets - set harvest to 0 for this scenario
saveDatasheet(myScenario,
data.frame(TransitionGroupID = "Harvest",
Amount = 0),
"stsim_TransitionTarget")Output Options regulates the model outputs and determines the frequency at which syncrosim saves the model outputs:
# Output options
# datasheet(myScenario, "stsim_OutputOptions")
sheetData <- data.frame(
SummaryOutputSC = T, SummaryOutputSCTimesteps = 1,
SummaryOutputTR = T, SummaryOutputTRTimesteps = 1
)
saveDatasheet(myScenario, sheetData, "stsim_OutputOptions")
sheetData <- data.frame(
RasterOutputSC = T, RasterOutputSCTimesteps = 1,
RasterOutputTR = T, RasterOutputTRTimesteps = 1,
RasterOutputAge = T, RasterOutputAgeTimesteps = 1
)
saveDatasheet(myScenario, sheetData, "stsim_OutputOptionsSpatial")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.
myScenarioHarvest <- scenario(myProject,
scenario = "Harvest",
sourceScenario = myScenario)
saveDatasheet(myScenarioHarvest, data.frame(TransitionGroupID = "Harvest",
Amount = 20),
"stsim_TransitionTarget")We can display the harvest levels for both scenarios.
We can now run the Scenario and look at the results. We run both Scenarios together; each Monte Carlo realization is run in parallel as a separate multiprocessing job. 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.
To look at the results we first need to retrieve the unique scenarioId for each child Result Scenario.
resultIDNoHarvest <- subset(resultSummary,
parentID == scenarioId(myScenario))$scenarioId
resultIDHarvest <- subset(resultSummary,
parentID == scenarioId(myScenarioHarvest))$scenarioIdWe then can retrieve tabular output regarding the projected State Class over time (for both scenarios combined) from the stsim_OutputStratumState Datasheet.
outputStratumState <- datasheet(myProject,
scenario = c(resultIDNoHarvest, resultIDHarvest),
name = "stsim_OutputStratumState")Finally, we can get the State Class raster output (here for the Harvest scenario only).
myRastersTimestep5 <- datasheetRaster(myProject,
scenario = resultIDHarvest,
"stsim_OutputSpatialState",
timestep = 5)
# Plot raster for timestep 5 (first realization only)
myRastersTimestep5
plot(myRastersTimestep5[[1]])