Stochastic outward spread, implemented in C++

A fast, deliberately narrow alternative to spread() for the common case: a per-cell probability of being spread to, one or more starting cells, an optional per-fire maximum size, and cell indices back.

spreadCpp(
  landscape,
  loci,
  spreadProb,
  maxSize = Inf,
  directions = 8L,
  iterations = Inf,
  minSize = 0,
  jumpTries = 0L,
  jumpMeanDist = 3,
  minSizeTries = 100L,
  stopCells = NULL,
  stopEvent = NULL,
  stopAt = NULL
)

Arguments

landscape: A SpatRaster; only its geometry (number of columns and cells) is used.
loci: Integer vector of starting cell indices, one per fire.
spreadProb: Numeric of length 1 or terra::ncell(landscape), or a SpatRaster or RasterLayer congruent with landscape: the probability that a cell is spread to.
maxSize: Numeric of length 1 or length(loci): the largest number of cells each fire may reach. Inf (the default) means no limit.
directions: 4 or 8. Default 8.
iterations: Maximum number of generations. Default Inf.
minSize: Numeric of length 1 or length(loci): a size each fire reaches before the normal generation rule applies. While a fire has fewer than minSize cells, all of its burning cells stay active from one generation to the next (persistence) and draw again, with the same spreadProb, against their still-unburned neighbours, so the fire grows into a patch shaped by spreadProb rather than dying out early. The generation that reaches minSize stops there, exactly. The next generation spreads from every cell of the patch that still has an unburned, burnable neighbour (its perimeter), not only from the cells that caught last; after that only the cells that caught in the last generation spread, as usual. A fire under minSize with no unburned neighbour whose spreadProb is above 0 (and not NA) stops where it is; otherwise only iterations limits how long a very low spreadProb takes to get there. Must not exceed maxSize. The default, 0, is the behaviour without it, with the same random draws. Used to start every escaped fire at the escape size, e.g. minSize = 9 cells for 50 ha at 5.76-ha cells.
jumpTries: Integer, default 0 (off). Only for a fire that is still under minSize and stuck – none of its burning cells has an unburned neighbour with spreadProb above 0 – up to jumpTries attempts are made to jump: a random cell of the fire, a distance drawn from an exponential with mean jumpMeanDist, truncated to 1.5-20 cells so it clears the adjacent ring, and a uniform direction. A target off the landscape, already burned or unburnable is rejected; otherwise it catches with its spreadProb. The first success joins the fire and the fire keeps growing from all its cells. If every attempt fails the fire stops where it is. With 0 no extra random draws are made. This is a fixed rule for fires trapped in small patches, not a fitted spotting process; spotting tied to fire weather is future work.
jumpMeanDist: Mean jump distance in cells (see jumpTries). Default 3.
minSizeTries: Integer, default 100. With minSize above 1, minSize is reached by rejection rather than persistence: each fire spreads normally from its one ignition cell, and a fire that dies below minSize is undone (its cells freed, its rows dropped) and started again from its ignition cell, up to minSizeTries times. A fire that reaches minSize is kept and goes on uncut, so its size, shape and burning front are those of a fire that got there under its own spreadProb – a sample from the fires that reach minSize. After minSizeTries rejections the fire falls back to the persistence rule of minSize (and jumpTries) from its ignition cell. The expected number of tries per fire is 1 / P(reaching minSize), so a low spreadProb or a small burnable patch costs more; the result carries an attribute "minSizeTries", a data.table with each fire's id, its rejected tries, and whether it used the fallback. minSizeTries = 0 gives persistence alone (the behaviour before this argument existed, with the same random draws). Without minSize it has no effect.
stopCells, stopEvent, stopAt: Optional per-fire stop rule; all NULL (the default) leaves the results and the random numbers used exactly as without them. stopCells is a vector of cell indices (indexed like loci), stopEvent the same length: the fire (an index into loci, 1-based) each stop cell belongs to, and stopAt an integer of length 1 or length(loci): the number of a fire's own stop cells at which it stops. A cell may be a stop cell of several fires, and repeated (cell, fire) pairs count once. Only a fire's own stop cells count towards its stopAt, and a fire with no stop cells never stops by this rule. The fire stops the moment it burns its stopAt-th stop cell, in the middle of a generation, in the same way as at maxSize: the cells already burned stay burned and the other fires carry on. It interacts with minSize and jumpTries as maxSize does: a fire that has stopped is never pushed on to minSize or made to jump, and a try rejected under minSizeTries is undone with its count. The ignition cell counts if it is a stop cell. Giving stopCells needs stopEvent and stopAt.

Value

A data.table keyed on id, with integer columns id, initialLocus and indices, and a logical active (always FALSE, since the spread has finished) – the same shape spread() returns for returnIndices = TRUE. When minSize is above 1 and minSizeTries above 0, it also carries an attribute "minSizeTries": a data.table with one row per fire (in loci order) and columns id, tries (its rejected attempts) and fallback (TRUE if it ran out of tries and used persistence). The attribute is transient: most data.table operations drop it, e.g. subsetting or dt[, .N, by = "id"], so read it from the returned object straight away, e.g. attr(out, "minSizeTries").

The rules it follows

Growth happens in generations. The cells that caught in one generation are the only ones that try to spread in the next, so fires grow outwards and a fire's extent is tied to how many generations it has had.
One draw per (burning cell, unburned neighbour) pair, compared against spreadProb of the neighbour, not of the cell doing the spreading. A cell with k burning neighbours therefore catches with probability 1 - (1 - p)^k in that generation, exactly as in spread().
A cell belongs to exactly one fire, and is never spread to twice. Once a cell has caught it is never reconsidered. Burning cells are visited in a random order within a generation, so when two fires reach the same cell neither is systematically favoured.
A fire stops at its own maxSize and never exceeds it. spread() adds a whole generation and then randomly drops the excess; spreadCpp() stops on the boundary, which lands on maxSize exactly.
NA in spreadProb means unburnable, as in spread().

Limitations

spreadCpp() is deliberately narrow. It covers one case – a per-cell probability of being spread to, one or more ignitions, an optional per-fire size cap – and takes only the arguments listed above. Every other argument spread() accepts is unsupported, and passing one is an error rather than being silently ignored:

allowOverlap: Not supported. Fires never overlap; a cell belongs to one fire.
returnDistances, circle, circleMaxRadius: Not supported. No distances are computed and spread is not constrained to a disc.
asymmetry, asymmetryAngle: Not supported. Spread is isotropic; there is no directional bias.
neighProbs, relativeSpreadProb: Not supported. Every eligible neighbour is drawn for independently; the number of neighbours a cell spreads to is not itself drawn from a distribution, and probabilities are absolute rather than rescaled within a cell's neighbourhood.
stopRule, stopRuleBehavior: Not supported. The only stopping conditions are maxSize, iterations, and having nowhere left to go.
exactSizes: Not supported, though maxSize is exact in the sense that a fire which reaches its cap stops there.
persistence: Not supported. A burned cell stays burned and does not re-burn.
mask: Not supported. Use NA (or 0) in spreadProb to make cells unburnable.
spreadState: Not supported. A spread cannot be continued from an earlier one; each call starts from loci.
spreadProbLater: Not supported. One spreadProb applies to every generation.
torus: Not supported. The landscape does not wrap; spread stops at the edges.
plot.it, id, returnIndices: Not supported as options. Nothing is plotted, and the return value is always the indices form.

Use spread() when any of these are needed.

Reproducibility

spreadCpp() is not a drop-in reimplementation of spread(). It makes its own random draws, so for a given seed it does not reproduce spread()'s cells, and swapping one for the other changes results and invalidates cached output. It is reproducible in the ordinary sense: the same seed and the same inputs give the same answer.

What the two do agree on is behaviour in aggregate. Over 300 seeds on a 1.44M-cell landscape the two produce the same total burned area (Wilcoxon p = 0.34) and the same fire-size distribution (Kolmogorov-Smirnov p = 0.56), with quantiles within 7% of each other from the 10th percentile to the 99th.

Performance

Roughly 2x faster than spread() across a range of fire regimes: 1.6-2.4x over nine scenarios spanning 596 to 188,762 burned cells and landscapes from 0.36M to 9M cells, median 2.1x.

The ratio understates the benefit on the expensive calls, so judge it in seconds. Across those nine scenarios the speedup RATIO falls slightly as more burns (correlation -0.37 with burned cells) while the seconds saved rise steeply (correlation +0.95). The scenario with the lowest ratio, 161,701 burned cells at 1.58x, saves 0.048 s per call – 24 times more than the 2.00x scenario at the other end, which saves 0.002 s. On real fire-model fits the same holds: two landscapes whose ratios differ by 30% (1.44x and 1.88x) save almost identical wall time per objective evaluation (9.5 s and 9.6 s).

The gain depends on what fraction of the landscape burns, and it can vanish. The per-cell state is allocated once per call and costs the same in both implementations, so it is pure overhead that neither avoids; what spreadCpp() makes cheaper is the spreading itself. When almost nothing burns, that overhead dominates and there is nothing left to win. Measured on one landscape of ~8.4M cells: with 0.13% of it burned, spreadCpp() was no faster at all (0.95x), where scenarios burning 0.7% and 4.2% of their landscapes gave 2.1x and 2.4x.

The same allocation makes cost grow with landscape size even when the fires do not: holding the fires fixed, going from 1.44M to 9M cells costs an extra 0.035 s per call. So if the fires occupy a small part of a large landscape, crop first and map the indices back – that both speeds up the call and restores the advantage. Cropping 1.44M to 300k cells is worth about 30%.

Examples

if (requireNamespace("terra", quietly = TRUE)) {
  r <- terra::rast(terra::ext(0, 100, 0, 100), resolution = 1, vals = 0)
  set.seed(1)
  out <- spreadCpp(r, loci = c(2500L, 7500L), spreadProb = 0.23,
                   maxSize = c(50, 100))
  out[, .N, by = "id"]
}
#> Key: <id>
#>       id     N
#>    <int> <int>
#> 1:     1     4
#> 2:     2    10