In LandSciTech/LSTDConnect: ECCC Landscape Science & Technology Division Connectivity Analysis Tools

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

Setup

devtools::install_github("LandSciTech/LSTDConnect")

We load the package along with raster and samc. We also load the microbenchmark package to benchmark code performance. We will use ggplot for some visualizations.

library(LSTDConnect)
library(raster)
library(samc)
library(microbenchmark)
library(ggplot2)

Then, access the buolt-in dataset, ghm. This raster represent habitat resistance for a small part of Canada (100 cells squares, at a resolution of 1 km).

res <- LSTDConnect::ghm
plot(res)

Running SAMC in paralell

LSTDconnect is optimized to obtain the transient distributions from a SAMC calculation. We will ise the ghm layer as resistance and will use the opposite layer for occupancy occ such as occ = 1-res.

occ <- 101 - res
plot(occ)

Mortality is built from resistance so that when resistance is larger than 80, mortality is high whereas under 80, it is low.

high_mort <- 0.25
low_mort <- 0.01
cutoff <- 80

mort <- res
mort[mort >= cutoff] <- high_mort
mort[mort < cutoff] <- low_mort

We first run LSTDConnect::samc to prepare the data for the computation.

samc_cache <- LSTDConnect::samc(resistance = res, absorption = mort, 
                                directions = 8)

We then use LSTDConnect::distribution, analogous to samc::distribution, to get the transient distribution of individuals for a given time step. We get a 2 elements list, which can be

dists <- LSTDConnect::distribution(samc = samc_cache, occ = occ, time = 1000)
plot(dists$occ)

Comparison with the original samc package

We can compare the results with those given by the samc package.

tr_list <- list(fun = function(x) 1 / mean(x),
                dir = 8, 
                sym = TRUE)

samc_cache_comp <- samc::samc(data = res, absorption = mort, tr_args = tr_list)
dists_comp <- samc::distribution(samc = samc_cache_comp, occ = occ, time = 1000)
dists_mapped <- samc::map(samc_cache_comp, dists_comp)
plot(dists_mapped)

Let's look at the differences. We can see the error is extremely small

diff <- dists_mapped - dists$occ
plot(diff)

hist(values(diff), main = "Error")

Let's now benchmark and compare the speed of these two packages by running each process 100 times

t <- 1000

samc_obj <- samc::samc(data = res,
                       absorption = mort,
                       tr_args = tr_list)  
samc_obj_custom <- LSTDConnect::samc(resistance = res, 
                                     absorption = mort, 
                                     directions = 8)

mbm <- microbenchmark(
  "samc" = {
    short_disp <- samc::distribution(samc = samc_obj,
                                     occ = mort,
                                     time = t)
    gc()
  },
  "LSTDConnect" = {
    short_disp_custom <- LSTDConnect::distribution(samc = samc_obj_custom,
                                                   occ = occ, 
                                                   time = t)
    gc()
  },
  times = 100,
  unit = "s", 
  control = list(order = "block")) 
mbm

boxplot(mbm)

LandSciTech/LSTDConnect documentation built on March 18, 2023, 9:15 p.m.