Type:	Package
Title:	Additional Tools for Developing Spatially Explicit Discrete Event Simulation (SpaDES) Models
Description:	Provides GIS and map utilities, plus additional modeling tools for developing cellular automata, dynamic raster models, and agent based models in 'SpaDES'. Included are various methods for spatial spreading, spatial agents, GIS operations, random map generation, and others. See '?SpaDES.tools' for an categorized overview of these additional tools. The suggested package 'NLMR' can be installed from the following repository: (https://PredictiveEcology.r-universe.dev).
URL:	https://spades-tools.predictiveecology.org, https://github.com/PredictiveEcology/SpaDES.tools
Date:	2024-05-14
Version:	2.0.7
Depends:	R (≥ 4.2)
Imports:	backports, checkmate (≥ 1.8.2), data.table (≥ 1.10.4), fpCompare (≥ 0.2.1), graphics, methods, parallel, Rcpp (≥ 0.12.12), reproducible (≥ 2.0.9), stats, terra
Suggests:	animation, bit (≥ 1.1-12), covr, DEoptim (≥ 2.2-4), dqrng, fastmatch, knitr, NLMR (≥ 1.1.1), quickPlot (≥ 1.0.2), raster (≥ 2.5-8), rmarkdown, sf, sp (≥ 1.2-4), testthat (≥ 3.0.0), tools, withr
Encoding:	UTF-8
Language:	en-CA
LinkingTo:	Rcpp
License:	GPL-3
Additional_repositories:	https://predictiveecology.r-universe.dev/
BugReports:	https://github.com/PredictiveEcology/SpaDES.tools/issues
ByteCompile:	yes
Collate:	'RcppExports.R' 'heading.R' 'SELES.R' 'distanceFromEachPoint.R' 'environment.R' 'helpers.R' 'initialize.R' 'mapReduce.R' 'mergeRaster.R' 'movement.R' 'neighbourhood.R' 'numerical-comparisons.R' 'probability.R' 'resample.R' 'rings.R' 'spades-tools-deprecated.R' 'spades-tools-package.R' 'splitRaster.R' 'spread.R' 'spread2.R' 'spread3.R' 'studyArea.R' 'zzz.R'
RoxygenNote:	7.3.1
Config/testthat/edition:	3
NeedsCompilation:	yes
Packaged:	2024-05-15 04:27:45 UTC; achubaty
Author:	Eliot J B McIntire [aut], Alex M Chubaty [aut, cre], Yong Luo [ctb], Ceres Barros [ctb], Steve Cumming [ctb], Jean Marchal [ctb], His Majesty the King in Right of Canada, as represented by the Minister of Natural Resources Canada [cph]
Maintainer:	Alex M Chubaty <achubaty@for-cast.ca>
Repository:	CRAN
Date/Publication:	2024-05-15 09:00:11 UTC

Categorized overview of the SpaDES.tools package

Description

1 Spatial spreading/distances methods

Spatial contagion is a key phenomenon for spatially explicit simulation models. Contagion can be modelled using discrete approaches or continuous approaches. Several functions assist with these:

2 Spatial agent methods

Agents have several methods and functions specific to them:

3 GIS operations

In addition to the vast amount of GIS operations available in R (mostly from contributed packages such as sp, raster, maps, maptools and many others), we provide the following GIS-related functions:

4 Map-reduce - type operations

These functions convert between reduced and mapped representations of the same data. This allows compact representation of, e.g., rasters that have many individual pixels that share identical information.

5 Random Map Generation

It is often useful to build dummy maps with which to build simulation models before all data are available. These dummy maps can later be replaced with actual data maps.

See the NLMR package for tools to generate random landscapes (rasters).

6 SELES-type approach to simulation

These functions are essentially skeletons and are not fully implemented. They are intended to make translations from SELES. You must know how to use SELES for these to be useful:

7 Package options

SpaDES packages use the following options() to configure behaviour:

• spades.lowMemory: If true, some functions will use more memory efficient (but slower) algorithms. Default FALSE.

Author(s)

Maintainer: Alex M Chubaty achubaty@for-cast.ca (ORCID)

Authors:

• Eliot J B McIntire eliot.mcintire@nrcan-rncan.gc.ca (ORCID)

Other contributors:

• Yong Luo Yong.Luo@gov.bc.ca [contributor]

• Ceres Barros ceres.barros@ubc.ca (ORCID) [contributor]

• Steve Cumming Steve.Cumming@sbf.ulaval.ca [contributor]

• Jean Marchal jean.d.marchal@gmail.com [contributor]

• His Majesty the King in Right of Canada, as represented by the Minister of Natural Resources Canada [copyright holder]

See Also

Useful links:

• https://spades-tools.predictiveecology.org

• https://github.com/PredictiveEcology/SpaDES.tools

• Report bugs at https://github.com/PredictiveEcology/SpaDES.tools/issues

Find factors

Description

Internal function (used in gaussMap()). Finds the integer factors of an integer.

Usage

.findFactors(x)

Arguments

Value

A vector of integer factors

Calculate matched point directions

Description

Internal function

Usage

.matchedPointDirection(to, from)

The SpaDES.tools package environment

Description

Environment used internally to store internal package objects and methods.

Usage

.pkgEnv

Format

An object of class environment of length 0.

Alternative point distance (and direction) calculations

Description

These have been written with speed in mind.

Usage

.pointDistance(
  from,
  to,
  angles = NA,
  maxDistance = NA_real_,
  otherFromCols = FALSE
)

Arguments

Fast adjacent function, and Just In Time compiled version

Description

Faster function for determining the cells of the 4, 8 or bishop neighbours of the cells. This is a hybrid function that uses matrix for small numbers of loci (<1e4) and data.table for larger numbers of loci

Usage

adj(
  x = NULL,
  cells,
  directions = 8,
  sort = FALSE,
  pairs = TRUE,
  include = FALSE,
  target = NULL,
  numCol = NULL,
  numCell = NULL,
  match.adjacent = FALSE,
  cutoff.for.data.table = 2000,
  torus = FALSE,
  id = NULL,
  numNeighs = NULL,
  returnDT = FALSE
)

Arguments

Details

Between 4x (large number loci) to 200x (small number loci) speed gains over adjacent in raster package. There is some extra speed gain if NumCol and NumCells are passed rather than a raster. Efficiency gains come from:

1. use data.table internally

   • no need to remove NAs because wrapped or outside points are just removed directly with data.table

   • use data.table to sort and fast select (though not fastest possible)

2. don't make intermediate objects; just put calculation into return statement

The steps used in the algorithm are:

1. Calculate indices of neighbouring cells

2. Remove "to" cells that are

   • ⁠< 1⁠ or ⁠> numCells⁠ (i.e., they are above or below raster), using a single modulo calculation

   • where the modulo of "to" cells is equal to 1 if "from" cells are 0 (wrapped right to left)

   • or where the modulo of the "to" cells is equal to 0 if "from" cells are 1 (wrapped left to right)

Value

Either a matrix (if more than 1 column, i.e., pairs = TRUE, and/or id is provided), a vector (if only one column), or a data.table (if cutoff.for.data.table is less than length(cells) and returnDT is TRUE. To get a consistent output, say a matrix, it would be wise to test the output for its class. The variable output is done to minimize coercion to maintain speed. The columns will be one or more of id, from, to.

Author(s)

Eliot McIntire

See Also

terra::adjacent()

Examples

library(terra)

origDTThreads <- data.table::setDTthreads(2L)
origNcpus <- options(Ncpus = 2L)

a <- rast(ext(0, 1000, 0, 1000), res = 1)
sam <- sample(1:ncell(a), 1e4)
numCol <- ncol(a)
numCell <- ncell(a)
adj.new <- adj(numCol = numCol, numCell = numCell, cells = sam, directions = 8)
adj.new <- adj(numCol = numCol, numCell = numCell, cells = sam, directions = 8,
               include = TRUE)

# clean up
data.table::setDTthreads(origDTThreads)
options(Ncpus = origNcpus)

SELES - Agent Location at initiation

Description

Sets the the location of the initiating agents. NOT YET FULLY IMPLEMENTED.

A SELES-like function to maintain conceptual backwards compatibility with that simulation tool. This is intended to ease transitions from SELES.

You must know how to use SELES for these to be useful.

Usage

agentLocation(map)

Arguments

Value

Object of same class as provided as input. If a ⁠Raster*⁠, then zeros are converted to NA.

Author(s)

Eliot McIntire

Identify pixels in a circle or ring (doughnut) around an object.

Description

Identify the pixels and coordinates that are at a (set of) buffer distance(s) of the objects passed into coords. This is similar to sf::st_buffer but much faster and without the georeferencing information. In other words, it can be used for similar problems, but where speed is important. This code is substantially adapted from PlotRegionHighlighter::createCircle.

Usage

cir(
  landscape,
  coords,
  loci,
  maxRadius = ncol(landscape)/4,
  minRadius = maxRadius,
  allowOverlap = TRUE,
  allowDuplicates = FALSE,
  includeBehavior = "includePixels",
  returnDistances = FALSE,
  angles = NA_real_,
  returnAngles = FALSE,
  returnIndices = TRUE,
  closest = FALSE,
  simplify = TRUE
)

Arguments

Details

This function identifies all the pixels as defined by a donut with inner radius minRadius and outer radius of maxRadius. The includeBehavior defines whether the cells that intersect the radii but whose centres are not inside the donut are included includePixels or not excludePixels in the returned pixels identified. If this is excludePixels, and if a minRadius and maxRadius are equal, this will return no pixels.

Value

A matrix with 4 columns, id, indices, x, y. The x and y indicate the exact coordinates of the indices (i.e., cell number) of the landscape associated with the ring or circle being identified by this function.

See Also

rings() which uses spread internally. cir tends to be faster when there are few starting points, rings tends to be faster when there are many starting points. cir scales with maxRadius^2 and coords. Another difference between the two functions is that rings takes the centre of the pixel as the centre of a circle, whereas cir takes the exact coordinates. See example. For the specific case of creating distance surfaces from specific points, see distanceFromEachPoint(), which is often faster. For the more general GIS buffering, see sf::st_buffer.

Examples

library(data.table)
library(terra)

origDTThreads <- data.table::setDTthreads(2L)
origNcpus <- options(Ncpus = 2L)
set.seed(1462)

# circle centred
ras <- rast(ext(0, 15, 0, 15), res = 1, val = 0)
middleCircle <- cir(ras)
ras[middleCircle[, "indices"]] <- 1
circlePoints <- vect(middleCircle[, c("x", "y")])
if (interactive()) {
  # clearPlot()
  terra::plot(ras)
  terra::plot(circlePoints, add = TRUE)
}

# circles non centred
ras <- randomPolygons(ras, numTypes = 4)
n <- 2
agent <- vect(cbind(x = stats::runif(n, xmin(ras), xmax(ras)),
                    y = stats::runif(n, xmin(ras), xmax(ras))))

cirs <- cir(ras, agent, maxRadius = 15, simplify = TRUE) ## TODO: empty with some seeds! e.g. 1642
cirsSP <- vect(cirs[, c("x", "y")]) ## TODO: error with some seeds! e.g. 1642
cirsRas <- rast(ras)
cirsRas[] <- 0
cirsRas[cirs[, "indices"]] <- 1

if (interactive()) {
  terra::plot(ras)
  terra::plot(cirsRas, add = TRUE, col = c("transparent", "#00000055"))
  terra::plot(agent, add = TRUE)
  terra::plot(cirsSP, add = TRUE)
}

# Example comparing rings and cir
hab <- rast(system.file("extdata", "hab1.tif", package = "SpaDES.tools"))
radius <- 4
n <- 2
coords <- vect(cbind(x = stats::runif(n, xmin(hab), xmax(hab)),
                     y = stats::runif(n, xmin(hab), xmax(hab))))

# cirs
cirs <- cir(hab, coords, maxRadius = rep(radius, length(coords)), simplify = TRUE)

ras1 <- rast(hab)
ras1[] <- 0
ras1[cirs[, "indices"]] <- cirs[, "id"]
if (interactive()) {
  terra::plot(ras1)
}

# rings
loci <- cellFromXY(hab, crds(coords))
cirs2 <- rings(hab, loci, maxRadius = radius, minRadius = radius - 1, returnIndices = TRUE)

ras2 <- rast(hab)
ras2[] <- 0
ras2[cirs2$indices] <- cirs2$id
if (interactive()) {
  terra::plot(c(ras1, ras2))
}

hab <- rast(system.file("extdata", "hab2.tif", package = "SpaDES.tools"))
cirs <- cir(hab, coords, maxRadius = 44, minRadius = 0)
ras1 <- rast(hab)
ras1[] <- 0
cirsOverlap <- data.table::data.table(cirs)[, list(sumIDs = sum(id)), by = indices]
ras1[cirsOverlap$indices] <- cirsOverlap$sumIDs
if (interactive()) {
  terra::plot(ras1)
}

# Provide a specific set of angles
ras <- rast(ext(0, 330, 0, 330), res = 1)
ras[] <- 0
n <- 2
coords <- cbind(x = stats::runif(n, xmin(ras), xmax(ras)),
                y = stats::runif(n, xmin(ras), xmax(ras)))
circ <- cir(ras, coords, angles = seq(0, 2 * pi, length.out = 21),
            maxRadius = 200, minRadius = 0, returnIndices = FALSE,
            allowOverlap = TRUE, returnAngles = TRUE)

# clean up
data.table::setDTthreads(origDTThreads)
options(Ncpus = origNcpus)

This is a very fast version of cir with allowOverlap = TRUE, allowDuplicates = FALSE, returnIndices = TRUE, returnDistances = TRUE, and includeBehavior = "excludePixels". It is used inside spread2, when asymmetry is active. The basic algorithm is to run cir just once, then add to the x,y coordinates of every locus.

Description

This is a very fast version of cir with allowOverlap = TRUE, allowDuplicates = FALSE, returnIndices = TRUE, returnDistances = TRUE, and includeBehavior = "excludePixels". It is used inside spread2, when asymmetry is active. The basic algorithm is to run cir just once, then add to the x,y coordinates of every locus.

Usage

.cirSpecialQuick(landscape, loci, maxRadius, minRadius)

Arguments

Calculate distances and directions between many points and many grid cells

Description

This is a modification of terra::distance() for the case of many points. This version can often be faster for a single point because it does not return a RasterLayer. This is different than terra::distance() because it does not take the minimum distance from the set of points to all cells. Rather this returns the every pair-wise point distance. As a result, this can be used for doing inverse distance weightings, seed rain, cumulative effects of distance-based processes etc. If memory limitation is an issue, maxDistance will keep memory use down, but with the consequences that there will be a maximum distance returned. This function has the potential to use a lot of memory if there are a lot of from and to points.

Internal function.

Usage

directionFromEachPoint(from, to = NULL, landscape)

.pointDirection(from, to)

Arguments

Details

directionFromEachPoint calls .pointDirection, which is not intended to be called directly by the user.

If knowing the which from cell matches with which to cell is important, put a column "id" (e.g., starting cell) in the from matrix.

Value

A sorted matrix on id with same number of rows as to, but with one extra column, angles indicating the angle in radians between from and to. For speed, this angle will be between -pi/2 and 3*pi/2. If the user wants this between say, 0 and 2*pi, then ⁠angles \%\% (2*pi)⁠ will do the trick. See example.

Author(s)

Eliot McIntire

See Also

distanceFromEachPoint(), which will also return directions if angles = TRUE.

Examples

library(terra)

origDTThreads <- data.table::setDTthreads(2L)
origNcpus <- options(Ncpus = 2L)

N <- 2
dirRas <- terra::rast(terra::ext(0,40,0,40), res = 1)
coords <- cbind(x = round(runif(N, xmin(dirRas), xmax(dirRas))) + 0.5,
                y = round(runif(N, xmin(dirRas), xmax(dirRas))) + 0.5,
                id = 1:N)

dirs1 <- directionFromEachPoint(from = coords, landscape = dirRas)

# helper function for converting radians to degrees
deg2 <- function(radian) (radian * 180)/pi
dirs1[, "angles"] <- deg2(dirs1[,"angles"] %% (2*pi))
indices <- cellFromXY(dirRas,dirs1[, c("x", "y")])
minDir <- tapply(dirs1[, "angles"], indices, function(x) min(x)) # minimum angle
dirRas[] <- as.vector(minDir)

if (interactive()) {
  terra::plot(dirRas)
  start <- terra::vect(coords[, c("x", "y"), drop = FALSE])
  terra::plot(start, add = TRUE)
}

# clean up
data.table::setDTthreads(origDTThreads)
options(Ncpus = origNcpus)

Calculate distances and directions between many points and many grid cells

Description

This is a modification of terra::distance() for the case of many points. This version can often be faster for a single point because it does not return a RasterLayer. This is different than terra::distance() because it does not take the minimum distance from the set of points to all cells. Rather this returns the every pair-wise point distance. As a result, this can be used for doing inverse distance weightings, seed rain, cumulative effects of distance-based processes etc. If memory limitation is an issue, maxDistance will keep memory use down, but with the consequences that there will be a maximum distance returned. This function has the potential to use a lot of memory if there are a lot of from and to points.

Usage

distanceFromEachPoint(
  from,
  to = NULL,
  landscape,
  angles = NA_real_,
  maxDistance = NA_real_,
  cumulativeFn = NULL,
  distFn = function(dist) 1/(1 + dist),
  cl,
  ...
)

Arguments

Details

This function is cluster aware if the raster package is available. If there is a cluster running, it will use it. To start a cluster use raster::beginCluster(), with N being the number of cores to use. See examples in SpaDES.core::experiment.

If the user requires an id (indicating the from cell for each to cell) to be returned with the function, the user must add an identifier to the from matrix, such as "id". Otherwise, the function will only return the coordinates and distances.

distanceFromEachPoint calls .pointDistance, which is not intended to be called directly by the user.

This function has the potential to return a very large object, as it is doing pairwise distances (and optionally directions) between from and to. If there are memory limitations because there are many from and many to points, then cumulativeFn and distFn can be used. These two functions together will be used iteratively through the from points. The distFn should be a transformation of distances to be used by the cumulativeFn function. For example, if distFn is 1 / (1+x), the default, and cumulativeFn is +, then it will do a sum of inverse distance weights. See examples.

Value

A sorted matrix on id with same number of rows as to, but with one extra column, "dists", indicating the distance between from and to.

See Also

rings(), cir(), terra::distance(), which can all be made to do the same thing, under specific combinations of arguments. But each has different primary use cases. Each is also faster under different conditions. For instance, if maxDistance is relatively small compared to the number of cells in the landscape, then cir() will likely be faster. If a minimum distance from all cells in the landscape to any cell in from, then distanceFromPoints will be fastest. This function scales best when there are many to points or all cells are used to = NULL (which is default).

Examples

library(terra)

origDTThreads <- data.table::setDTthreads(2L)
origNcpus <- options(Ncpus = 2L)

n <- 2
distRas <- rast(ext(0, 40, 0, 40), res = 1)
coords <- cbind(x = round(runif(n, xmin(distRas), xmax(distRas))) + 0.5,
                y = round(runif(n, xmin(distRas), xmax(distRas))) + 0.5)

# inverse distance weights
dists1 <- distanceFromEachPoint(coords, landscape = distRas)
indices <- cellFromXY(distRas, dists1[, c("x", "y")])
invDist <- tapply(dists1[, "dists"], indices, function(x) sum(1 / (1 + x))) # idw function
distRas[] <- as.vector(invDist)
if (interactive()) {
  # clearPlot()
  terra::plot(distRas)
}

# With iterative summing via cumulativeFn to keep memory use low, with same result
dists1 <- distanceFromEachPoint(coords[, c("x", "y"), drop = FALSE],
                                landscape = distRas, cumulativeFn = `+`)
idwRaster <- rast(distRas)
idwRaster[] <- dists1[, "dists"]
if (interactive()) terra::plot(idwRaster)

all(idwRaster[] == distRas[]) # TRUE

# A more complex example of cumulative inverse distance sums, weighted by the value
#  of the origin cell
ras <- rast(ext(0, 34, 0, 34), res = 1, val = 0)
rp <- randomPolygons(ras, numTypes = 10) ^ 2
n <- 15
cells <- sample(ncell(ras), n)
coords <- xyFromCell(ras, cells)
distFn <- function(landscape, fromCell, dist) landscape[fromCell] / (1 + dist)

#beginCluster(3) # can do parallel
dists1 <- distanceFromEachPoint(coords[, c("x", "y"), drop = FALSE],
                                landscape = rp, distFn = distFn, cumulativeFn = `+`)
#endCluster() # if beginCluster was run

idwRaster <- rast(ras)
idwRaster[] <- dists1[, "dists"]
if (interactive()) {
  # clearPlot()
  terra::plot(rp)
  sp1 <- vect(coords)
  terra::plot(sp1, add = TRUE)
  terra::plot(idwRaster)
  terra::plot(sp1, add = TRUE)
}

# On linux; can use numeric passed to cl; will use mclapply with mc.cores = cl
if (identical(Sys.info()["sysname"], "Linux")) {
  dists1 <- distanceFromEachPoint(coords[, c("x", "y"), drop = FALSE],
                                  landscape = rp, distFn = distFn,
                                  cumulativeFn = `+`, cl = 2)
}

# clean up
data.table::setDTthreads(origDTThreads)
options(Ncpus = origNcpus)

Rcpp duplicated on integers using Rcpp Sugar

Description

.duplicatedInt does same as duplicated in R, but only on integers, and faster. It uses Rcpp sugar

Usage

duplicatedInt(x)

Arguments

Value

A logical vector, as per duplicated

Vectorized wrapped normal density function

Description

This is a modified version of dwrpnorm in the CircStats package to allow for multiple angles at once (i.e., vectorized on theta and mu).

Usage

dwrpnorm2(theta, mu, rho, sd = 1, acc = 1e-05, tol = acc)

Arguments

Author(s)

Eliot McIntire

Examples

# Values for which to evaluate density
theta <- c(1:500) * 2 * pi / 500
# Compute wrapped normal density function
density <- c(1:500)
for(i in 1:500) {
  density[i] <- dwrpnorm2(theta[i], pi, .75)
}

if (interactive()) {
  plot(theta, density)
}

# Approximate area under density curve
sum(density * 2 * pi / 500)

fastCrop is deprecated.

Description

fastCrop is deprecated.

Usage

fastCrop(x, y, ...)

Arguments

See Also

velox::VeloxRaster_crop

Produce a raster of a random Gaussian process.

Description

Defunct.

Usage

gaussMap(
  x,
  scale = 10,
  var = 1,
  speedup = 1,
  method = "RMexp",
  alpha = 1,
  inMemory = FALSE,
  ...
)

Arguments

Value

A raster map with same extent as x, with a Gaussian random pattern.

Heading between spatial points.

Description

Determines the heading between spatial points.

Usage

heading(from, to)

Arguments

Value

The heading between the points, in degrees.

Author(s)

Eliot McIntire

Examples

library(terra)

origDTThreads <- data.table::setDTthreads(2L)
origNcpus <- options(Ncpus = 2L)

N <- 10L                       # number of agents
x1 <- stats::runif(N, -50, 50) # previous X location
y1 <- stats::runif(N, -50, 50) # previous Y location
x0 <- stats::rnorm(N, x1, 5)   # current X location
y0 <- stats::rnorm(N, y1, 5)   # current Y location

# using SpatVector
prev <- terra::vect(cbind(x = x1, y = y1))
curr <- terra::vect(cbind(x = x0, y = y0))
heading(prev, curr)

# using matrix
prev <- matrix(c(x1, y1), ncol = 2, dimnames = list(NULL, c("x","y")))
curr <- matrix(c(x0, y0), ncol = 2, dimnames = list(NULL, c("x","y")))
heading(prev, curr)

#using both
prev <- terra::vect(cbind(x = x1, y = y1))
curr <- matrix(c(x0, y0), ncol = 2, dimnames = list(NULL, c("x","y")))
heading(prev, curr)

prev <- matrix(c(x1, y1), ncol = 2, dimnames = list(NULL, c("x","y")))
curr <- terra::vect(cbind(x = x0, y = y0))
heading(prev, curr)

# clean up
data.table::setDTthreads(origDTThreads)
options(Ncpus = origNcpus)

Test whether a number lies within range ⁠[a,b]⁠

Description

Default values of ⁠a=0; b=1⁠ allow for quick test if x is a probability.

Usage

inRange(x, a = 0, b = 1)

Arguments

Value

Logical vectors. NA values in x are retained.

Author(s)

Alex Chubaty

Examples

set.seed(100)
x <- stats::rnorm(4) # -0.50219235  0.13153117 -0.07891709  0.88678481
inRange(x, 0, 1)

SELES - Initiate agents

Description

Sets the the number of agents to initiate. THIS IS NOT FULLY IMPLEMENTED.

A SELES-like function to maintain conceptual backwards compatibility with that simulation tool. This is intended to ease transitions from SELES.

You must know how to use SELES for these to be useful.

Usage

initiateAgents(map, numAgents, probInit, asSpatialPoints = TRUE, indices)

Arguments

Value

A SpatialPointsDataFrame, with each row representing an individual agent

Author(s)

Eliot McIntire

Examples

if (require("sf", quietly = TRUE)) {
  library(data.table)
  library(terra)

  origDTThreads <- data.table::setDTthreads(2L)
  origNcpus <- options(Ncpus = 2L)

  map <- rast(system.file("extdata", "map.tif", package = "SpaDES.tools"))
  names(map) <- "layer"
  pr <- probInit(map, p = (map[] / terra::minmax(map)[2])^2)
  agents <- initiateAgents(map, 100, pr, asSpatialPoints = "sf")
  if (interactive()) {
    terra::plot(map)
    terra::plot(agents, add = TRUE)
  }
  # Test that they are indeed selecting according to probabilities in pr
  dt1 <- data.table(table(round(extract(map, agents), 0)[, "layer"]))
  setnames(dt1, old = "N", new = "count")
  dt2 <- data.table(table(round(map[], 0)))
  setnames(dt2, old = "N", new = "available")
  dt <- dt1[dt2, on = "V1"]  # join the counts and available data.tables
  setnames(dt, old = "V1", new = "mapValue")
  dt[, selection := count / available]
  dt[is.na(selection), selection := 0]
  if (interactive()) {
    with(dt, plot(mapValue, selection))
  }
  #'
  # Note, can also produce a Raster representing agents,
  # then the number of points produced can't be more than
  # the number of pixels:
  agentsRas <- initiateAgents(map, 30, pr, asSpatialPoints = FALSE)
  if (interactive()) {
    terra::plot(agentsRas)
  }
  #'
  # clean up
  data.table::setDTthreads(origDTThreads)
  options(Ncpus = origNcpus)
}

Split and re-merge RasterLayer(s)

Description

splitRaster divides up a raster into an arbitrary number of pieces (tiles). Split rasters can be recombined using do.call(merge, y) or mergeRaster(y), where y <- splitRaster(x).

Usage

mergeRaster(x, fun = NULL)

## S4 method for signature 'list'
mergeRaster(x, fun = NULL)

splitRaster(
  r,
  nx = 1,
  ny = 1,
  buffer = c(0, 0),
  path = NA,
  cl,
  rType = "FLT4S",
  fExt = ".tif"
)

Arguments

Details

mergeRaster differs from merge in how overlapping tile regions are handled: merge retains the values of the first raster in the list. This has the consequence of retaining the values from the buffered region in the first tile in place of the values from the neighbouring tile. On the other hand, mergeRaster retains the values of the tile region, over the values in any buffered regions. This is useful for reducing edge effects when performing raster operations involving contagious processes.

This function is parallel-aware using the same mechanism as used in raster: NOTE: This may not work as expected as we transition away from raster. Specifically, if you start a cluster using raster::beginCluster(), then this function will automatically use that cluster. It is always a good idea to stop the cluster when finished, using raster::endCluster().

Value

mergeRaster returns a RasterLayer object.

splitRaster returns a list (length nx*ny) of cropped raster tiles.

Author(s)

Yong Luo, Alex Chubaty, Tati Micheletti & Ian Eddy

Alex Chubaty and Yong Luo

See Also

terra::merge(), terra::mosaic()

do.call(), terra::merge().

Examples

library(terra)

origDTThreads <- data.table::setDTthreads(2L)
origNcpus <- options(Ncpus = 2L)
set.seed(1462)

## an example with dimensions: nrow = 77,  ncol = 101, nlayers = 3
b <- rast(system.file("ex/logo.tif", package = "terra"))
r <- b[[1]] # use first layer only
nx <- 3
ny <- 4

tmpdir <- dir.create(file.path(tempdir(), "splitRaster-example"), showWarnings = FALSE)

y0 <- splitRaster(r, nx, ny, path = file.path(tmpdir, "y0")) # no buffer

## buffer: 10 pixels along both axes
y1 <- splitRaster(r, nx, ny, c(10, 10), path = file.path(tmpdir, "y1"))

## buffer: half the width and length of each tile
y2 <- splitRaster(r, nx, ny, c(0.5, 0.5), path = file.path(tmpdir, "y2"))

## the original raster:
if (interactive()) plot(r) # may require a call to `dev()` if using RStudio

## the split raster:
layout(mat = matrix(seq_len(nx * ny), ncol = nx, nrow = ny))
plotOrder <- unlist(lapply(split(1:12, rep(1:nx, each = ny)), rev))

if (interactive()) {
  invisible(lapply(y0[plotOrder], terra::plot))
}

## parallel splitting
if (requireNamespace("raster", quietly = TRUE) &&
    requireNamespace("parallel")) {
  if (interactive()) {
    n <- pmin(parallel::detectCores(), 4) # use up to 4 cores
    raster::beginCluster(n, type = "PSOCK")
    y3 <- splitRaster(r, nx, ny, c(0.7, 0.7), path = file.path(tmpdir, "y3"))
    raster::endCluster()
    if (interactive()) {
      invisible(lapply(y3[plotOrder], terra::plot))
    }
  }
}

## can be recombined using `terra::merge`
m0 <- do.call(merge, y0)
all.equal(m0, r) ## TRUE

m1 <- do.call(merge, y1)
all.equal(m1, r) ## TRUE

m2 <- do.call(merge, y2)
all.equal(m2, r) ## TRUE

## or recombine using mergeRaster
n0 <- mergeRaster(y0)
all.equal(n0, r) ## TRUE

n1 <- mergeRaster(y1)
all.equal(n1, r) ## TRUE

n2 <- mergeRaster(y2)
all.equal(n2, r) ## TRUE

# clean up
data.table::setDTthreads(origDTThreads)
options(Ncpus = origNcpus)
unlink(tmpdir, recursive = TRUE)

Return the (approximate) middle pixel on a raster

Description

Return the (approximate) middle pixel on a raster

Usage

middlePixel(ras)

Arguments

Value

This calculation is slightly different depending on whether the nrow(ras) and ncol(ras) are even or odd. It will return the exact middle pixel if these are odd, and the pixel just left and/or above the middle pixel if either dimension is even, respectively.

Move

Description

Wrapper for selecting different animal movement methods.

This version uses just turn angles and step lengths to define the correlated random walk.

Usage

move(hypothesis = "crw", ...)

crw(
  agent,
  extent,
  stepLength,
  stddev,
  lonlat = FALSE,
  torus = FALSE,
  returnMatrix = FALSE
)

Arguments

Details

This simple version of a correlated random walk is largely the version that was presented in Turchin 1998, but it was also used with bias modifications in McIntire, Schultz, Crone 2007.

Value

A SpatVector points object with updated spatial position defined by a single occurrence of step length(s) and turn angle(s).

Author(s)

Eliot McIntire

References

Turchin, P. 1998. Quantitative analysis of movement: measuring and modeling population redistribution in animals and plants. Sinauer Associates, Sunderland, MA.

McIntire, E. J. B., C. B. Schultz, and E. E. Crone. 2007. Designing a network for butterfly habitat restoration: where individuals, populations and landscapes interact. Journal of Applied Ecology 44:725-736.

See Also

terra::distance()

wrap()

Examples

origDTThreads <- data.table::setDTthreads(2L)
origNcpus <- options(Ncpus = 2L)

# using just matrix
N <- 10
xrange <- yrange <- c(-50, 50)
starts <- cbind(x = stats::runif(N, xrange[1], xrange[2]),
                y = stats::runif(N, yrange[1], yrange[2]))
moved <- crw(starts, stepLength = 5, stddev = 10)
plot(starts, col = rainbow(10), pch = 19)
points(moved, col = rainbow(10))

# as SpatVector
agent <- terra::vect(starts)
moved <- crw(agent, stepLength = 5, stddev = 10)
movedAgain <- crw(moved, stepLength = 5, stddev = 10)
terra::plot(agent)
terra::plot(moved, add = TRUE, col = "red")
terra::plot(movedAgain, add = TRUE, col = "green")

# 1000x faster!! -- returnMatrix = TRUE
agentOrig <- agent
reps <- 1e2
system.time({
  for (i in 1:reps) agent <- crw(agent, stepLength = 5, stddev = 10, returnMatrix = TRUE)
})
agent <- agentOrig
system.time({
  for (i in 1:reps) agent <- crw(agent, stepLength = 5, stddev = 10)
})

# as sp
if (requireNamespace("sp")) {
  agent <- sp::SpatialPoints(starts)
  spdf <- crw(agent, stepLength = 5, stddev = 10)
  spdfNew <- crw(spdf, stepLength = 5, stddev = 10)
  terra::plot(spdf, pch = 19)
  terra::points(spdfNew, col = "blue", pch = 19)
}

# clean up
data.table::setDTthreads(origDTThreads)
options(Ncpus = origNcpus)

Produce a neutral landscape using a midpoint displacement algorithm

Description

This is a wrapper for the nlm_mpd function in the NLMR package. The main addition is that it makes sure that the output raster conforms in extent with the input raster x, since nlm_mpd can output a smaller raster.

Usage

neutralLandscapeMap(
  x,
  pad = 10L,
  type = c("nlm_mpd", "nlm_gaussianfield", "nlm_distancegradient", "nlm_edgegradient",
    "nlm_fbm", "nlm_mosaicfield", "nlm_mosaicgibbs", "nlm_mosaictess", "nlm_neigh",
    "nlm_percolation", "nlm_planargradient", "nlm_random",
    "nlm_randomrectangularcluster"),
  ...
)

Arguments

Value

A RasterLayer/SpatRaster with same extent as x, with randomly generated values.

See Also

nlm_mpd

Examples

  if (requireNamespace("NLMR", quietly = TRUE) &&
      requireNamespace("raster", quietly = TRUE)) {
    library(terra)
    nx <- ny <- 100L
    r <- rast(nrows = ny, ncols = nx,
              xmin = -nx/2, xmax = nx/2,
              ymin = -ny/2, ymax = ny/2)
    ## or with raster package:
    # r <- raster::raster(nrows = ny, ncols = nx,
    #                     xmn = -nx/2, xmx = nx/2,
    #                     ymn = -ny/2, ymx = ny/2)
    map1 <- neutralLandscapeMap(r,
                                type = "nlm_mpd",
                                roughness = 0.65,
                                rand_dev = 200,
                                rescale = FALSE,
                                verbose = FALSE)
    if (interactive()) plot(map1)
  }

SELES - Number of Agents to initiate

Description

Sets the the number of agents to initiate. THIS IS NOT YET FULLY IMPLEMENTED.

A SELES-like function to maintain conceptual backwards compatibility with that simulation tool. This is intended to ease transitions from SELES.

You must know how to use SELES for these to be useful.

Usage

numAgents(N, probInit)

Arguments

Value

A numeric, indicating number of agents to start

Author(s)

Eliot McIntire

Patch size

Description

Patch size

Usage

patchSize(patches)

Arguments

SELES - Probability of Initiation

Description

Describes the probability of initiation of agents or events. THIS IS NOT FULLY IMPLEMENTED.

A SELES-like function to maintain conceptual backwards compatibility with that simulation tool. This is intended to ease transitions from SELES.

You must know how to use SELES for these to be useful.

Usage

probInit(map, p = NULL, absolute = NULL)

Arguments

Value

A RasterLayer with probabilities of initialization. There are several combinations of inputs possible and they each result in different behaviours.

If p is numeric or Raster and between 0 and 1, it is treated as an absolute probability, and a map will be produced with the p value(s) everywhere.

If p is numeric or Raster and not between 0 and 1, it is treated as a relative probability, and a map will be produced with p/max(p) value(s) everywhere.

If absolute is provided, it will override the previous statements, unless absolute = TRUE and p is not between 0 and 1 (i.e., is not a probability).

Author(s)

Eliot McIntire

Produce a SpatRaster of random polygons

Description

These are built with the spread() function internally.

Produces a SpatVector polygons object with 1 feature that will have approximately an area equal to area (expecting area in hectares), #' and a centre at approximately x.

Usage

randomPolygons(
  ras = rast(ext(0, 15, 0, 15), res = 1, vals = 0),
  numTypes = 2,
  ...
)

randomPolygon(x, hectares, area)

## Default S3 method:
randomPolygon(x, hectares, area)

Arguments

Value

A map of extent ext with random polygons.

A SpatVector polygons object, with approximately the area request, centred approximately at the coordinates requested, in the projection of x.

See Also

spread(), randomPolygons()

Examples

origDTThreads <- data.table::setDTthreads(2L)
origNcpus <- options(Ncpus = 2L)

set.seed(1234)
Ras <- randomPolygons(numTypes = 5)
if (interactive() ) {
  terra::plot(Ras, col = c("yellow", "dark green", "blue", "dark red"))
}

# more complex patterning, with a range of patch sizes
r <- terra::rast(terra::ext(0, 50, 0, 50), resolution = 1, vals = 0)
a <- randomPolygons(numTypes = 400, r)
a[a < 320] <- 0
a[a >= 320] <- 1
clumped <- terra::patches(a)
if (interactive()) {
  terra::plot(a)
}

# clean up
data.table::setDTthreads(origDTThreads)
options(Ncpus = origNcpus)

a1 <- terra::vect(cbind(-110, 59), crs = "epsg:4326")
a2 <- randomPolygon(a1, area = 1e7)

if (interactive()) {
  terra::plot(a1)
  terra::points(a2, pch = 19)
}

if (require("sf", quietly = TRUE)) {
  b1 <- list(cbind(
    x = c(-122.98, -116.1, -99.2, -106, -122.98),
    y = c(59.9, 65.73, 63.58, 54.79, 59.9)
  )) |>
    sf::st_polygon() |>
    sf::st_sfc() |>
    sf::st_sf(geometry = _, ID = 1L, shinyLabel = "zone2", crs = "epsg:4326")
  b2 <- randomPolygon(b1, area = 1e10)

  if (interactive()) {
    plot(sf::st_geometry(b1))
    plot(sf::st_geometry(b2), add = TRUE)
  }
}

Create default study areas for use with SpaDES modules

Description

Create default study areas for use with SpaDES modules

Usage

randomStudyArea(center = NULL, size = 10000, seed = NULL)

Arguments

Value

SpatVector

Examples

a <- randomStudyArea(seed = 123)
if (interactive()) {
  terra::plot(a)
}

Convert reduced representation to full raster

Description

Convert reduced representation to full raster

Usage

rasterizeReduced(
  reduced,
  fullRaster,
  newRasterCols,
  mapcode = names(fullRaster),
  ...
)

Arguments

Value

A RasterLayer/SpatRaster or list of RasterLayer/SpatRaster of with same dimensions as fullRaster representing newRasterCols spatially, according to the join between the mapcode contained within reduced and fullRaster

Author(s)

Eliot McIntire

See Also

terra::rast()

Examples

library(data.table)
library(terra)

origDTThreads <- data.table::setDTthreads(2L)
origNcpus <- options(Ncpus = 2L)

ras <- rast(ext(0, 15, 0, 15), res = 1)
fullRas <- randomPolygons(ras, numTypes = 2)
names(fullRas) <- "mapcodeAll"
uniqueComms <- unique(fullRas)
reducedDT <- data.table(uniqueComms,
                        communities = sample(1:1000, length(uniqueComms)),
                        biomass = rnbinom(length(uniqueComms), mu = 4000, 0.4))
biomass <- rasterizeReduced(reducedDT, fullRas, "biomass")

# The default key is the layer name of the fullRas, so rekey incase of miskey
setkey(reducedDT, biomass)

communities <- rasterizeReduced(reducedDT, fullRas, "communities")
coltab(communities) <- c("blue", "orange", "red")
if (interactive()) {
  terra::plot(c(biomass, communities, fullRas))
}

## with a factor SpatRaster, the mapcode should correspond to the
## active category (not the ids)
cls <- data.frame(id = sort(unique(as.vector(fullRas[]))))
cls$Bclass <- LETTERS[cls$id]
levels(fullRas) <- cls
is.factor(fullRas)

clsDT <- as.data.table(cls)
reducedDT <- reducedDT[clsDT, on = "mapcodeAll==id"]
reducedDT[, mapcodeAll := Bclass]

biomass2 <- rasterizeReduced(reducedDT, fullRas, "biomass")

# clean up
data.table::setDTthreads(origDTThreads)
options(Ncpus = origNcpus)

Internal helpers

Description

Not for users. A function to setnames and rbindlist that is used in spread2.

Usage

rbindlistDtDtpot(
  dt,
  dtPotential,
  returnFrom,
  needDistance,
  dtPotentialColNames
)

reorderColsWDistance(needDistance, dtPotential, dtPotentialColNames)

angleQuality(from, to, landscape, actualAsymmetryAngle)

asymmetryAdjust(angleQualities, quantity, actualAsymmetry)

Arguments

Adapted directly from the base::sample() help file.

Description

Adapted directly from the base::sample() help file.

resampleZeroProof is a version that works even if sum of all probabilities passed to sample.int is zero. This causes an error in sample.int.

Usage

resample(x, ...)

resampleZeroProof(spreadProbHas0, x, n, prob)

Arguments

Note

Intended for internal use only.

Identifies all cells within a ring around the focal cells

Description

Identifies the cell numbers of all cells within a ring defined by minimum and maximum distances from focal cells. Uses spread() under the hood, with specific values set. Under many situations, this may be faster than using sf::st_buffer twice (once for smaller ring and once for larger ring, then removing the smaller ring cells).

Usage

rings(
  landscape,
  loci = NA_real_,
  id = FALSE,
  minRadius = 2,
  maxRadius = 5,
  allowOverlap = FALSE,
  returnIndices = FALSE,
  returnDistances = TRUE,
  ...
)

Arguments

Value

This will return a data.table with columns as described in spread when returnIndices = TRUE.

Author(s)

Eliot McIntire

See Also

cir() which uses a different algorithm. cir tends to be faster when there are few starting points, rings tends to be faster when there are many starting points. Another difference between the two functions is that rings takes the centre of the pixel as the centre of a circle, whereas cir takes the exact coordinates. See example.

sf::st_buffer

Examples

library(terra)

origDTThreads <- data.table::setDTthreads(2L)
origNcpus <- options(Ncpus = 2L)
set.seed(1462)

# Make random forest cover map
emptyRas <- terra::rast(terra::ext(0, 1e2, 0, 1e2), res = 1)

# start from two cells near middle
loci <- (ncell(emptyRas) / 2 - ncol(emptyRas)) / 2 + c(-3, 3)

# No overlap is default, occurs randomly
emptyRas[] <- 0
rngs <- rings(emptyRas, loci = loci, minRadius = 7, maxRadius = 9, returnIndices = TRUE)
emptyRas[rngs$indices] <- rngs$id
if (interactive()) {
  terra::plot(emptyRas)
}

# Variable ring widths, including centre cell for smaller one
emptyRas[] <- 0
rngs <- rings(emptyRas, loci = loci, minRadius = c(0, 7), maxRadius = c(8, 18),
              returnIndices = TRUE)
emptyRas[rngs$indices] <- rngs$id
if (interactive()) {
  terra::plot(emptyRas)
}

# clean up
data.table::setDTthreads(origDTThreads)
options(Ncpus = origNcpus)

Rcpp Sugar version of runif

Description

Slightly faster than runif, and used a lot

Usage

runifC(N)

Arguments

Value

A vector of uniform random numbers as per runif

Resample

Description

A version of sample that doesn't have awkward behaviour when length(x) == 1. Vectorized version of sample using Vectorize

Usage

sampleV(x, size, replace = FALSE, prob = NULL)

Arguments

Details

Intended for internal use only. size is vectorized.

Value

A random permutation, as in sample, but with size vectorized.

Initiate a specific number of agents in a map of patches

Description

Instantiate a specific number of agents per patch. The user can either supply a table of how many to initiate in each patch, linked by a column in that table called pops.

Usage

specificNumPerPatch(patches, numPerPatchTable = NULL, numPerPatchMap = NULL)

Arguments

Value

A raster with 0s and 1s, where the 1s indicate starting locations of agents following the numbers above.

Examples

library(data.table)

origDTThreads <- data.table::setDTthreads(2L)
origNcpus <- options(Ncpus = 2L)

set.seed(1234)
Ntypes <- 4
ras <- randomPolygons(numTypes = Ntypes)
if (interactive()) {
  terra::plot(ras)
}

# Use numPerPatchTable
patchDT <- data.table(pops = 1:Ntypes, num.in.pop = c(1, 3, 5, 7))
rasAgents <- specificNumPerPatch(ras, patchDT)
rasAgents[is.na(rasAgents)] <- 0

if (require(testthat))
  expect_true(all(unname(table(ras[rasAgents])) == patchDT$num.in.pop))

# Use numPerPatchMap
rasPatches <- ras
for (i in 1:Ntypes) {
  rasPatches[rasPatches==i] <- patchDT$num.in.pop[i]
}
if (interactive()) {
  terra::plot(c(ras, rasPatches))
}
rasAgents <- specificNumPerPatch(ras, numPerPatchMap = rasPatches)
rasAgents[is.na(rasAgents)] <- 0
if (interactive()) {
  terra::plot(rasAgents)
}

# clean up
data.table::setDTthreads(origDTThreads)
options(Ncpus = origNcpus)

Identify outward radiating spokes from initial points

Description

This is a generalized version of a notion of a viewshed. The main difference is that there can be many "viewpoints".

Usage

spokes(
  landscape,
  coords,
  loci,
  maxRadius = ncol(landscape)/4,
  minRadius = maxRadius,
  allowOverlap = TRUE,
  stopRule = NULL,
  includeBehavior = "includePixels",
  returnDistances = FALSE,
  angles = NA_real_,
  nAngles = NA_real_,
  returnAngles = FALSE,
  returnIndices = TRUE,
  ...
)

Arguments

Value

A matrix containing columns id (representing the row numbers of coords), angles (from coords to each point along the spokes), x and y coordinates of each point along the spokes, the corresponding indices on the landscape Raster, dists (the distances between each coords and each point along the spokes), and stop, indicating if it was a point that caused a spoke to stop going outwards due to stopRule.

Author(s)

Eliot McIntire

Examples

library(data.table)
library(terra)

origDTThreads <- data.table::setDTthreads(2L)
origNcpus <- options(Ncpus = 2L)
set.seed(1234)

ras <- terra::rast(terra::ext(0, 10, 0, 10), res = 1, val = 0)
rp <- randomPolygons(ras, numTypes = 10)

if (interactive())
  terra::plot(rp)

angles <- seq(0, pi * 2, length.out = 17)
angles <- angles[-length(angles)]
n <- 2
loci <- sample(terra::ncell(rp), n)
coords <- terra::vect(terra::xyFromCell(rp, loci))
stopRule <- function(landscape) landscape < 3
d2 <- spokes(rp, coords = coords, stopRule = stopRule,
             minRadius = 0, maxRadius = 50,
             returnAngles = TRUE, returnDistances = TRUE,
             allowOverlap = TRUE, angles = angles, returnIndices = TRUE)

# Assign values to the "patches" that were in the viewshed of a ray
rasB <- terra::rast(ras)
rasB[] <- 0
rasB[d2[d2[, "stop"] == 1, "indices"]] <- 1

if (interactive()) {
  rasB[rasB == 0] <- NA
  terra::plot(rasB, add = TRUE, col = "red", legend = FALSE)
}

if (NROW(d2) > 0) {
  sp1 <- terra::vect(d2[, c("x", "y")])
  if (interactive())
    terra::plot(sp1, add = TRUE, pch = 19)
}
if (interactive())
  terra::plot(coords, add = TRUE, pch = 19, col = "blue")

# clean up
data.table::setDTthreads(origDTThreads)
options(Ncpus = origNcpus)

Simulate a spread process on a landscape.

Description

This can be used to simulate fires, seed dispersal, calculation of iterative, concentric landscape values (symmetric or asymmetric) and many other things. Essentially, it starts from a collection of cells (loci) and spreads to neighbours, according to the directions and spreadProb arguments. This can become quite general, if spreadProb is 1 as it will expand from every loci until all cells in the landscape have been covered. With id set to TRUE, the resulting map will be classified by the index of the cell where that event propagated from. This can be used to examine things like fire size distributions. NOTE: See also spread2(), which is more robust and can be used to build custom functions. However, under some conditions, this spread function is faster. The two functions can accomplish many of the same things, and key differences are internal.

Usage

spread(
  landscape,
  loci = NA_real_,
  spreadProb = 0.23,
  persistence = 0,
  mask = NA,
  maxSize = 100000000L,
  directions = 8L,
  iterations = 1000000L,
  lowMemory = NULL,
  returnIndices = FALSE,
  returnDistances = FALSE,
  mapID = NULL,
  id = FALSE,
  plot.it = FALSE,
  spreadProbLater = NA_real_,
  spreadState = NA,
  circle = FALSE,
  circleMaxRadius = NA_real_,
  stopRule = NA,
  stopRuleBehavior = "includeRing",
  allowOverlap = FALSE,
  asymmetry = NA_real_,
  asymmetryAngle = NA_real_,
  quick = FALSE,
  neighProbs = NULL,
  exactSizes = FALSE,
  relativeSpreadProb = FALSE,
  ...
)

Arguments

Details

For large rasters, a combination of lowMemory = TRUE and returnIndices = TRUE or returnIndices = 2 will be fastest and use the least amount of memory. 2022-07-25: lowMemory = TRUE is deprecated due to removal of package ffbase from CRAN.

This function can be interrupted before all active cells are exhausted if the iterations value is reached before there are no more active cells to spread into. If this is desired, returnIndices should be TRUE and the output of this call can be passed subsequently as an input to this same function. This is intended to be used for situations where external events happen during a spread event, or where one or more arguments to the spread function change before a spread event is completed. For example, if it is desired that the spreadProb change before a spread event is completed because, for example, a fire is spreading, and a new set of conditions arise due to a change in weather.

asymmetry is currently used to modify the spreadProb in the following way. First for each active cell, spreadProb is converted into a length 2 numeric of Low and High spread probabilities for that cell: ⁠spreadProbsLH <- (spreadProb*2) // (asymmetry+1)*c(1,asymmetry)⁠, whose ratio is equal to asymmetry. Then, using asymmetryAngle, the angle between the initial starting point of the event and all potential cells is found. These are converted into a proportion of the angle from -asymmetryAngle to asymmetryAngle using: angleQuality <- (cos(angles - rad2(asymmetryAngle))+1)/2 where rad2 <- function (degree) (degree * pi)/180

These are then converted to multiple spreadProbs by spreadProbs <- lowSpreadProb + (angleQuality * diff(spreadProbsLH)) To maintain an expected spreadProb that is the same as the asymmetric spreadProbs, these are then rescaled so that the mean of the asymmetric spreadProbs is always equal to spreadProb at every iteration: spreadProbs <- spreadProbs - diff(c(spreadProb, mean(spreadProbs)))

Value

Either a RasterLayer indicating the spread of the process in the landscape or a data.table if returnIndices is TRUE. If a RasterLayer, then it represents every cell in which a successful spread event occurred. For the case of, say, a fire this would represent every cell that burned. If allowOverlap is TRUE, This RasterLayer will represent the sum of the individual event ids (which are numerics seq_along(loci). This will generally be of minimal use because it won't be possible to distinguish if event 2 overlapped with event 5 or if it was just event 7.

If returnIndices is TRUE, then this function returns a data.table with columns:

This will generally be more useful when allowOverlap is TRUE.

Breaking out of spread events

There are 4 ways for the spread to "stop" spreading. Here, each "event" is defined as all cells that are spawned from a single starting loci. So, one spread call can have multiple spreading "events". The ways outlines below are all acting at all times, i.e., they are not mutually exclusive. Therefore, it is the user's responsibility to make sure the different rules are interacting with each other correctly. Using spreadProb or maxSize are computationally fastest, sometimes dramatically so.

The spread function does not return the result of this stopRule. If, say, an event has both circleMaxRadius and stopRule, and it is the circleMaxRadius that caused the event spreading to stop, there will be no indicator returned from this function that indicates which rule caused the stop.

stopRule has many use cases. One common use case is evaluating a neighbourhood around a focal set of points. This provides, therefore, an alternative to the terra::buffer() function or terra::focal() function. In both of those cases, the window/buffer size must be an input to the function. Here, the resulting size can be emergent based on the incremental growing and calculating of the landscape values underlying the spreading event.

stopRuleBehavior

This determines how the stopRule should be implemented. Because spreading occurs outwards in concentric circles or shapes, one cell width at a time, there are 4 possible ways to interpret the logical inequality defined in stopRule. In order of number of cells included in resulting events, from most cells to fewest cells:

Note

dqrng version 0.4.0 changed the default RNG. If backwards compatibility is needed, set dqrng::dqRNGkind("Xoroshiro128+") before running spread to ensure numerical reproducibility with previous versions.

Author(s)

Eliot McIntire and Steve Cumming

See Also

spread2() for a different implementation of the same algorithm. It is more robust, meaning, there will be fewer unexplainable errors, and the behaviour has been better tested, so it is more likely to be exactly as described under all argument combinations. Also, rings() which uses spread but with specific argument values selected for a specific purpose. terra::distance(). cir() to create "circles"; it is fast for many small problems.

Examples

library(terra)
origDTThreads <- data.table::setDTthreads(2L)
origNcpus <- options(Ncpus = 2L)

# Make random forest cover map
set.seed(123)
emptyRas <- rast(ext(0, 1e2, 0, 1e2), res = 1)
hab <- randomPolygons(emptyRas, numTypes = 40)
names(hab) <- "hab"
mask <- rast(emptyRas)
values(mask) <- 0
mask[1:5000] <- 1
numCol <- ncol(emptyRas)
numCell <- ncell(emptyRas)
directions <- 8

# Can use transparent as a colour
coltab(hab) <- paste(c("transparent", grey(0:40/40)))

terra::plot(hab)

# initiate 10 fires
startCells <- as.integer(sample(1:ncell(emptyRas), 100))
fires <- spread(hab, loci = startCells, 0.235, persistence = 0, numNeighs = 2,
                mask = NULL, maxSize = 1e8, directions = 8, iterations = 1e6, id = TRUE)

terra::plot(hab, type = "classes", legend = FALSE)
fires[fires == 0] <- NA
terra::plot(fires, add = TRUE, col = "red", type = "continuous", legend = FALSE)

# Instead, to give a colour to the zero values, use \code{zero.color=}
coltab(fires) <- NULL
# need to specify "type" to get correct legend
terra::plot(fires,  col = c(colorRampPalette(c("blue", "green"))(100)),
            type = "continuous")

##------------------------------------------------------------------------------
## Continue event by passing interrupted object into spreadState
##------------------------------------------------------------------------------

## Interrupt a spread event using iterations - need `returnIndices = TRUE` to
##  use outputs as new inputs in next iteration
fires <- spread(hab, loci = as.integer(sample(1:ncell(hab), 10)),
                returnIndices = TRUE, 0.235, 0, NULL, 1e8, 8, iterations = 3, id = TRUE)
fires[, list(size = length(initialLocus)), by = id]  # See sizes of fires

fires2 <- spread(hab, loci = NA_real_, returnIndices = TRUE, 0.235, 0, NULL,
                 1e8, 8, iterations = 2, id = TRUE, spreadState = fires)
# NOTE events are assigned arbitrary IDs, starting at 1


## Use data.table and loci...
fires <- spread(hab, loci = as.integer(sample(1:ncell(hab), 10)),
                returnIndices = TRUE, 0.235, 0, NULL, 1e8, 8,
                iterations = 2, id = TRUE)
fullRas <- rast(hab)
fullRas[] <- 1:ncell(hab)
burned <- fires[active == FALSE]
burnedMap <- rasterizeReduced(burned, fullRas, "id", "indices")

terra::plot(burnedMap, type = "classes")

####################
## stopRule examples
####################

# examples with stopRule, which means that the eventual size is driven by the values on the raster
#  passed in to the landscape argument. It won't be exact because the pixel values
#  will likely not allow it
stopRule22 <- function(landscape) sum(landscape) > 100

set.seed(1234)
stopRule1 <- function(landscape) sum(landscape) > 50
stopRuleA <- spread(hab, loci = as.integer(sample(1:ncell(hab), 10)), 1, 0, NULL,
                    maxSize = 1e6, 8, 1e6, id = TRUE, circle = TRUE, stopRule = stopRule1,
                    stopRuleBehavior = "excludePixel")
tapply(hab[], stopRuleA[], sum) # all below 50

set.seed(1234)
# using stopRuleBehavior = "excludePixel"
stopRuleB <- spread(hab, loci = as.integer(sample(1:ncell(hab), 10)), 1, 0, NULL,
                    maxSize = 1e6, 8, 1e6, id = TRUE, circle = TRUE, stopRule = stopRule22,
                    stopRuleBehavior = "excludePixel")
tapply(hab[], stopRuleB[], sum) # all below 100

if (interactive())
  terra::plot(c(stopRuleA, stopRuleB))
# Cellular automata shapes
# Diamonds - can make them with: a boolean raster, directions = 4,
#    stopRule in place, spreadProb = 1
diamonds <- spread(hab > 0, spreadProb = 1, directions = 4, id = TRUE, stopRule = stopRule22)

terra::plot(diamonds)

# Squares - can make them with: a boolean raster, directions = 8,
#    stopRule in place, spreadProb = 1
squares <- spread(hab > 0, spreadProb = 1, directions = 8, id = TRUE, stopRule = stopRule22)
terra::plot(squares)

# Interference shapes - can make them with: a boolean raster, directions = 8,
#    stopRule in place, spreadProb = 1
stopRule2 <- function(landscape) sum(landscape) > 200
squashedDiamonds <- spread(hab > 0, spreadProb = 1,
                           loci = (ncell(hab) - ncol(hab)) / 2 + c(4, -4),
                           directions = 4, id = TRUE, stopRule = stopRule2)
terra::plot(squashedDiamonds)

# Circles with spreadProb < 1 will give "more" circular shapes, but definitely not circles
stopRule2 <- function(landscape) sum(landscape) > 200
seed <- sample(1e4, 1)
set.seed(seed)

circlish <- spread(hab > 0, spreadProb = 1, iterations = 10,
                   loci = (ncell(hab) - ncol(hab)) / 2 + c(4, -4),
                   directions = 8, id = TRUE, circle = TRUE)#, stopRule = stopRule2)
if (interactive())
  terra::plot(c(circlish))

set.seed(seed)
regularCA <- spread(hab > 0, spreadProb = 1, iterations = 10,
                    loci = (ncell(hab) - ncol(hab)) / 2 + c(4, -4),
                    directions = 8, id = TRUE)#, stopRule = stopRule2)

if (interactive()) # compare to circlish
  terra::plot(regularCA)

####################
# complex stopRule
####################

initialLoci <- sample(seq_len(ncell(hab)), 2)
endSizes <- seq_along(initialLoci) * 200

# Can be a function of landscape, id, and/or any other named
#   variable passed into spread
stopRule3 <- function(landscape, id, endSizes) sum(landscape) > endSizes[id]

set.seed(1)
twoCirclesDiffSize <- spread(hab, spreadProb = 1, loci = initialLoci,
                             circle = TRUE, directions = 8, id = TRUE,
                             stopRule = stopRule3, endSizes = endSizes,
                             stopRuleBehavior = "excludePixel")

# or using named list of named elements:
set.seed(1)
twoCirclesDiffSize2 <- spread(hab, spreadProb = 1, loci = initialLoci,
                              circle = TRUE, directions = 8, id = TRUE,
                              stopRule = stopRule3,
                              vars = list(endSizes = endSizes),
                              stopRuleBehavior = "excludePixel")

compareGeom(twoCirclesDiffSize, twoCirclesDiffSize2, res = TRUE,
            stopOnError = FALSE)

terra::plot(c(twoCirclesDiffSize , twoCirclesDiffSize2))

cirs <- values(twoCirclesDiffSize)
vals <- tapply(hab[][cirs > 0], cirs[cirs > 0], sum) # one is <200, other is <400 as per endSizes

# Stop if sum of landscape is big or mean of quality is too small
quality <- rast(hab)
quality[] <- runif(ncell(quality), 0, 1)
stopRule4 <- function(landscape, quality, cells) {
  (sum(landscape) > 20) | (mean(values(quality)[cells]) < 0.3)
}

twoCirclesDiffSize <- spread(hab, spreadProb = 1, loci = initialLoci, circle = TRUE,
                             directions = 8, id = TRUE, stopRule = stopRule4,
                             quality = quality, stopRuleBehavior = "excludePixel")


## Using alternative algorithm, not probabilistic diffusion
## Will give exactly correct sizes, yet still with variability
## within the spreading (i.e., cells with and without successes)
seed <- sample(1e6, 1)
set.seed(seed)
startCells <- startCells[1:4]
maxSizes <- rexp(length(startCells), rate = 1 / 500)
fires <- spread(hab, loci = startCells, 1, persistence = 0,
                neighProbs = c(0.5, 0.5, 0.5) / 1.5,
                mask = NULL, maxSize = maxSizes, directions = 8,
                iterations = 1e6, id = TRUE, plot.it = FALSE, exactSizes = TRUE)
all(table(fires[fires > 0][]) == floor(maxSizes))

terra::plot(fires)
hist(fires[][fires[] > 0], main = "fire size distribution")

## Example with relativeSpreadProb ... i.e., a relative probability spreadProb
##  (shown here because because spreadProb raster is not a probability).
##  Here, we force the events to grow, choosing always 2 neighbours,
##  according to the relative probabilities contained on hab layer.
##
## Note: `neighProbs = c(0,1)` forces each active pixel to move to 2 new pixels
## (`prob = 0` for 1 neighbour, `prob = 1` for 2 neighbours)
##
## Note: set hab3 to be very distinct probability differences, to detect spread
##  differences
hab3 <- (hab < 20) * 200 + 1
seed <- 643503
set.seed(seed)
sam <- sample(which(hab3[] == 1), 1)
set.seed(seed)
events1 <- spread(hab3, spreadProb = hab3, loci = sam, directions = 8,
                  neighProbs = c(0, 1), maxSize = c(70), exactSizes = TRUE)

# Compare to absolute probability version
set.seed(seed)
events2 <- spread(hab3, id = TRUE, loci = sam, directions = 8,
                  neighProbs = c(0, 1), maxSize = c(70), exactSizes = TRUE)

terra::plot(events1)

terra::plot(events2, col = c("white", "red", "red"))

hist(events1[], breaks = 30, main = "Event size distribution") ## TODO: fix this plot
# Compare outputs -- should be more high value hab pixels spread to in event1
#  (randomness may prevent this in all cases)
sum(hab3[events1[] > 0]) >= sum(hab3[events2[] > 0]) ## should be usually TRUE

# clean up
data.table::setDTthreads(origDTThreads)
options(Ncpus = origNcpus)

Simulate a contagious spread process on a landscape, with data.table internals

Description

This can be used to simulate fires, seed dispersal, calculation of iterative, concentric, symmetric (currently) landscape values and many other things. Essentially, it starts from a collection of cells (start, called "events") and spreads to neighbours, according to the directions and spreadProb with modifications due to other arguments. NOTE: the spread function is similar, but sometimes slightly faster, but less robust, and more difficult to use iteratively.

Usage

spread2(
  landscape,
  start = ncell(landscape)/2 - ncol(landscape)/2,
  spreadProb = 0.23,
  persistProb = NA_real_,
  asRaster = TRUE,
  maxSize,
  exactSize,
  directions = 8L,
  iterations = 1000000L,
  returnDistances = FALSE,
  returnDirections = FALSE,
  returnFrom = FALSE,
  maxRetriesPerID = 10,
  spreadProbRel = NA_real_,
  plot.it = FALSE,
  circle = FALSE,
  asymmetry = NA_real_,
  asymmetryAngle = NA_real_,
  allowOverlap = 0,
  neighProbs = NA_real_,
  oneNeighbourOnly = FALSE,
  skipChecks = FALSE
)

Arguments

Details

There are 2 main underlying algorithms for active cells to "spread" to nearby cells (adjacent cells): spreadProb and neighProb. Using spreadProb, every "active" pixel will assess all neighbours (either 4 or 8, depending on directions), and will "activate" whichever neighbours successfully pass independent calls to runif(1,0,1) < spreadProb. The algorithm will iterate again and again, each time starting from the newly "activated" cells. Several built-in decisions are as follows.

1. no active cell can activate a cell that was already activated by the same event (i.e., "it won't go backwards"). 2. If allowOverlap is FALSE, then the previous rule will also apply, regardless of which "event" caused the pixels to be previously active.

This function can be interrupted before all active cells are exhausted if the iterations value is reached before there are no more active cells to spread2 into. The interrupted output (a data.table) can be passed subsequently as an input to this same function (as start). This is intended to be used for situations where external events happen during a spread2 event, or where one or more arguments to the spread2 function change before a spread2 event is completed. For example, if it is desired that the spreadProb change before a spread2 event is completed because, for example, a fire is spreading, and a new set of conditions arise due to a change in weather.

asymmetry here is slightly different than in the spread function, so that it can deal with a RasterLayer of asymmetryAngle. Here, the spreadProb values of a given set of neighbours around each active pixel are adjusted to create adjustedSpreadProb which is calculated maintain the following two qualities:

mean(spreadProb) = mean(ajustedSpreadProb)

and

max(spreadProb)/min(spreadProb) = asymmetry

along the axis of asymmetryAngle. NOTE: this means that the 8 neighbours around an active cell may not fulfill the preceeding equality if asymmetryAngle is not exactly one of the 8 angles of the 8 neighbours. This means that

max(spreadProb)/min(spreadProb)

will generally be less than asymmetry, for the 8 neighbours. The exact adjustment to the spreadProb is calculated with:

angleQuality <- (cos(angles - rad2(asymmetryAngle))+1)/2

which is multiplied to get an angle-adjusted spreadProb:

spreadProbAdj <- actualSpreadProb * angleQuality

which is then rescaled:

adjustedSpreadProb = (spreadProbAdj - min(spreadProbAdj)) * par2 + par1

, where par1 and par2 are parameters calculated internally to make the 2 conditions above true.

If maxSize or exactSize are used, then spreading will continue and stop before or at maxSize or at exactSize, respectively. If iterations is specified, then the function will end, and the returned data.table may (if maxSize) or will (if exactSize) have at least one active cell per event that did not already achieve maxSize or exactSize. This will be very useful to build new, customized higher-level wrapper functions that iteratively call spread2.

Value

Either a data.table (asRaster=FALSE) or a RasterLayer (asRaster=TRUE, the default). The data.table will have one attribute named spreadState, which is a list containing a data.table of current cluster-level information about the spread events. If asRaster=TRUE, then the data.table (with its spreadState attribute) will be attached to the Raster as an attribute named pixel as it provides pixel-level information about the spread events.

The RasterLayer represents every cell in which a successful spread2 event occurred. For the case of, say, a fire this would represent every cell that burned. If allowOverlap is TRUE, the return will always be a data.table.

If asRaster is FALSE, then this function returns a data.table with 3 (or 4 if returnFrom is TRUE) columns:

The attribute saved with the name "spreadState" (e.g., attr(output, "spreadState")) includes a data.table with columns:

and several other objects that provide significant speed ups in iterative calls to spread2. If the user runs spread2 iteratively, there will likely be significant speed gains if the data.table passed in to start should have the attribute attached, or re-attached if it was lost, e.g., via setattr(outInput, "spreadState", attr(out, "spreadState")), where out is the returned data.table from the previous call to spread2, and outInput is the modified data.table. Currently, the modified data.table must have the same order as out.

Breaking out of spread2 events

There are 3 ways for the spread2 to "stop" spreading. Here, each "event" is defined as all cells that are spawned from each unique start location. The ways outlined below are all acting at all times, i.e., they are not mutually exclusive. Therefore, it is the user's responsibility to make sure the different rules are interacting with each other correctly.

Building custom spreading events

This function can be used iteratively, with relatively little overhead compared to using it non-iteratively. In general, this function can be called with arguments set as user needs, and with specifying e.g., iterations = 1. This means that the function will spread outwards 1 iteration, then stop. The returned object will be a data.table or RasterLayer that can be passed immediately back as the start argument into a subsequent call to spread2. This means that every argument can be updated at each iteration.

When using this function iteratively, there are several things to keep in mind. The output will likely be sorted differently than the input (i.e., the order of start, if a vector, may not be the same order as that returned). This means that when passing the same object back into the next iteration of the function call, maxSize or exactSize may not be in the same order. To get the same order, the easiest thing to do is sort the initial start objects by their pixel location, increasing. Then, of course, sorting any vectorized arguments (e.g., maxSize) accordingly.

NOTE: the data.table or RasterLayer should not be altered when passed back into spread2.

allowOverlap

If 1 (or TRUE), then individual events can overlap with one another, i.e., allow overlap between events. If 2 (or NA), then each pixel is essentially independent, allowing overlap between and within events. This likely requires a user to intervene as it is possible to spread back onto itself. If 3 (did not exist previously), individual events can overlap, and there can be overlap within an event, but only within an iteration, i.e., once an iteration is finished, and a pixel was activated, then the spreading will not return onto these pixels. If 0 (or FALSE), then once a pixel is activated, it cannot be re-activated, within or between event. This allows events to not interfere with one another i.e., they do not interact (this is slower than if allowOverlap = FALSE). Default is 0. In the case of 2 or 3, this would be, perhaps, useful for dispersal of, say, insect swarms.

Note

exactSize may not be achieved if there aren't enough cells in the map. Also, exactSize may not be achieved because the active cells are "stuck", i.e., they have no inactivated cells to move to; or the spreadProb is low. In the latter two cases, the algorithm will retry again, but it will only retry from the last iteration's active cells. The algorithm will only retry 10 times before quitting. Currently, there will also be an attempt to "jump" up to four cells away from the active cells to try to continue spreading.

A common way to use this function is to build wrappers around this, followed by iterative calls in a while loop. See example.

Author(s)

Eliot McIntire and Steve Cumming

See Also

spread() for a different implementation of the same algorithm. spread is less robust but it is often slightly faster.

Examples

library(terra)

origDTThreads <- data.table::setDTthreads(2L)
origNcpus <- options(Ncpus = 2L)

a <- rast(ext(0, 10, 0, 10), res = 1)
sams <- sort(sample(ncell(a), 3))

# Simple use -- similar to spread(...)
out <- spread2(a, start = sams, 0.225)
if (interactive()) {
  terra::plot(out)
}

# Use maxSize -- this gives an upper limit
maxSizes <- sort(sample(1:10, size = length(sams)))
out <- spread2(a, start = sams, 0.225, maxSize = maxSizes, asRaster = FALSE)
# check TRUE using data.table .N
out[, .N, by = "initialPixels"]$N <= maxSizes

# Use exactSize -- gives an exact size, if there is enough space on the Raster
exactSizes <- maxSizes
out <- spread2(a, start = sams, spreadProb = 0.225,
               exactSize = exactSizes, asRaster = FALSE)
out[, .N, by = "initialPixels"]$N == maxSizes # should be TRUE TRUE TRUE

# Use exactSize -- but where it can't be achieved
exactSizes <- sort(sample(100:110, size = length(sams)))
out <- spread2(a, start = sams, 1, exactSize = exactSizes)

# Iterative calling -- create a function with a high escape probability
spreadWithEscape <- function(ras, start, escapeProb, spreadProb) {
  out <- spread2(ras, start = sams, spreadProb = escapeProb, asRaster = FALSE)
  while (any(out$state == "sourceActive")) {
    # pass in previous output as start
    out <- spread2(ras, start = out, spreadProb = spreadProb,
                    asRaster = FALSE, skipChecks = TRUE) # skipChecks for speed
  }
  out
}

set.seed(421)
out1 <- spreadWithEscape(a, sams, escapeProb = 0.25, spreadProb = 0.225)
set.seed(421)
out2 <- spread2(a, sams, 0.225, asRaster = FALSE)
# The one with high escape probability is larger (most of the time)
NROW(out1) > NROW(out2) ## TODO: not true

## Use neighProbs, with a spreadProb that is a RasterLayer
# Create a raster of different values, which will be the relative probabilities
#   i.e., they are rescaled to relative probabilities within the 8 neighbour choices.
#   The neighProbs below means 70% of the time, 1 neighbour will be chosen,
#   30% of the time 2 neighbours.
#   The cells with spreadProb of 5 are 5 times more likely than cells with 1 to be chosen,
#   when they are both within the 8 neighbours
sp <- rast(ext(0, 3, 0, 3), res = 1, vals = 1:9) #small raster, simple values
# Check neighProbs worked
out <- list()

# enough replicates to see stabilized probabilities
for (i in 1:100) {
  out[[i]] <- spread2(sp, spreadProbRel = sp, spreadProb = 1,
                      start = 5, iterations = 1,
                      neighProbs = c(1), asRaster = FALSE)
}
out <- data.table::rbindlist(out)[pixels != 5] # remove starting cell
table(sp[out$pixels])
# should be non-significant -- note no 5 because that was the starting cell
#  This tests whether the null model is true ... there should be proportions
#  equivalent to 1:2:3:4:6:7:8:9 ... i.e,. cell 9 should have 9x as many events
#  spread to it as cell 1. This comes from sp object above which is providing
#  the relative spread probabilities
keep <- c(1:4, 6:9)
chisq.test(keep, unname(tabulate(sp[out$pixels]$lyr.1, 9)[keep]),
           simulate.p.value = TRUE)

## Example showing asymmetry
sams <- ncell(a) / 4 - ncol(a) / 4 * 3
circs <- spread2(a, spreadProb = 0.213, start = sams,
                 asymmetry = 2, asymmetryAngle = 135,
                 asRaster = TRUE)
## ADVANCED: Estimate spreadProb when using asymmetry, such that the expected
##   event size is the same as without using asymmetry

  if (interactive()) {
    # Still using `raster` as it works more easily with parallelism due to not using pointers
    #   This will updated at a later release
    if (requireNamespace("raster", quietly = TRUE)) {
      ras <- raster::raster(a)
      ras[] <- 1

      n <- 100
      sizes <- integer(n)
      for (i in 1:n) {
        circs <- spread2(ras, spreadProb = 0.225,
                         start = round(ncell(ras) / 4 - ncol(ras) / 4 * 3),
                         asRaster = FALSE)
        sizes[i] <- circs[, .N]
      }
      goalSize <- mean(sizes)

      if (requireNamespace("DEoptim", quietly = TRUE)) {
        library(parallel)
        library(DEoptim)

        # need 10 cores for 10 populations in DEoptim
        cl <- makeCluster(pmin(10, detectCores() - 2))
        parallel::clusterEvalQ(cl, {
          library(SpaDES.tools)
          library(terra)
          library(raster)
          library(fpCompare)
        })

        objFn <- function(sp, n = 20, ras, goalSize) {
          sizes <- integer(n)
          for (i in 1:n) {
            circs <- spread2(ras, spreadProb = sp, start = ncell(ras) / 4 - ncol(ras) / 4 * 3,
                             asymmetry = 2, asymmetryAngle = 135,
                             asRaster = FALSE)
            sizes[i] <- circs[, .N]
          }
          abs(mean(sizes) - goalSize)
        }
        aa <- DEoptim(objFn, lower = 0.2, upper = 0.23,
                      control =
                        DEoptim.control(
                          cluster = cl, NP = 10, VTR = 0.02,
                          # imposing itermax simply for example; should let go to completion
                          itermax = 5,
                          initialpop = as.matrix(rnorm(10, 0.213, 0.001))),
                      ras = ras, goalSize = goalSize)

        # The value of spreadProb that will give the
        #    same expected event sizes to spreadProb = 0.225 is:
        sp <- aa$optim$bestmem
        circs <- spread2(ras, spreadProb = sp, start = ncell(ras) / 4 - ncol(ras) / 4 * 3,
                         asymmetry = 2, asymmetryAngle = 135, asRaster = FALSE)

        stopCluster(cl)
      }
    }
  }

  # clean up
  data.table::setDTthreads(origDTThreads)
  options(Ncpus = origNcpus)

An alternative spread function, conceived for insects

Description

This is built with spread2() and is still experimental. This one differs from other attempts in that it treats the advection and dispersal as mathematical vectors that are added together. They are "rounded" to pixel centres.

Usage

spread3(
  start,
  rasQuality,
  rasAbundance,
  advectionDir,
  advectionMag,
  meanDist,
  dispersalKernel = "exponential",
  sdDist = 1,
  plot.it = 2,
  minNumAgents = 50,
  verbose = getOption("LandR.verbose", 0),
  saveStack = NULL,
  skipChecks = FALSE
)

Arguments

Value

A data.table with all information used during the spreading

Examples

## these tests are fairly heavy, so don't run during automated tests
#########################################################
# Simple case, no variation in rasQuality, numeric advectionDir and advectionMag
#########################################################

  library(terra)

  origDTThreads <- data.table::setDTthreads(2L)
  origNcpus <- options(Ncpus = 2L)

  maxDim <- 10000
  ras <- terra::rast(terra::ext(c(0, maxDim, 0, maxDim)), res = 100, vals = 0)
  rasQuality <- terra::rast(ras)
  rasQuality[] <- 1
  rasAbundance <- terra::rast(rasQuality)
  rasAbundance[] <- 0
  # startPixel <- middlePixel(rasAbundance)
  startPixel <- sample(seq(terra::ncell(rasAbundance)), 30)
  rasAbundance[startPixel] <- 1000
  advectionDir <- 70
  advectionMag <- 4 * res(rasAbundance)[1]
  meanDist <- 2600

  # Test the dispersal kernel -- create a function
  plotDispersalKernel <- function(out, meanAdvectionMag) {
    out[, disGroup := round(distance / 100) * 100]
    freqs <- out[, .N, by = "disGroup"]
    freqs[, `:=`(cumSum = cumsum(N), N = N)]
    plot(freqs$disGroup, freqs$cumSum, # addTo = "CumulativeNumberSettled",
         main = "Cumulative Number Settled") # can plot the distance X number
    abline(v = meanAdvectionMag + meanDist)
    newTitle <- "Number Settled By Distance"
    plot(freqs$disGroup, freqs$N, # addTo = gsub(" ", "", newTitle),
         main = newTitle) # can plot the distance X number
    abline(v = meanAdvectionMag + meanDist)
    # should be 0.63:
    freqs[disGroup == meanAdvectionMag + meanDist, cumSum] / tail(freqs, 1)[, cumSum]
    mtext(side = 3, paste("Average habitat quality: ",
                          round(mean(rasQuality[], na.rm = TRUE), 2)),
          outer = TRUE, line = -2, cex = 2)
  }
  out <- spread3(rasAbundance = rasAbundance,
                 rasQuality = rasQuality,
                 advectionDir = advectionDir,
                 advectionMag = advectionMag,
                 meanDist = meanDist, verbose = 2,
                 plot.it = interactive())

  plotDispersalKernel(out, advectionMag)

  # The next examples are potentially time consuming; avoid on automated testing
  if (interactive()) {
    #########################################################
    ### The case of variable quality raster
    #########################################################
    rasQuality <- terra::rast(system.file("extdata", "rasQuality.tif",
                                          package = "SpaDES.tools"))
    terra::crs(rasQuality) <- system.file("extdata", "targetCRS.rds", package = "SpaDES.tools") |>
      readRDS() |>
      slot("projargs")
    mask <- rasQuality < 5
    rasQuality[mask[] %in% TRUE] <- 0
    # rescale so min is 0.75 and max is 1
    rasQuality[] <- rasQuality[] / (reproducible::maxFn(rasQuality) * 4) + 1 / 4
    rasAbundance <- terra::rast(rasQuality)
    rasAbundance[] <- 0
    startPixel <- sample(seq(ncell(rasAbundance)), 300)
    rasAbundance[startPixel] <- 1000
    advectionDir <- 75
    advectionMag <- 4 * res(rasAbundance)[1]
    meanDist <- 2600
    out <- spread3(rasAbundance = rasAbundance,
                   rasQuality = rasQuality,
                   advectionDir = advectionDir,
                   advectionMag = advectionMag,
                   meanDist = meanDist, verbose = 2,
                   plot.it = interactive())
    if (interactive()) {
      plotDispersalKernel(out, advectionMag)
    }

    ###############################################################################
    ### The case of variable quality raster, raster for advectionDir & advectionMag
    ###############################################################################
    maxDim <- 10000
    ras <- terra::rast(terra::ext(c(0, maxDim, 0, maxDim)), res = 100, vals = 0)
    rasQuality <- terra::rast(ras)
    rasQuality[] <- 1
    rasAbundance <- terra::rast(rasQuality)
    rasAbundance[] <- NA
    # startPixel <- middlePixel(rasAbundance)
    startPixel <- sample(seq(ncell(rasAbundance)), 25)
    rasAbundance[startPixel] <- 1000

    # raster for advectionDir
    advectionDir <- terra::rast(system.file("extdata", "advectionDir.tif",
                                            package = "SpaDES.tools"))
    crs(advectionDir) <- crs(rasQuality)
    # rescale so min is 0.75 and max is 1
    advectionDir[] <- advectionDir[] / (reproducible::maxFn(advectionDir)) * 180

    # raster for advectionMag
    advectionMag <- terra::rast(system.file("extdata", "advectionMag.tif",
                                            package = "SpaDES.tools"))
    crs(advectionMag) <- crs(rasQuality)
    # rescale so min is 0.75 and max is 1
    advectionMag[] <- advectionMag[] / (reproducible::maxFn(advectionMag)) * 600

    out <- spread3(rasAbundance = rasAbundance,
                   rasQuality = rasQuality,
                   advectionDir = advectionDir,
                   advectionMag = advectionMag,
                   meanDist = meanDist, verbose = 2,
                   plot.it = interactive())

    if (interactive()) {
      names(advectionDir) <- "Wind direction"
      names(advectionMag) <- "Wind speed"
      names(rasAbundance) <- "Initial abundances"
      terra::plot(c(advectionDir, advectionMag, rasAbundance))

      plotDispersalKernel(out, mean(advectionMag[]))
    }

    #########################################
    # save iterations to a stack to make animated GIF
    ########################################
    tmpStack <- tempfile(pattern = "stackToAnimate", fileext = ".tif")
    out <- spread3(rasAbundance = rasAbundance,
                   rasQuality = rasQuality,
                   advectionDir = advectionDir,
                   advectionMag = advectionMag,
                   meanDist = 2600, verbose = 2,
                   plot.it = interactive(), saveStack = tmpStack)

    ## This animates the series of images into an animated GIF
    if (require(animation, quietly = TRUE)) {
      out2 <- terra::rast(tmpStack)
      gifName <- file.path(tempdir(), "animation.gif")

      # Only works on some systems; may need to configure
      # Works on Windows without system adjustments
      if (identical(.Platform$OS.type, "windows"))
        saveGIF(interval = 0.1, movie.name = gifName, expr = {
          for (i in seq(length(names(out2)))) terra::plot(out2[[i]])
        })
    }
  }

  # clean up
  data.table::setDTthreads(origDTThreads)
  options(Ncpus = origNcpus)

Test that metadata of 2 or more objects is the same

Description

Currently, only Raster class has a useful method. Defaults to all(sapply(list(...)[-1], function(x) identical(list(...)[1], x)))

Usage

testEquivalentMetadata(...)

Arguments

SELES - Transitioning to next time step

Description

Describes the probability of an agent successfully persisting until next time step. THIS IS NOT YET FULLY IMPLEMENTED.

A SELES-like function to maintain conceptual backwards compatibility with that simulation tool. This is intended to ease transitions from SELES.

You must know how to use SELES for these to be useful.

Usage

transitions(p, agent)

Arguments

Value

Returns new ⁠SpatialPoints*⁠ object with potentially fewer agents.

Author(s)

Eliot McIntire

Wrap coordinates or pixels in a torus-like fashion

Description

Generally useful for model development purposes. Primarily used internally in e.g., crw if torus = TRUE.

Usage

wrap(X, bounds, withHeading = FALSE)

Arguments

Details

If withHeading used, then X must be an sf or SpatVector object that contains two columns, x1 and y1, with the immediately previous agent locations.

Value

Object of the same class as X, but with coordinates updated to reflect the wrapping.

Author(s)

Eliot McIntire

Examples

origDTThreads <- data.table::setDTthreads(2L)
origNcpus <- options(Ncpus = 2L)

xrange <- yrange <- c(-50, 50)
hab <- terra::rast(terra::ext(c(xrange, yrange)))
hab[] <- 0

# initialize agents
N <- 10

# previous points
x1 <- y1 <- rep(0, N)
# initial points
starts <- cbind(x = stats::runif(N, xrange[1], xrange[2]),
                y = stats::runif(N, yrange[1], yrange[2]))

# create the agent object # the x1 and y1 are needed for "previous location"
agent <- terra::vect(data.frame(x1, y1, starts), geom = c("x", "y"))

ln <- rlnorm(N, 1, 0.02) # log normal step length
sd <- 30 # could be specified globally in params

if (interactive()) {
  # clearPlot()
  terra::plot(hab, col = "white")
}

for (i in 1:10) {
  agent <- crw(agent = agent, extent = terra::ext(hab), stepLength = ln,
               stddev = sd, lonlat = FALSE, torus = FALSE) # don't wrap
  if (interactive()) terra::plot(agent[, 1], add = TRUE, col = 1:10)
}
terra::crds(agent) # many are "off" the map, i.e., beyond the extent of hab
agent <- SpaDES.tools::wrap(agent, bounds = terra::ext(hab))
terra::plot(agent, add = TRUE, col = 1:10) # now inside the extent of hab

# clean up
data.table::setDTthreads(origDTThreads)
options(Ncpus = origNcpus)