R/spread2.R
spread2.Rd
This can be used to simulate fires, seed dispersal, calculation of iterative,
concentric, symmetric (currently) landscape values and many other things.
Essentially, it starts from a collection of cells (start
, called "events")
and spreads to neighbours, according to the directions
and spreadProb
with modifications due to other arguments. NOTE:
the spread
function is similar, but sometimes slightly faster, but less
robust, and more difficult to use iteratively.
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_, skipChecks = FALSE )
landscape  Required. A 

start  Required. Either a vector of pixel numbers to initiate spreading, or a
data.table that is the output of a previous 
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 spread2. Recycled to match 
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 spread2.
Leaving this 
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 TRUE, then plot the raster at every iteration, so one can watch the spread2 event grow. 
circle  Logical. If TRUE, then outward spread2 will be by equidistant rings,
rather than solely by adjacent cells (via 
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 
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. 
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 preceeding
"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
.
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.
1. no active cell can active 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  rad(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 exactSize
or maxSize
are used, then spreading will continue and stop
before or at maxSize
or at exactSize
. If iterations
is specified,
then the function will end, and the returned data.table
will still
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
.
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 unactivated 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 iterations 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.
There are 3 ways for the spread2 to "stop" spreading.
Here, each "event" is defined as all cells that are spawned from each unique
start
location.
The ways outlined below are all acting at all times, i.e., they are not
mutually exclusive.
Therefore, it is the user's responsibility to make sure the different rules
are interacting with each other correctly.
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 iterations = 1 (say). 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 use 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.
spread
for a different implementation of the same algorithm.
spread
is less robust but it is often slightly faster.
Eliot McIntire and Steve Cumming
library(data.table) library(raster) library(quickPlot) data.table::setDTthreads(1) a < raster(extent(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()) { clearPlot() 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#> logical(0)# 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#> logical(0)# 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)#> [1] 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 < raster(extent(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])#> #> 1 2 3 4 6 7 8 9 #> 3 4 5 10 18 11 22 27# 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], 9)[keep]), simulate.p.value = TRUE)#> #> Pearson's Chisquared test with simulated pvalue (based on 2000 #> replicates) #> #> data: keep and unname(tabulate(sp[out$pixels], 9)[keep]) #> Xsquared = 56, df = NA, pvalue = 1 #>## 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 (FALSE) { if (interactive()) { ras < 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) library(parallel) library(DEoptim) cl < makeCluster(pmin(10, detectCores()  2)) # need 10 cores for 10 populations in DEoptim parallel::clusterEvalQ(cl, { library(SpaDES.tools) 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, initialpop = as.matrix(rnorm(10, 0.213, 0.001))), ras = a, 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) } }