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.runiverse.dev>). 
Authors:  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 <[email protected]> 
License:  GPL3 
Version:  2.0.7.9000 
Built:  20240912 06:08:53 UTC 
Source:  https://github.com/PredictiveEcology/SpaDES.tools 
SpaDES.tools
packageSpatial 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:
adj() 
An optimized (i.e., faster) version of
terra::adjacent()

cir() 
Identify pixels in a circle around a
SpatialPoints* object 
directionFromEachPoint() 
Fast calculation of direction and distance surfaces 
distanceFromEachPoint() 
Fast calculation of distance surfaces 
rings() 
Identify rings around focal cells (e.g., buffers and donuts) 
spokes() 
TO DO: need description 
spread() 
Contagious cellular automata 
wrap() 
Create a torus from a grid 
Agents have several methods and functions specific to them:
crw() 
Simple correlated random walk function 
heading() 
Determines the heading between SpatialPoints*

quickPlot::makeLines() 
Makes SpatialLines object for, e.g., drawing arrows 
move() 
A meta function that can currently only take "crw" 
specificNumPerPatch() 
Initiate a specific number of agents per patch 
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 GISrelated functions:
quickPlot::equalExtent() 
Assess whether a list of extents are all equal 
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.
rasterizeReduced() 
Convert reduced representation to full raster 
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.
randomPolygons() 
Creates a random polygon with specified number of classes. 
See the NLMR package for tools to generate random landscapes (rasters).
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:
agentLocation() 
Agent location 
initiateAgents() 
Initiate agents into a SpatialPointsDataFrame

numAgents() 
Number of agents 
probInit() 
Probability of initiating an agent or event 
transitions() 
Transition probability 
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
.
Maintainer: Alex M Chubaty [email protected] (ORCID)
Authors:
Eliot J B McIntire [email protected] (ORCID)
Other contributors:
Yong Luo [email protected] [contributor]
Ceres Barros [email protected] (ORCID) [contributor]
Steve Cumming [email protected] [contributor]
Jean Marchal [email protected] [contributor]
His Majesty the King in Right of Canada, as represented by the Minister of Natural Resources Canada [copyright holder]
Useful links:
Report bugs at https://github.com/PredictiveEcology/SpaDES.tools/issues
These have been written with speed in mind.
.pointDistance( from, to, angles = NA, maxDistance = NA_real_, otherFromCols = FALSE )
.pointDistance( from, to, angles = NA, maxDistance = NA_real_, otherFromCols = FALSE )
from 
Numeric matrix with 2 or 3 or more columns. They must include x and y,
representing x and y coordinates of "from" cell. If there is a column
named "id", it will be "id" from 
to 
Numeric matrix with 2 or 3 columns (or optionally more, all of which
will be returned),
x and y, representing x and y coordinates of "to" cells, and
optional "id" which will be matched with "id" from 
angles 
Logical. If 
maxDistance 
Numeric in units of number of cells. The algorithm will build
the whole surface (from 
otherFromCols 
other columns to use as 'from' 
adjacent
function, and Just In Time compiled versionFaster 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
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 )
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 )
x 

cells 
vector of cell numbers for which adjacent cells should be found. Cell numbers start with 1 in the upperleft corner and increase from left to right and from top to bottom. 
directions 
the number of directions in which cells should be connected:
4 (rook's case), 8 (queen's case), or 
sort 
logical. Whether the outputs should be sorted or not, using cell ids
of the 
pairs 
logical. If 
include 
logical. Should the focal cells be included in the result? 
target 
a vector of cells that can be spread to. This is the inverse of a mask. 
numCol 
numeric indicating number of columns in the raster.
Using this with 
numCell 
numeric indicating number of cells in the raster.
Using this with 
match.adjacent 
logical. Should the returned object be the same as

cutoff.for.data.table 
numeric. If the number of cells is above this value, the function uses data.table which is faster with large numbers of cells. Default is 5000, which appears to be the turning point where data.table becomes faster. 
torus 
Logical. Should the spread event wrap around to the other side of the raster?
Default is 
id 
numeric If not 
numNeighs 
A numeric scalar, indicating how many neighbours to return. Must be
less than or equal to 
returnDT 
A logical. If TRUE, then the function will return the result
as a 
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:
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)
don't make intermediate objects; just put calculation into return statement
The steps used in the algorithm are:
Calculate indices of neighbouring cells
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)
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
.
Eliot McIntire
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)
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 initiationSets 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.
agentLocation(map)
agentLocation(map)
map 
A 
Object of same class as provided as input.
If a Raster*
, then zeros are converted to NA
.
Eliot McIntire
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
.
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 )
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 )
landscape 
Raster on which the circles are built. 
coords 
Either a matrix with 2 (or 3) columns, 
loci 
Numeric. An alternative to 
maxRadius 
Numeric vector of length 1 or same length as 
minRadius 
Numeric vector of length 1 or same length as 
allowOverlap 
Logical. Should duplicates across id be removed or kept. Default TRUE. 
allowDuplicates 
Logical. Should duplicates within id be removed or kept. Default FALSE. This is useful if the actual x, y coordinates are desired, rather than the cell indices. This will increase the size of the returned object. 
includeBehavior 
Character string. Currently accepts only 
returnDistances 
Logical. If 
angles 
Numeric. Optional vector of angles, in radians, to use. This will create
"spokes" outward from 
returnAngles 
Logical. If 
returnIndices 
Logical or numeric. If 
closest 
Logical. When determining nonoverlapping circles, should the function
give preference to the closest 
simplify 
logical. If 
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.
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.
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
.
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)
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)
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.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.
.cirSpecialQuick(landscape, loci, maxRadius, minRadius)
.cirSpecialQuick(landscape, loci, maxRadius, minRadius)
landscape 
Raster on which the circles are built. 
loci 
Numeric. An alternative to 
maxRadius 
Numeric vector of length 1 or same length as 
minRadius 
Numeric vector of length 1 or same length as 
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 pairwise point distance.
As a result, this can be used for doing inverse distance weightings, seed rain,
cumulative effects of distancebased 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.
distanceFromEachPoint( from, to = NULL, landscape, angles = NA_real_, maxDistance = NA_real_, cumulativeFn = NULL, distFn = function(dist) 1/(1 + dist), cl, ... )
distanceFromEachPoint( from, to = NULL, landscape, angles = NA_real_, maxDistance = NA_real_, cumulativeFn = NULL, distFn = function(dist) 1/(1 + dist), cl, ... )
from 
Numeric matrix with 2 or 3 or more columns. They must include x and y,
representing x and y coordinates of "from" cell. If there is a column
named "id", it will be "id" from 
to 
Numeric matrix with 2 or 3 columns (or optionally more, all of which
will be returned),
x and y, representing x and y coordinates of "to" cells, and
optional "id" which will be matched with "id" from 
landscape 

angles 
Logical. If 
maxDistance 
Numeric in units of number of cells. The algorithm will build
the whole surface (from 
cumulativeFn 
A function that can be used to incrementally accumulate
values in each 
distFn 
A function. This can be a function of 
cl 
A cluster object. Optional. This would generally be created using

... 
Any additional objects needed for 
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.
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
.
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).
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)
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)
.duplicatedInt
does same as duplicated
in R, but only on integers, and faster.
It uses Rcpp sugar
duplicatedInt(x)
duplicatedInt(x)
x 
Integer Vector 
A logical vector, as per duplicated
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
).
dwrpnorm2(theta, mu, rho, sd = 1, acc = 1e05, tol = acc)
dwrpnorm2(theta, mu, rho, sd = 1, acc = 1e05, tol = acc)
theta 
value at which to evaluate the density function, measured in radians. 
mu 
mean direction of distribution, measured in radians. 
rho 
mean resultant length of distribution. 
sd 
different way of select rho, see details below. 
acc 
parameter defining the accuracy of the estimation of the density. Terms are added to the infinite summation that defines the density function until successive estimates are within acc of each other. 
tol 
same as 
Eliot McIntire
# 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)
# 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.fastCrop
is deprecated.
fastCrop(x, y, ...)
fastCrop(x, y, ...)
x 
Raster to crop 
y 
Raster to crop with 
... 
other 
velox::VeloxRaster_crop
raster
of a random Gaussian process.Defunct.
gaussMap( x, scale = 10, var = 1, speedup = 1, method = "RMexp", alpha = 1, inMemory = FALSE, ... )
gaussMap( x, scale = 10, var = 1, speedup = 1, method = "RMexp", alpha = 1, inMemory = FALSE, ... )
x 
A spatial object (e.g., a 
scale 
The spatial scale in map units of the Gaussian pattern. 
var 
Spatial variance. 
speedup 
An numeric value indicating how much faster than 'normal' to generate maps. It may be necessary to give a value larger than 1 for large maps. Default is 1. 
method 
The type of model used to produce the Gaussian pattern.
Should be one of 
alpha 
A required parameter of the 
inMemory 
Should the 
... 
Additional arguments to 
A raster map with same extent as x
, with a Gaussian random pattern.
Determines the heading between spatial points.
heading(from, to)
heading(from, to)
from 
The starting position; an object of class 
to 
The ending position; an object of class 
The heading between the points, in degrees.
Eliot McIntire
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)
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)
SELES
 Initiate agentsSets 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.
initiateAgents(map, numAgents, probInit, asSpatialPoints = TRUE, indices)
initiateAgents(map, numAgents, probInit, asSpatialPoints = TRUE, indices)
map 

numAgents 
numeric resulting from a call to 
probInit 
a 
asSpatialPoints 
logical or 
indices 
numeric. Indices of where agents should start 
A SpatialPointsDataFrame
, with each row representing an individual agent
Eliot McIntire
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) }
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) }
[a,b]
Default values of a=0; b=1
allow for quick test if
x
is a probability.
inRange(x, a = 0, b = 1)
inRange(x, a = 0, b = 1)
x 
values to be tested 
a 
lower bound (default 0) 
b 
upper bound (default 1) 
Logical vectors. NA
values in x
are retained.
Alex Chubaty
set.seed(100) x < stats::rnorm(4) # 0.50219235 0.13153117 0.07891709 0.88678481 inRange(x, 0, 1)
set.seed(100) x < stats::rnorm(4) # 0.50219235 0.13153117 0.07891709 0.88678481 inRange(x, 0, 1)
RasterLayer
(s)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)
.
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" )
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" )
x 
A list of split raster tiles (i.e., from 
fun 
Function (e.g. 
r 
The raster to be split. 
nx 
The number of tiles to make along the xaxis. 
ny 
The number of tiles to make along the yaxis. 
buffer 
Numeric vector of length 2 giving the size of the buffer along the x and y axes.
If values greater than or equal to 
path 
Character specifying the directory to which the split tiles will be saved. If missing, the function will write to memory. 
cl 
A cluster object. Optional. This would generally be created using

rType 
Data type of the split rasters. Defaults to FLT4S. 
fExt 
file extension (e.g., 
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 parallelaware 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()
.
mergeRaster
returns a RasterLayer
object.
splitRaster
returns a list (length nx*ny
) of cropped raster tiles.
Yong Luo, Alex Chubaty, Tati Micheletti & Ian Eddy
Alex Chubaty and Yong Luo
terra::merge()
, terra::mosaic()
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(), "splitRasterexample"), 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)
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(), "splitRasterexample"), 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
middlePixel(ras)
middlePixel(ras)
ras 
A 
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.
Wrapper for selecting different animal movement methods.
This version uses just turn angles and step lengths to define the correlated random walk.
move(hypothesis = "crw", ...) crw( agent, extent, stepLength, stddev, lonlat = FALSE, torus = FALSE, returnMatrix = FALSE )
move(hypothesis = "crw", ...) crw( agent, extent, stepLength, stddev, lonlat = FALSE, torus = FALSE, returnMatrix = FALSE )
hypothesis 
Character vector, length one, indicating which movement
hypothesis/method to test/use. Currently defaults to
'crw' (correlated random walk) using 
... 
arguments passed to the function in 
agent 
A 
extent 
An optional 
stepLength 
Numeric vector of length 1 or number of agents describing step length. 
stddev 
Numeric vector of length 1 or number of agents describing standard deviation of wrapped normal turn angles. 
lonlat 
Logical. If 
torus 
Logical. Should the movement be wrapped to the opposite
side of the map, as determined by the 
returnMatrix 
If 
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.
A SpatVector
points object with updated spatial position defined
by a single occurrence of step length(s) and turn angle(s).
Eliot McIntire
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:725736.
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)
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)
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.
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"), ... )
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"), ... )
x 
A 
pad 
Integer. Number of cells by which to pad 
type 
One of the supported 
... 
Further arguments passed to 
A RasterLayer
/SpatRaster
with same extent as x
, with randomly generated values.
nlm_mpd
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) }
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) }
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.
numAgents(N, probInit)
numAgents(N, probInit)
N 
Number of agents to initiate (integer scalar). 
probInit 
Probability of initializing an agent at the location. 
A numeric, indicating number of agents to start
Eliot McIntire
Patch size
patchSize(patches)
patchSize(patches)
patches 
Number of patches. 
SELES
 Probability of InitiationDescribes 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.
probInit(map, p = NULL, absolute = NULL)
probInit(map, p = NULL, absolute = NULL)
map 
A 
p 
probability, provided as a numeric or raster 
absolute 
logical. Is 
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).
Eliot McIntire
SpatRaster
of random polygonsThese 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
.
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)
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)
ras 
A raster that whose extent will be used for the random polygons. 
numTypes 
Numeric value. The number of unique polygon types to use. 
... 
Other arguments passed to 
x 
Either a 
hectares 
Deprecated. Use 
area 
A numeric, the approximate area in 
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
.
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) } }
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) } }
SpaDES
modulesCreate default study areas for use with SpaDES
modules
randomStudyArea(center = NULL, size = 10000, seed = NULL)
randomStudyArea(center = NULL, size = 10000, seed = NULL)
center 

size 
Numeric specifying the approximate size of the area in m^2.
Default 
seed 
Numeric indicating the random seed to set internally (useful for ensuring the same study area is produced each time). 
SpatVector
a < randomStudyArea(seed = 123) if (interactive()) { terra::plot(a) }
a < randomStudyArea(seed = 123) if (interactive()) { terra::plot(a) }
Convert reduced representation to full raster
rasterizeReduced( reduced, fullRaster, newRasterCols, mapcode = names(fullRaster), ... )
rasterizeReduced( reduced, fullRaster, newRasterCols, mapcode = names(fullRaster), ... )
reduced 

fullRaster 

newRasterCols 
Character vector, length 1 or more, with the name(s) of
the column(s) in 
mapcode 
a character, length 1, with the name of the column in 
... 
Other arguments. None used yet. 
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
Eliot McIntire
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)
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)
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).
rings( landscape, loci = NA_real_, id = FALSE, minRadius = 2, maxRadius = 5, allowOverlap = FALSE, returnIndices = FALSE, returnDistances = TRUE, ... )
rings( landscape, loci = NA_real_, id = FALSE, minRadius = 2, maxRadius = 5, allowOverlap = FALSE, returnIndices = FALSE, returnDistances = TRUE, ... )
landscape 
A 
loci 
A vector of locations in 
id 
Logical. If 
minRadius 
Numeric. Minimum radius to be included in the ring.
Note: this is inclusive, i.e., 
maxRadius 
Numeric. Maximum radius to be included in the ring.
Note: this is inclusive, i.e., 
allowOverlap 
Logical. If 
returnIndices 
Logical or numeric. If 
returnDistances 
Logical. Should the function include a column with the
individual cell distances from the locus where that event
started. Default is 
... 
Any other argument passed to 
This will return a data.table
with columns as described in
spread
when returnIndices = TRUE
.
Eliot McIntire
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
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)
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)
Slightly faster than runif, and used a lot
runifC(N)
runifC(N)
N 
Integer Vector 
A vector of uniform random numbers as per runif
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
.
specificNumPerPatch(patches, numPerPatchTable = NULL, numPerPatchMap = NULL)
specificNumPerPatch(patches, numPerPatchTable = NULL, numPerPatchMap = NULL)
patches 

numPerPatchTable 
A 
numPerPatchMap 
A 
A raster with 0s and 1s, where the 1s indicate starting locations of agents following the numbers above.
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)
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)
This is a generalized version of a notion of a viewshed. The main difference is that there can be many "viewpoints".
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, ... )
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, ... )
landscape 
Raster on which the circles are built. 
coords 
Either a matrix with 2 (or 3) columns, 
loci 
Numeric. An alternative to 
maxRadius 
Numeric vector of length 1 or same length as 
minRadius 
Numeric vector of length 1 or same length as 
allowOverlap 
Logical. Should duplicates across id be removed or kept. Default TRUE. 
stopRule 
A function. If the spokes are to stop. This can be a function
of 
includeBehavior 
Character string. Currently accepts only 
returnDistances 
Logical. If 
angles 
Numeric. Optional vector of angles, in radians, to use. This will create
"spokes" outward from 
nAngles 
Numeric, length one. Alternative to angles. If provided, the function
will create a sequence of angles from 
returnAngles 
Logical. If 
returnIndices 
Logical or numeric. If 
... 
Objects to be used by 
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
.
Eliot McIntire
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)
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)
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.
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, ... )
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, ... )
landscape 
A 
loci 
A vector of locations in 
spreadProb 
Numeric, 
persistence 
A length 1 probability that an active cell will continue to burn, per time step. 
mask 

maxSize 
Numeric. Maximum number of cells for a single or all events to be spread.
Recycled to match 
directions 
The number of adjacent cells in which to look; default is 8 (Queen case). Can only be 4 or 8. 
iterations 
Number of iterations to spread.
Leaving this 
lowMemory 
Deprecated. 
returnIndices 
Logical or numeric. If 
returnDistances 
Logical. Should the function include a column with the
individual cell distances from the locus where that event
started. Default is 
mapID 
Deprecated. Use 
id 
Logical. If 
plot.it 
If 
spreadProbLater 
Numeric, or 
spreadState 

circle 
Logical. If 
circleMaxRadius 
Numeric. A further way to stop the outward spread of events.
If 
stopRule 
A function which will be used to assess whether each
individual cluster should stop growing.
This function can be an argument of 
stopRuleBehavior 
Character. Can be one of 
allowOverlap 
Logical. If 
asymmetry 
A numeric indicating the ratio of the asymmetry to be used.
Default is 
asymmetryAngle 
A numeric indicating the angle in degrees (0 is "up",
as in North on a map), that describes which way the 
quick 
Logical. If 
neighProbs 
A numeric vector, whose sum is 1.
It indicates the probabilities an individual spread iteration
spreading to 
exactSizes 
Logical. If 
relativeSpreadProb 
Logical. If 
... 
Additional named vectors or named list of named vectors required for 
For large rasters, a combination of lowMemory = TRUE
and
returnIndices = TRUE
or returnIndices = 2
will be fastest and use the least amount of memory.
20220725: 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)))
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:
id 
an arbitrary ID 1:length(loci) identifying
unique clusters of spread events, i.e., all cells
that have been spread into that have a
common initial cell. 
initialLocus 
the initial cell number of that particular spread event. 
indices 
The cell indices of cells that have been touched by the spread algorithm. 
active 
a logical indicating whether the cell is active (i.e., could still be a source for spreading) or not (no spreading will occur from these cells). 
This will generally be more useful when allowOverlap
is TRUE
.
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.
spreadProb 
Probabilistically, if spreadProb is low enough, active spreading events will stop. In practice, active spreading events will stop. In practice, this number generally should be below 0.3 to actually see an event stop 
maxSize 
This is the number of cells that are "successfully" turned on during a spreading event. This can be vectorized, one value for each event 
circleMaxRadius 
If circle is TRUE, then this will be the maximum
radius reached, and then the event will stop. This is
vectorized, and if length is >1, it will be matched
in the order of loci

stopRule 
This is a function that can use "landscape", "id", "cells",
or any named vector passed into spread in the ... .
This can take on relatively complex functions.
Passing in, say, a RasterLayer to spread
can access the individual values on that arbitrary
RasterLayer using "cells".
These will be calculated within all the cells of the individual
event (equivalent to a "group_by(event)" in dplyr .
So, sum(arbitraryRaster[cells]) would sum up all
the raster values on the arbitraryRaster raster
that are overlaid by the individual event.
This can then be used in a logical statement. See examples.
To confirm the cause of stopping, the user can assess the values
after the function has finished. 
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:
"includeRing" 
Will include the entire ring of cells that, as a group,
caused stopRule to be TRUE . 
"includePixel" 
Working backwards from the entire ring that caused the
stopRule to be TRUE , this will iteratively
random cells in the final ring
until the stopRule is FALSE . This will add back
the last removed cell and include it in the return result
for that event. 
"excludePixel" 
Like "includePixel" , but it will not add back the cell
that causes stopRule to be TRUE

"excludeRing" 
Analogous to "excludePixel" , but for the entire final
ring of cells added. This will exclude the entire ring of cells
that caused the stopRule to be TRUE

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.
Eliot McIntire and Steve Cumming
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.
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)
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)
data.table
internalsThis 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.
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 )
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 )
landscape 
Required. A 
start 
Required. Either a vector of pixel numbers to initiate spreading, or a

spreadProb 
Numeric of length 1 or length 
persistProb 
Numeric of length 1 or 
asRaster 
Logical, length 1. If 
maxSize 
Numeric. Maximum number of cells for a single or all events to be 
exactSize 
Numeric vector, length 1 or 
directions 
The number adjacent cells in which to look; default is 8 (Queen case). Can only be 4 or 8. 
iterations 
Number of iterations to 
returnDistances 
Logical. Should the function include a column with the individual cell distances from the locus where that event started. Default is FALSE. See Details. 
returnDirections 
Logical. Should the function include a column with the individual directions (in radians) from the locus where that event started. Default is FALSE. 
returnFrom 
Logical. Should the function return a column with the source, i.e, the lag 1 "from" pixel, for each iteration. 
maxRetriesPerID 
Only active if 
spreadProbRel 
Optional 
plot.it 
If 
circle 
Logical. If TRUE, then outward 
asymmetry 
A numeric or 
asymmetryAngle 
A numeric or 
allowOverlap 

neighProbs 
An optional numeric vector, whose sum is 1.
It indicates the probabilities that an individual spread iteration
will spread to 
oneNeighbourOnly 
Logical. Default is 
skipChecks 
Logical. If TRUE, the argument checking (i.e., assertions) will be skipped. This should likely only be used once it is clear that the function arguments are well understood and function speed is of the primary importance. This is likely most useful in repeated iteration cases i.e., if this call is using the previous output from this same function. 
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 builtin decisions are as follows.
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 angleadjusted 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 higherlevel wrapper functions that
iteratively call spread2
.
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 clusterlevel 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 pixellevel 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:
initialPixels 
the initial cell number of that particular
spread2 event. 
pixels 
The cell indices of cells that have
been touched by the spread2 algorithm. 
state 
a logical indicating whether the cell is active (i.e., could still be a source for spreading) or not (no spreading will occur from these cells). 
from 
The pixel indices that were the immediately preceding
"source" for each pixels , i.e., the lag 1 pixels.
Only returned if returnFrom is TRUE 
The attribute saved with the name "spreadState"
(e.g., attr(output, "spreadState")
)
includes a data.table
with columns:
id 
An arbitrary code, from 1 to length(start) for each "event". 
initialPixels 
the initial cell number of that particular
spread2 event. 
numRetries 
The number of restarts the event did because it got
stuck (normally only because exactSize was used
and was not achieved. 
maxSize 
The number of pixels that were provided as inputs via
maxSize or exactSize . 
size 
The current size, in pixels, of each event. 
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 reattached 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
.
spread2
eventsThere 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.
spreadProb 
Probabilistically, if spreadProb is low enough,
active spreading events will stop.
In practice, this number generally should be below 0.3
to actually see an event stop. 
maxSize 
This is the number of cells that are "successfully" turned
on during a spreading event. spreadProb will still
be active, so, it is possible that the end size of each event
is smaller than maxSize , but they will not be greater
than maxSize

exactSize 
This is the number of cells that are "successfully" turned
on during a spreading event. This will override an event that
stops probabilistically via spreadProb , but forcing
its last set of active cells to try again to find neighbours.
It will try maxRetriesPerID times per event, before giving up.
During those maxRetriesPerID times, it will try to "jump" up to
4 cells outwards from each of the active cells, every 5 retries. 
iterations 
This is a hard cap on the number of internal iterations to
complete before returning the current state of the system
as a data.table . 
This function can be used iteratively, with relatively little overhead compared to using
it noniteratively. 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 reactivated, 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.
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.
Eliot McIntire and Steve Cumming
spread()
for a different implementation of the same algorithm.
spread
is less robust but it is often slightly faster.
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 nonsignificant  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)
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 nonsignificant  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)
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.
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 )
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 )
start 
Raster indices from which to initiate dispersal 
rasQuality 
A raster with habitat quality. Currently, must be scaled from 0 to 1, i.e., a probability of "settling" 
rasAbundance 
A raster where each pixel represents the number of "agents" or pseudoagents contained. This number of agents, will be spread horizontally, and distributed from each pixel that contains a nonzero non NA value. 
advectionDir 
A single number or 
advectionMag 
A single number or 
meanDist 
A single number indicating the mean distance parameter in map units
(not pixels), for a negative exponential distribution
dispersal kernel (e.g., 
dispersalKernel 
One of either 
sdDist 
A single number indicating the 
plot.it 
Numeric. With increasing numbers above 0, there will be plots produced during iterations. Currently, only 0, 1, or 2+ are distinct. 
minNumAgents 
Single numeric indicating the minimum number of agents to consider all dispersing finished. Default is 50. 
verbose 
Numeric. With increasing numbers above 0, there will be more messages produced. Currently, only 0, 1, or 2+ are distinct. 
saveStack 
If provided as a character string, it will save each iteration
as part of a 
skipChecks 
Logical. If 
A data.table
with all information used during the spreading
## 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)
## 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)
Currently, only Raster class has a useful method. Defaults to
all(sapply(list(...)[1], function(x) identical(list(...)[1], x)))
testEquivalentMetadata(...)
testEquivalentMetadata(...)
... 
2 or more of the same type of object to test for equivalent metadata. 
SELES
 Transitioning to next time stepDescribes 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.
transitions(p, agent)
transitions(p, agent)
p 
realized probability of persisting (i.e., either 0 or 1). 
agent 

Returns new SpatialPoints*
object with potentially fewer agents.
Eliot McIntire
Generally useful for model development purposes. Primarily used internally
in e.g., crw
if torus = TRUE
.
wrap(X, bounds, withHeading = FALSE)
wrap(X, bounds, withHeading = FALSE)
X 

bounds 
Either a 
withHeading 
logical. If 
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.
Object of the same class as X
, but with coordinates updated to reflect the wrapping.
Eliot McIntire
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)
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)