In andrewmarx/samc: Spatial Absorbing Markov Chains

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

required <- c("viridis")
if (!all(sapply(required, requireNamespace, quietly = TRUE))) {
  knitr::opts_chunk$set(eval = FALSE)
}
options(max.print = 20)

Introduction

This tutorial shows the basics of working with the short-term temporal metrics and how to plot them using the base plotting tools in R.

Setup

# First step is to load the libraries. Not all of these libraries are stricly
# needed; some are used for convenience and visualization for this tutorial.
library("terra")
library("samc")
library("viridisLite")


# "Load" the data. In this case we are using data built into the package.
# In practice, users will likely load raster data using the raster() function
# from the raster package.
res_data <- samc::example_split_corridor$res
abs_data <- samc::example_split_corridor$abs

# To make things easier for plotting later, convert the matrices to rasters
res_data <- rasterize(res_data)
abs_data <- rasterize(abs_data)

# Setup the details for our transition function
rw_model <- list(fun = function(x) 1/mean(x), # Function for calculating transition probabilities
                 dir = 8, # Directions of the transitions. Either 4 or 8.
                 sym = TRUE) # Is the function symmetric?

# Create a samc object using the resistance and absorption data. We use the 
# recipricol of the arithmetic mean for calculating the transition matrix. Note,
# the input data here are matrices, not RasterLayers.
samc_obj <- samc(res_data, abs_data, model = rw_model)

Temporal Analysis

Given a starting location, we're going to calculate the probabilities of where an individual will be at varying time steps. The starting location is going to be cell 1 in the landscape, which is the first non-NA cell going in a left-to-right and top-to-bottom order. The actual analysis and visualization will be shown two ways.

# First, let's specify the time steps we are interested in
time_steps <- c(10, 100, 1000, 10000)

The easiest, most efficient, and safest approach to doing this is using the time step vector inputs added in version 1.1.0.

# First, calculate the metric for the entire vector of time steps
results <- distribution(samc_obj, origin = 1, time = time_steps)

# The result is a list of vectors. Note that the list is named with the time steps
str(results)

# We can take this list, and use map() to convert it to a list of RasterLayers.
results_map <- map(samc_obj, results)
str(results_map, max.level = 1) # max.level is to hide a lot of gory details

# A list of SpatRasters can be turned
results_stack <- rast(results_map)

# Let's look at the names
names(results_stack)

# RasterStacks are convenient for a lot of different things, like processing the
# rasters all at once using spatial methods. But we're just going to plot them
plot(results_stack, xlab = "x", ylab = "y", col = viridis(256))

If you want more control, the following shows how to loop through the results using two different options

# The results of individual time steps can be retrieved by either index or by name 
results[[3]]
results[["1000"]]

# The latter is particularly useful in making reliable for loops because indexed
# for loops can be harder to troubleshoot if something goes wrong. 
for (ts in time_steps) {
  name <- as.character(ts)
  r <- results[[name]]
  r_map <- map(samc_obj, r)
  plot(r_map, main = paste("Individual Location at Time", ts), xlab = "x", ylab = "y", col = viridis(256))
}

# For comparison, here is an indexed for loop that does the same thing
for (i in 1:length(time_steps)) {
  name <- as.character(time_steps[i])
  r <- results[[i]] # or we could replace `i` with `name` on this line
  # ...
}

Previous version required individually calculating the metric for each time step. This is ultimately slower computationally and more error prone for users that are not comfortable with for loops, but is included here for posterity.

for (ts in time_steps) {
  dist <- distribution(samc_obj, origin = 1, time = ts)
  dist_map <- map(samc_obj, dist)
  plot(dist_map, main = paste("Individual Location at Time", ts), xlab = "x", ylab = "y", col = viridis(256))
}