Skip to main content

LandWeb Manual: v. 3.0.0

8 LandMine Module

This documentation is work in progress. Please report any discrepancies or omissions at https://github.com/PredictiveEcology/LandMine/issues.

8.0.0.1 Authors:

Eliot J B McIntire [aut, cre], Alex M. Chubaty [ctb]

8.1 Module Overview

8.1.1 Module summary

Landmine is a model created for simulating the natural range of variation for landscapes in the boreal forest (Andison 1996; Andison 1998). It has been widely used by the public and the private sector for various purposes. This SpaDES module is a rewrite of the fire component in native R.

8.1.2 Model Differences

The current version has not yet been fully tested and compared with the original version, but there are currently several known differences:

  1. Fire sizes are taken from a Truncated Pareto distribution, resulting in numerous very small fires, and few large fires;
  2. Parameters have been fitted to the landscapes that are under study in the LandWeb project.

8.1.3 Known Species

Landmine requires the following codes as inputs (the genus and species codes below), which converts and groups species as follows. Each of the species groups has its own Rate of Spread (ROS) for fire spreading:

Table 8.1: LandMine species codes.
Species Group Code
Jack pine Pine (PINU) Pinu_ban
Lodgepole pine Pine (PINU) Pinu_con
Unspecified pine species Pine (PINU) Pinu_sp
Paper birch Deciduous (DECI) Betu_pap
Balsam poplar Deciduous (DECI) Popu_bal
Trembling aspen Deciduous (DECI) Popu_tre
Larch/Tamarack Deciduous (DECI) Lari_lar
Black spruce Black spruce (PICE_MAR) Pice_mar
White spruce White spruce (PICE_GLA) Pice_gla
Fir species Fir (ABIE) Abie_sp

8.1.4 Module inputs and parameters

Table 8.2 shows the full list of module inputs.

Table 8.2: List of LandMine input objects and their description.
objectName objectClass desc sourceURL
cohortData data.table Columns: B, pixelGroup, speciesCode (as a factor of the names), age. indicating several features about the current vegetation of stand. NA
fireReturnInterval Raster A raster layer that is a factor raster, with at least 1 column called fireReturnInterval, representing the fire return interval in years. NA
pixelGroupMap RasterLayer Pixels with identical values share identical stand features NA
rasterToMatch RasterLayer DESCRIPTION NEEDED NA
ROSTable data.table A data.table with 3 columns, ‘age’, ‘leading’, and ‘ros’. The values under the ‘age’ column can be ‘mature’, ‘immature’, ‘young’ and compound versions of these, e.g., ‘immature_young’ which can be used when 2 or more age classes share same ‘ros’. ‘leading’ should be vegetation type. ‘ros’ gives the rate of spread values for each age and type. NA
rstFlammable Raster A raster layer, with 0, 1 and NA, where 1 indicates areas that are flammable, 0 not flammable (e.g., lakes) and NA not applicable (e.g., masked) NA
rstTimeSinceFire Raster a time since fire raster layer NA
species data.table Columns: species, speciesCode, Indicating several features about species NA
sppColorVect character named character vector of hex colour codes corresponding to each species NA
sppEquiv data.table Multi-columned data.table indicating species name equivalencies. Default taken from LandR::sppEquivalencies_CA which has names for species of trees in Canada NA
studyArea SpatialPolygonsDataFrame multipolygon, typically buffered around an area of interest (i.e., studyAreaReporting) to use for simulation. Defaults to an area in Southwestern Alberta, Canada. NA
studyAreaReporting SpatialPolygonsDataFrame multipolygon (typically smaller/unbuffered than studyArea) to use for plotting/reporting. Defaults to an area in Southwestern Alberta, Canada. NA

Summary of user-visible parameters (Table 8.3).

Table 8.3: List of LandMine parameters and their description.
paramName paramClass default min max paramDesc
biggestPossibleFireSizeHa numeric 1e+06 10000 2e+06 An upper limit, in hectares, of the truncated Pareto distribution of fire sizes
burnInitialTime numeric 1 NA NA This describes the simulation time at which the first burn event should occur
fireTimestep numeric 1 NA NA This describes the simulation time interval between burn events
flushCachedRandomFRI logical FALSE NA NA If no Fire Return Interval map is supplied, then a random one will be created and cached. Use this to make a new one.
maxReburns integer 10 0 20 Number of attempts to reburn fires that don’t reach their target fire size. Unlike maxRetriesPerID, which tries spreading an ongoing fire to nearby cells, this reignites the fire from a new pixel, so it’s less likely to continue being stuck in a region with sinuous fires.
maxRetriesPerID integer 4 0 20 Number of attempts (‘jumps’) that will be made per event ID, before abandoning. See ?SpaDES.tools::spread2. NOTE: increasing this value results in drastically longer simulation times because firelets get ‘stuck’.
minPropBurn numeric 0.9 0 1 Minimum proportion burned pixels to use when triggering warnings about simulated fires.
mixedType numeric 2 1 2 How to define mixed stands: 1 for any species admixture; 2 for deciduous > conifer. See ?vegTypeMapGenerator.
mode character single NA NA use ‘single’ to run part of a landscape simulation; use ‘multi’ to run as part of postprocessing multiple simulation runs.
reps integer NA 1 NA number of replicates/runs per study area when running in ‘multi’ mode.
ROSother integer 30 NA NA default ROS value for non-forest vegetation classes.this is needed when passing a modified ROSTable, e.g. using log-transformed values.
ROStype character default NA NA One of ‘burny’, ‘equal’, ‘log’, or ‘default’.
sppEquivCol character LandR NA NA The column in sim$specieEquivalency data.table to use as a naming convention.
useSeed integer NA NA Only used for creating a starting cohortData dataset. If NULL, then it will be randomly generated; If non-NULL, will pass this value to set.seed and be deterministic and identical each time. WARNING: setting the seed to a specific value will cause all simulations to be identical!
vegLeadingProportion numeric 0.8 0 1 a number that define whether a species is leading for a given pixel
.plotInitialTime numeric 1 NA NA This describes the simulation time at which the first plot event should occur
.plotInterval numeric 1 NA NA This describes the simulation time interval between plot events
.plots character png, screen NA NA Passed to types in Plots (see ?Plots). There are a few plots that are made within this module, if set. Note that plots (or their data) saving will ONLY occur at end(sim). If NA, plotting is turned off completely (this includes plot saving).
.saveInitialTime numeric NA NA NA This describes the simulation time at which the first save event should occur
.saveInterval numeric NA NA NA This describes the simulation time interval between save events
.studyAreaName character NA NA NA Human-readable name for the study area used - e.g., a hash of the studyarea obtained using reproducible::studyAreaName()
.useCache logical FALSE NA NA Should this entire module be run with caching activated? This is generally intended for data-type modules, where stochasticity and time are not relevant
.unitTest logical TRUE NA NA Some functions can have internal testing. This will turn those on or off, if any exist.
.useParallel numeric 2 NA NA Used in burning. Will be passed to data.table::setDTthreads(). NOTE: use .useParallel <= 2 as the additonal RAM overhead too high given marginal speedup.

8.1.5 Module outputs

Description of the module outputs (Table 8.4).

Table 8.4: List of LandMine outputs and their description.
objectName objectClass desc
fireInitialTime numeric The initial event time of the burn event. This is simply a reassignment from P(sim)$burnInitialTime.
fireSizes list A list of data.tables, one per burn event, each with two columns, size and maxSize. These indicate the actual sizes and expected sizes burned, respectively. These can be put into a single data.table with rbindlist(sim$fireSizes, idcol = 'year')
fireReturnInterval RasterLayer A Raster map showing the fire return interval. This is created from the rstCurrentBurn.
fireReturnIntervalsByPolygonNumeric numeric A vector of the fire return intervals, ordered by the numeric representation of polygon ID
fireTimestep numeric The number of time units between successive fire events in a fire module.
friSummary data.table summary fire return interval table
kBest numeric A numeric scalar that is the optimal value of K in the Truncated Pareto distribution (rtruncpareto)
numFiresPerYear numeric The average number of fires per year, by fire return interval level on rstCurrentBurn.
rstCurrentBurn RasterLayer A raster layer, produced at each timestep, where each pixel is either 1 or 0 indicating burned or not burned.
rstCurrentBurnCumulative RasterLayer Cumulative number of times a pixel has burned
sppEquiv data.table Same as input, but with new column, LandMine.

8.2 Usage

To run this Landmine module alone (i.e., for fitting), the following should work (iff raster inputs for rstStudyRegion and rstFlammable are available), assuming all R packages are available.

NB: Paths will have be changed for a different user.

8.2.1 Package dependencies

To determine which packages are used by LandMine, use:

SpaDES.core::packages(modules = "LandMine", paths = "..")[[1]]
##  [1] "SpaDES.core"                                        
##  [2] "assertthat"                                         
##  [3] "data.table"                                         
##  [4] "ggplot2"                                            
##  [5] "grDevices"                                          
##  [6] "gridExtra"                                          
##  [7] "magrittr"                                           
##  [8] "PredictiveEcology/LandR@development (>= 1.1.0.9003)"
##  [9] "PredictiveEcology/pemisc@development"               
## [10] "PredictiveEcology/SpaDES.tools@development"         
## [11] "raster"                                             
## [12] "RColorBrewer"                                       
## [13] "VGAM"

8.2.2 Module usage 

studyArea <- SpaDES.tools::randomStudyArea(seed = 1234, size = 1e+10)
rasterToMatch <- raster(studyArea, res = 250)
times <- list(start = 0, end = 13)

parameters <- list(# LandMine = list(flushCachedRandomFRI = TRUE) ) parameters
parameters <- list(# LandMine = list(flushCachedRandomFRI = TRUE) ) <-
parameters <- list(# LandMine = list(flushCachedRandomFRI = TRUE) ) list(#
parameters <- list(# LandMine = list(flushCachedRandomFRI = TRUE) ) LandMine
parameters <- list(# LandMine = list(flushCachedRandomFRI = TRUE) ) =
parameters <- list(# LandMine = list(flushCachedRandomFRI = TRUE) ) list(flushCachedRandomFRI
parameters <- list(# LandMine = list(flushCachedRandomFRI = TRUE) ) =
parameters <- list(# LandMine = list(flushCachedRandomFRI = TRUE) ) TRUE)
parameters <- list(# LandMine = list(flushCachedRandomFRI = TRUE) ) )

modules <- list("LandMine")

objects <- list(studyArea = studyArea, rasterToMatch = rasterToMatch)

paths <- list(cachePath = cacheDir, modulePath = moduleDir, inputPath = inputDir,
    outputPath = outputDir)
mySim <- simInit(times = times, params = parameters, modules = modules,
    objects = objects, paths = paths)

dev()
mySimOut <- spades(mySim, .plotInitialTime = times$start, debug = TRUE)

8.3 Testing the burn algorithm 

Require(c("data.table", "DEoptim", "parallel", "SDMTools"))

s <- simInit(times = times, params = parameters, modules = modules,
    objects = objects, paths = paths)

ros <- raster(extent(0, 250000, 0, 250000), res = 250, vals = 0)
ros <- ros == 0
fireSize <- 1e+05
source(file.path(moduleDir, "LandMine", "R", "burn.R"))

burnFun <- function(ros, centreCell, fireSize, spawnNewActive,
    sizeCutoffs, burn1, spreadProb) {
    burned <- burn1(landscape = ros, startCells = centreCell,
        fireSizes = fireSize, spreadProb = spreadProb, spreadProbRel = ros,
        spawnNewActive = spawnNewActive, sizeCutoffs = sizeCutoffs)
    burnedMap <- raster(ros)
    burnedMap[] <- NA

    burnedMap[burned$pixels] <- burned$initialPixels
    LM <- SDMTools::PatchStat(burnedMap, cellsize = res(burnedMap)[1])
    list(burnedMap = burnedMap, LM = LM)
}

makeParallel <- function(wantParallel, numClus, cl, funs, addlPkgs) {
    if (wantParallel) {
        library(parallel)
        if (is.null(cl)) {
            cl <- parallel::makeCluster(numClus)
        }
        funs <- c(funs)
        parallel::clusterExport(cl, funs)
        # env <- environment()
        parallel::clusterExport(cl, "addlPkgs", envir = parent.frame())
        parallel::clusterEvalQ(cl, {
            lapply(addlPkgs, library, character.only = TRUE)
        })
        return(cl)
    }
}

wrap <- function(cl = NULL, reps, par, burn1, burnFun, objs,
    addlPkgs) {
    funNames <- c(deparse(substitute(burn1)), deparse(substitute(burnFun)))
    cl <- makeParallel(TRUE, numClus = if (length(cl))
        length(cl) else detectCores() - 1, cl, c(funNames, objs), addlPkgs)
    objs <- append(mget(objs, envir = parent.frame()), list(burn1 = burn1))
    parallel::clusterExport(cl, "objs", envir = parent.frame())
    burnMapList <- parallel::clusterApplyLB(cl, reps, function(r) {
        do.call("burnFun", objs)
    })
    return(list(cl = cl, out = burnMapList))
}

fitSN <- function(sna, ros, centreCell, fireSizes = 10^(2:5),
    desiredPerimeterArea = 0.004, burnFun, burn1) {
    sizeCutoffs <- 10^sna[5:6]
    spreadProb <- sna[7]
    sna <- c(10^(sna[1]), 10^(sna[2]), 10^(sna[3]), 10^(sna[4]))
    bfs1 <- lapply(fireSizes, function(fireSize) {
        burnFun(ros, centreCell, fireSize, sna, sizeCutoffs,
            burn1 = burn1, spreadProb = spreadProb)
    })
    res <- lapply(seq(bfs1), function(bfCount) {
        abs(log(bfs1[[bfCount]]$LM[1, "perim.area.ratio"]) -
            log(desiredPerimeterArea)) + 100 * (sum(bfs1[[bfCount]]$burnedMap[],
            na.rm = TRUE) < fireSizes[bfCount])
        # it needs to get to above 90,000 HA for it to
        # count
    })
    a <- sum(unlist(res))  # * log10(fireSizes)) # weigh larger ones more
    attr(a, "bfs1") <- bfs1
    a
}

fitSN2 <- function(par, ros, centreCell, fireSizes = 10^(2:5),
    desiredPerimeterArea = 0.003, spreadProb = 0.9) {
    sizeCutoffs <- 10^c(par[4], par[5])
    bfs1 <- lapply(fireSizes, function(fireSize) {
        sna <- min(-0.15, par[1] + par[2] * log10(fireSize))
        sna <- 10^c(sna * par[3], sna * 2 * par[3], sna * 3 *
            par[3], sna * 4 * par[3])
        # sna <- -1
        burnFun(ros, centreCell, fireSize, sna, sizeCutoffs,
            burn1, spreadProb)
    })
    res <- lapply(seq(bfs1), function(bfCount) {
        abs(log(bfs1[[bfCount]]$LM[1, "perim.area.ratio"]) -
            log(desiredPerimeterArea)) + 100 * (sum(bfs1[[bfCount]]$burnedMap[],
            na.rm = TRUE) < fireSizes[bfCount])
        # it needs to get to above 90,000 HA for it to
        # count
    })
    a <- sum(unlist(res))
    attr(a, "bfs1") <- bfs1
    a
}

8.4 Optimizing parameters

The following code chunk tries to find values of spawnNewActive that creates “reasonable” fire shapes at all sizes.

wantParallel <- TRUE
maxRetriesPerID <- 4  ## 4 retries (5 attempts total)
spreadProb <- 0.9
spawnNewActive <- c(0.46, 0.2, 0.26, 0.11)
sizeCutoffs <- c(8000, 20000)
# spawnNewActive <- c(0.1, 0.04, 0.025, 0.01) sizeCutoffs
# <- c(8e3, 2e4)

NineCorners <- cellFromRowCol(ros, row = nrow(ros)/4 * rep(1:3,
    3), col = ncol(ros)/4 * rep(1:3, each = 3))

centreCell <- cellFromRowCol(ros, row = nrow(ros)/2, col = ncol(ros)/2)

## Set variables
objs <- c("ros", "centreCell", "fireSize", "spawnNewActive",
    "sizeCutoffs", "spreadProb")
funs <- c("burnFun", "burn1")
addlPkgs <- c("data.table", "raster", "SDMTools", "SpaDES.tools")
### SET UP CLUSTER FOR PARALLEL
numCores <- min(detectCores()/2, 120L)
if (Sys.info()["sysname"] == "Windows") {
    cl <- parallel::makeCluster(numCores)
} else if (grepl("W-VIC-A1053|pinus[.]for-cast[.]ca", Sys.info()["nodename"])) {
    Require("future")

    ## TODO: customize this to make a cluster on several
    ## machines via ssh
    clNames <- c(rep("localhost", 0), rep("picea.for-cast.ca",
        25), rep("pseudotsuga.for-cast.ca", 80))

    ## check the cluster: 1. at least 10 populations per
    ## parameter; 7 params; 2. ensure it's a multiple of
    ## the number of params (7); 3. R can't create more
    ## than ~125 socket connections.
    stopifnot(length(clNames) >= 70, length(clNames)%%7 == 0,
        length(clNames) <= 120)
    cl <- makeClusterPSOCK(clNames, homogeneous = FALSE, verbose = TRUE)

    clusterSetRNGStream(cl, sample(1e+08, 1))
} else {
    cl <- parallel::makeCluster(numCores, type = "FORK")
    clusterSetRNGStream(cl, sample(1e+08, 1))
}

if (!inherits(cl[[1]], "forknode")) {
    parallel::clusterExport(cl, funs)
    parallel::clusterExport(cl, objs)
    # env <- environment()
    parallel::clusterExport(cl, "addlPkgs", envir = parent.frame())
    parallel::clusterEvalQ(cl, {
        lapply(addlPkgs, library, character.only = TRUE)
    })
}
opt_sn <- DEoptim(fitSN, lower = c(-2, -3, -3, -3, 1, 3.5, 0.75),
    upper = c(-0.1, -0.5, -0.5, -1, 3.5, 5, 1), control = DEoptim.control(VTR = 0.001,
        NP = as.integer(length(clNames)), itermax = 200, cluster = cl,
        strategy = 6), ros = ros, centreCell = centreCell, fireSizes = c(10,
        100, 1000, 10000, 1e+05), desiredPerimeterArea = 0.003,
    burnFun = burnFun, burn1 = burn1)

opt_sn$optim$bestmem  ## best param values
# dput(unname(opt_sn$optim$bestmem))

## assign with suffix to facilitate multiple DEoptim runs
timestamp <- format(Sys.time(), "%Y%M%d%H%M%S")
assign(paste0("opt_sn_", timestamp), opt_sn)
saveRDS(opt_sn, file.path(moduleDir, "LandMine", "data", paste0(timestamp,
    "_DEoptim_250m.rds")))
fs_sn <- c(10, 100, 1000, 10000, 1e+05)
fit_sn <- fitSN(sna = c(-1, -1, -1, -2, 2, 4, 0.9), ros = ros,
    centreCell = centreCell, fireSizes = fs_sn, desiredPerimeterArea = 0.003,
    burnFun = burnFun, burn1 = burn1)
bfs1_sn <- purrr::transpose(attr(fit_sn, "bfs1"))
LM_sn <- do.call(rbind, bfs1_sn$LM)
plot(fs_sn, LM_sn[, "perim.area.ratio"])  ## NOTE: visual inspection - not too round; not too sinuous

A second (alternative) version tries the optimization using fewer parameters, to test whether a simpler version gets better/different results. Although this version was not used for the final module, we preserve it here for posterity.

fs_optim2 <- c(0.2, 1:8) * 10000
opt_sn2 <- DEoptim(fitSN2, lower = c(1, -1, 1, 3, 4), upper = c(3,
    -0.3, 3, 4, 5), control = DEoptim.control(VTR = 0.001, itermax = 40,
    cluster = cl, strategy = 6), ros = ros, centreCell = centreCell,
    fireSizes = fs_optim2, desiredPerimeterArea = 0.003)
fs_sn2 <- round(runif(10, 10, 4000))
fit_sn2 <- fitSN2(par = c(2, -0.63333, 1, 3.2, 4.4), ros = ros,
    centreCell = centreCell, fireSizes = fs_sn2, desiredPerimeterArea = 0.003,
    spreadProb = 0.9)
bfs1_sn2 <- purrr::transpose(attr(fit_sn2, "bfs1"))
LM_sn2 <- do.call(rbind, bfs1_sn2$LM)
plot(fs_sn2, LM_sn2[, "perim.area.ratio"])  ## NOTE: visual inspection - not too round; not too sinuous

8.5 Manual inspection of optimization results

8.5.1 Original (2018) version

The original version was run using 100m pixels, despite the simulations being run using 250m pixels. This was corrected and rerun below.

## 10,000 hectares burns gave this
spawnNewActive[2:3] <- c(0.0235999945606232, 0.0263073265505955)

# 100,000 hectare burns gave this spawnNewActive <-
# 10^c(-1.264780, -1.970946, -1.288213, -2.202580)
spawnNewActive <- 10^c(-0.586503645288758, -1.08108837273903,
    -2.14391896536108, -1.00221184641123)
sizeCutoffs <- 10^c(3.37711253212765, 4.52040993282571)

sns <- c(-1.733262, -0.933318, -2.562183, -2.493687, 3.064458,
    4.812305)
spawnNewActive <- 10^sns[1:4]
sizeCutoffs <- 10^sns[5:6]
# spawnNewActive <- 10^c(-1.646419, -1.815395, -2.809013,
# -2.613337) sizeCutoffs <- 10^c(3.888317, 4.641961)

## 100,000 pixel fires -- the next worked, but I think we
## can get better sns <- structure(
## c(-1.652459,-0.962121,-0.964879,-2.304902, 3.522345,
## 4.173242), .Names = c('par1',
sns <- structure(c(-1.641197, -1.152821, -0.697335, -1.751917,
    3.720378, 4.034059), .Names = c("par1", "par2", "par3", "par4",
    "par5", "par6"))
spawnNewActive <- 10^sns[1:4]
sizeCutoffs <- 10^(sns[5:6])
fireSize <- 30000

## 100
sns <- c(-0.77716149196811, -0.769325340166688, -1.2772046867758,
    -1.99332102853805, 3.14260408212431, 4.46155184064992)

## 1000
sns <- c(-0.775107, -1.03176, -0.599669, -1.958105, 3.048958,
    4.275831)

## seemed good for 100,000, pretty good for 1e3
sns <- c(-1.54885, -0.97052, -1.38305, -1.93759, 3.20379, 4.13237)

## good for 100 000, 10 000 ha -- too sinuous for 1000 and
## 100 ha
sns <- c(-1.537203, -1.462981, -0.524957, -1.002567, 3.642046,
    4.501754)

## good for 100 000, 10 000 ha (except some fires @ 1e5
## don't make it to full size) -- too sinuous for smaller
sns <- c(-1.484338, -1.22044, -2.948275, -2.15594, 3.945281,
    4.904893)

sns <- c(-1.495247, -0.800494, -1.58235, -2.270646, 3.530671,
    4.663245)

## final optimization after 75 iterations, Good: 1e5, 1e4
sns <- c(-1.47809, -0.86224, -1.34532, -1.93568, 3.27149, 4.20741)

## based on equal weights 10^(1:5)
sns <- c(-0.923528, -1.804549, -1.760455, -1.793594, 1.683355,
    4.466668)

## With spreadProb = 0.9 # Pretty GOOD!
sns <- c(-0.73152, -0.501823, -0.605968, -1.809726, 2.202732,
    4.69606, 0.9)  ## used in module

## With spreadProb = 0.9 # Optimal
sns <- c(-0.978947, -0.540946, -0.790736, -1.583039, 2.532013,
    4.267547, 0.94673)

spawnNewActive <- 10^sns[1:4]
sizeCutoffs <- 10^(sns[5:6])
if (length(sns) == 7) spreadProb <- sns[7]

# from linear model version
par <- c(1.548899, -0.396904, 2.191424, 3.903082, 4.854002)
sizeCutoffs <- 10^c(par[4], par[5])
sna <- min(-0.15, par[1] + par[2] * log10(fireSize))
sna <- 10^c(sna * par[3], sna * 2 * par[3], sna * 3 * par[3],
    sna * 4 * par[3])
spawnNewActive <- sna

###########################
clearPlot()
dev()

for (i in 1:5) {
    fireSize <- 10^i
    dim <- round(sqrt(fireSize) * 5 * 250)
    ros <- raster(extent(0, dim, 0, dim), res = 250, vals = 1)
    centreCell <- cellFromRowCol(ros, rownr = nrow(ros)/2, colnr = ncol(ros)/2)
    reps <- paste0("rep", 1:4 + (log10(fireSize) - 1) * 4)
    burnedMapList <- wrap(cl = cl, reps = reps, par = TRUE, burn1 = burn1,
        burnFun = burnFun, objs = objs, libs = addlPkgs)
    names(burnedMapList$out) <- reps
    burnedMapList <- purrr::transpose(burnedMapList$out)
    do.call(rbind, burnedMapList$LM)

    Plot(burnedMapList$burnedMap, cols = c("red"), new = FALSE,
        na.color = "white", legend = FALSE, visualSqueeze = 0.7,
        title = paste0("Fire size: ", fireSize))
}
reps <- paste0("rep", 1:1)
perims <- list()
perm <- list()
mod <- list()
dev()
clearPlot()
fireSizes <- 10^(4)

for (fs in fireSizes) {
    for (i in 1:1) {
        ros <- raster(extent(0, 2e+05, 0, 2e+05), res = 250,
            vals = 1)
        NineCorners <- cellFromRowCol(ros, rownr = nrow(ros)/4 *
            rep(1:3, 3), colnr = ncol(ros)/4 * rep(1:3, each = 3))
        centreCell <- NineCorners
        ran <- runif(4, -3, -1)
        spawnNewActive <- 10^ran
        # spawnNewActive <- 10^c(-0.1, -0.75, -1.2,
        # ran*2.5)
        fireSize = rep(fs, length(centreCell))
        sizeCutoffs <- 10^c(1, 3)
        burnedMapList <- wrap(cl = cl, reps = reps, par = TRUE,
            burn1 = burn1, burnFun = burnFun, objs = objs, libs = addlPkgs)
        names(burnedMapList$out) <- reps
        burnedMapList <- purrr::transpose(burnedMapList$out)
        do.call(rbind, burnedMapList$LM)
        Plot(burnedMapList$burnedMap, new = TRUE, zero.color = "white")
        perims[[i]] <- data.frame(perim = burnedMapList$LM$rep1$perim.area.ratio,
            spawnNewActive = mean(spawnNewActive), others = t(spawnNewActive))
    }
    fsChar <- as.character(fs)
    perm[[fsChar]] <- rbindlist(perims)
    perm[[fsChar]]$perim <- log10(perm[[fsChar]]$perim)
    perm[[fsChar]]$spawnNewActive <- perm[[fsChar]]$spawnNewActive
    Plot(perm[[fsChar]]$perim, perm[[fsChar]]$spawnNewActive,
        new = TRUE, addTo = paste0("fs", fsChar))
    mod[[fsChar]] <- lm(spawnNewActive ~ perim, data = perm[[fsChar]])
}

(predict(mod[["10"]], data.frame(perim = log10(0.003))))
(predict(mod[["100"]], data.frame(perim = log10(0.003))))
log10(predict(mod[["1000"]], data.frame(perim = log10(0.003))))

8.5.2 Current (2022) version

The original version was run using 100m pixels, despite the simulations being run using 250m pixels. This version uses 250m pixels.

## final optimization after 170 iterations (fitSN)
sns <- unname(opt_sn$optim$bestmem)

spawnNewActive <- 10^sns[1:4]
sizeCutoffs <- 10^(sns[5:6])
spreadProb <- sns[7]

## from linear model version
par <- c(1.548899, -0.396904, 2.191424, 3.903082, 4.854002)
sizeCutoffs <- 10^c(par[4], par[5])
sna <- min(-0.15, par[1] + par[2] * log10(fireSize))
sna <- 10^c(sna * par[3], sna * 2 * par[3], sna * 3 * par[3],
    sna * 4 * par[3])
spawnNewActive <- sna

clearPlot()

for (i in 1:5) {
    fireSize <- 10^i
    dim <- round(sqrt(fireSize) * 5 * 250)
    ros <- raster(extent(0, dim, 0, dim), res = 250, vals = 1)
    centreCell <- cellFromRowCol(ros, row = nrow(ros)/2, col = ncol(ros)/2)
    reps <- paste0("rep", 1:4 + (log10(fireSize) - 1) * 4)
    burnedMapList <- wrap(cl = cl, reps = reps, par = TRUE, burn1 = burn1,
        burnFun = burnFun, objs = objs, addlPkgs = addlPkgs)
    names(burnedMapList$out) <- reps
    burnedMapList <- purrr::transpose(burnedMapList$out)
    do.call(rbind, burnedMapList$LM)

    Plot(burnedMapList$burnedMap, cols = c("red"), new = FALSE,
        na.color = "white", legend = FALSE, visualSqueeze = 0.7,
        title = paste0("Fire size: ", fireSize))
}

8.5.3 Cleaning up 

parallel::stopCluster(cl)
unlink(cacheDir, recursive = TRUE)

8.5.4 Code and data availability

Code available from https://github.com/PredictiveEcology/LandMine.