Make every fitted ELF's fuel covariates for every pixel

NEWS.md

@@ -1,3 +1,7 @@
+# fireSense_dataPrepPredict (development version)
+
+- Several fitted ELFs in one study area: with `sppEquivs`, `nonForestedLCCGroupsList` and `missingLCCgroupList` from `fireSense_dataPrepFit` (one element per ELF), every ELF's fuel covariates are made for every pixel, each with its own species table, non-forest groups and `landcoverDT`, and merged by `pixelID` (spread) or by layer (ignition). `fireSense_SpreadPredict` and `fireSense_IgnitionPredict` then apply each ELF's model to its own columns. One ELF works as before.
+
 # fireSense_dataPrepPredict 1.0.2
 
 First release from `development` since `main` was last updated (2022-02-28). Full history: https://github.com/PredictiveEcology/fireSense_dataPrepPredict/compare/9b276ed...v1.0.2

fireSense_dataPrepPredict.R

@@ -10,7 +10,7 @@ defineModule(sim, list(
     person("Alex M", "Chubaty", role = "ctb", email = "achubaty@for-cast.ca")
   ),
   childModules = character(0),
-  version = list(fireSense_dataPrepPredict = "1.0.4.9000"),
+  version = list(fireSense_dataPrepPredict = "1.0.4.9001"),
   timeframe = as.POSIXlt(c(NA, NA)),
   timeunit = "year",
   citation = list("citation.bib"),
@@ -120,6 +120,13 @@ defineModule(sim, list(
       desc = paste(
         "Optional named list of landcover `SpatRaster`s, one per data year, as produced by",
         "`fireSense_dataPrepFit`. If supplied, the last element is used as `rstLCC_RTM`.")),
+    expectsInput("sppEquivs", "list", sourceURL = NA,
+      desc = paste("Only with several fitted ELFs: one `sppEquiv` per ELF, from `fireSense_dataPrepFit`. Each ELF's",
+                   "fuel covariates are made with its own, for every pixel.")),
+    expectsInput("nonForestedLCCGroupsList", "list", sourceURL = NA,
+      desc = "Only with several fitted ELFs: one `nonForestedLCCGroups` per ELF, in the order of `sppEquivs`."),
+    expectsInput("missingLCCgroupList", "list", sourceURL = NA,
+      desc = "Only with several fitted ELFs: one `missingLCCgroup` per ELF, in the order of `sppEquivs`."),
     expectsInput("sppEquiv", "data.table", sourceURL = NA,
       desc = "Table of LandR species equivalencies; must have columns `sppEquivCol` and `fuelClassCol`."),
     expectsInput("standAgeMap", "SpatRaster", sourceURL = NA,
@@ -378,28 +385,27 @@ prepare_IgnitionAndEscapePredict <- function(sim) {
 
   # Coming out of the CacheGeo, this is unreliably a data.frame instead of a data.table
   if (!data.table::is.data.table(sim$sppEquiv)) data.table::setDT(sim$sppEquiv)
-  if (is.null(mod$requiredFuelClasses))
-    mod$requiredFuelClasses <- sim$sppEquiv[[P(sim)$fuelClassCol]]
-  # if ignition and spread fuel classes are the same, this should use Mod
-  # to avoid doing it twice (in spreadFit, assuming people run both events)
-  fuelCovsCoarse <- prepare_FuelCovsCoarse(
+  ## one fuel set per fitted ELF (one, as before, when there is one ELF); the covariate tables are merged,
+  ## each ELF's columns alongside the others', for fireSense_IgnitionPredict to pick its own
+  fuelSets <- ELFfuelSets(sim)
+  fuelCovsCoarse <- mergeCovariateTables(lapply(fuelSets, function(fs) prepare_FuelCovsCoarse(
     cohortData = sim$cohortData,
     pixelGroupMap = sim$pixelGroupMap,
     flammableRTM = sim$flammableRTM,
-    landcoverDT = sim$landcoverDT,
+    landcoverDT = fs$landcoverDT,
     nonForest_timeSinceDisturbance = sim$nonForest_timeSinceDisturbance,
-    sppEquiv = sim$sppEquiv,
+    sppEquiv = fs$sppEquiv,
     sppEquivCol = P(sim)$sppEquivCol,
     fuelClassCol = P(sim)$fuelClassCol,
-    requiredFuelClasses = mod$requiredFuelClasses,
+    requiredFuelClasses = fs$requiredFuelClasses,
     cutoffForYoungAge = P(sim)$cutoffForYoungAge,
-    missingLCCgroup = sim$missingLCCgroup,
-    nonForestedLCCGroups = sim$nonForestedLCCGroups,
+    missingLCCgroup = fs$missingLCCgroup,
+    nonForestedLCCGroups = fs$nonForestedLCCGroups,
     nonForestCanBeYoungAge = P(sim)$nonForestCanBeYoungAge,
     studyAreaName = P(sim)$.studyAreaName,
     rasTemplate = sim$flammableRTM, fact = Par$igAggFactor,
     useCache = FALSE
-  )
+  )))
 
   ignitionClimateCoarse <- prepare_ignitionClimate(
     ignitionClimateList = as.list(ignitionClimate) |> setNames(names(ignitionClimate)), 
@@ -412,7 +418,7 @@ prepare_IgnitionAndEscapePredict <- function(sim) {
     mergePreparedCovs(years = list(yrLab) |> setNames(time(sim)), 
                       list(fuelCovsCoarse) |> setNames(yrLab), 
                       ignitionFirePoints = NULL, 
-                      sim$nonForestedLCCGroups,
+                      unionLCCGroups(fuelSets),
                       ignitionClimateCoarseList, 
                       sim$lightningMaps["lightningDays"], 
                       digest = append(NULL, list(P(sim)$igAggFactor)),
@@ -435,32 +441,47 @@ prepare_SpreadPredict <- function(sim) {
   if (is.null(spreadClimate))
     stop("spreadClimate is NULL; there is a problem to debug")
 
-  if (is.null(mod$requiredFuelClasses))
-    mod$requiredFuelClasses <- sim$sppEquiv[[P(sim)$fuelClassCol]]
-
-  ## much of this chunk can now be combined into a function, called for both ig and spread prep
-  ## this fits cohortData into fuel classes
-  ##  if pixels are missing/absent but are able to be forested as determined by landcoverDT,
-  ##  they receive 0 values - e.g. pixelGroup zero
-  spreadCovariates <- fireSenseUtils:::fireSenseCovariatesCreate(
-    cohortData = sim$cohortData,
-    pixelGroupMap = sim$pixelGroupMap,
-    flammableRTM = sim$flammableRTM,
-    landcoverDT = sim$landcoverDT,
-
-    sppEquiv = sim$sppEquiv,
-    sppEquivCol = P(sim)$sppEquivCol,
-    fuelClassCol = P(sim)$fuelClassCol,
-    requiredFuelClasses = mod$requiredFuelClasses,
-    cutoffForYoungAge = P(sim)$cutoffForYoungAge,
-    missingLCCgroup = sim$missingLCCgroup,
-    nonForestedLCCGroups = sim$nonForestedLCCGroups,
-    nonForestCanBeYoungAge = P(sim)$nonForestCanBeYoungAge,
-    nonForest_timeSinceDisturbance = sim$nonForest_timeSinceDisturbance,
-    studyAreaName = P(sim)$.studyAreaName,
-    useCache = FALSE # predict is annual, no point in caching
-  ) 
-
+  ## one fuel set per fitted ELF (one, as before, when there is one ELF). Every ELF's covariates are made for
+  ## every pixel and the tables merged, so fireSense_SpreadPredict can apply each ELF's model wherever it
+  ## predicts, including the blend zone around its own pixels. Column names say what they hold (fuel class,
+  ## non-forest LCC codes), so a column two ELFs share means the same thing in both.
+  fuelSets <- ELFfuelSets(sim)
+  spreadCovariates <- mergeCovariateTables(lapply(fuelSets, function(fs) {
+    ## this fits cohortData into fuel classes
+    ##  if pixels are missing/absent but are able to be forested as determined by landcoverDT,
+    ##  they receive 0 values - e.g. pixelGroup zero
+    covs <- fireSenseUtils:::fireSenseCovariatesCreate(
+      cohortData = sim$cohortData,
+      pixelGroupMap = sim$pixelGroupMap,
+      flammableRTM = sim$flammableRTM,
+      landcoverDT = fs$landcoverDT,
+
+      sppEquiv = fs$sppEquiv,
+      sppEquivCol = P(sim)$sppEquivCol,
+      fuelClassCol = P(sim)$fuelClassCol,
+      requiredFuelClasses = fs$requiredFuelClasses,
+      cutoffForYoungAge = P(sim)$cutoffForYoungAge,
+      missingLCCgroup = fs$missingLCCgroup,
+      nonForestedLCCGroups = fs$nonForestedLCCGroups,
+      nonForestCanBeYoungAge = P(sim)$nonForestCanBeYoungAge,
+      nonForest_timeSinceDisturbance = sim$nonForest_timeSinceDisturbance,
+      studyAreaName = P(sim)$.studyAreaName,
+      useCache = FALSE # predict is annual, no point in caching
+    )
+    # Sanity check - make sure the nonForest pixels have no forest fuels
+    nfCols <- setdiff(names(fs$landcoverDT), "pixelID")
+    df <- copy(covs)
+    fcs <- setdiff(colnames(df), c("pixelID", "youngAge", nfCols))
+    if (length(fcs)) {
+      df[df - min(df[[fcs[[1]]]], na.rm = TRUE) == 0] <- NA
+      fuelInSamePixelAsNonForest <- any(rowSums(!is.na(df[, ..fcs])) > 0 &
+                                          (rowSums(df[, ..nfCols]) > 0))
+      if (fuelInSamePixelAsNonForest)
+        stop("Flammable fuels (i.e., trees) are present in pixels that are identified as non-fore")
+    }
+    covs
+  }))
+
   # TODO: this chunk is untested 18/12/2024
   climateCovariates <- spreadClimate |> as.data.frame(cells = TRUE)
   climateCovariates <- na.omit(climateCovariates) |> as.data.table()
@@ -471,22 +492,79 @@ prepare_SpreadPredict <- function(sim) {
   nonFuelNames <- c("pixelID", names(spreadClimate), "youngAge")
   setcolorder(spreadCovariates, neworder = nonFuelNames)
   sim$fireSense_SpreadCovariates <- spreadCovariates
-
-  # Sanity check - make sure the nonForest pixels have no forest fuels
-  nfCols <- setdiff(names(sim$landcoverDT), "pixelID")
-  df <- sim$fireSense_SpreadCovariates
-  fcs <- setdiff(colnames(df), c(nonFuelNames, nfCols))
-  df[df - min(df[[fcs[[1]]]], na.rm = TRUE) == 0] <- NA
-  fuelInSamePixelAsNonForest <- any(rowSums(!is.na(df[, ..fcs])) > 0 & 
-                                      (rowSums(df[, ..nfCols]) > 0))
-
-  if (fuelInSamePixelAsNonForest)
-    stop("Flammable fuels (i.e., trees) are present in pixels that are identified as non-fore")
 
   return(invisible(sim))
 }
 
 
+#' The fuel sets to build covariates with: one per fitted ELF
+#'
+#' With several fitted ELFs, `fireSense_dataPrepFit` supplies one species table, non-forest grouping and
+#' missing-LCC group per ELF (`sppEquivs`, `nonForestedLCCGroupsList`, `missingLCCgroupList`); each gets its
+#' own `landcoverDT`, made once and kept in `mod`. With one ELF, the single set of objects, as before.
+#'
+#' @param sim A `simList`.
+#' @return list of lists, each with `sppEquiv`, `nonForestedLCCGroups`, `missingLCCgroup`, `landcoverDT` and
+#'   `requiredFuelClasses`.
+ELFfuelSets <- function(sim) {
+  fcc <- P(sim)$fuelClassCol
+  if (length(sim$sppEquivs) > 1L) {
+    n <- length(sim$sppEquivs)
+    if (length(sim$nonForestedLCCGroupsList) != n || length(sim$missingLCCgroupList) != n)
+      stop("fireSense_dataPrepPredict: sppEquivs, nonForestedLCCGroupsList and missingLCCgroupList must have one ",
+           "element per ELF")
+    if (is.null(mod$ELFlandcoverDTs)) {
+      rstLCC <- if (!LandR::.compareRas(sim$flammableRTM, sim$rstLCC_RTM, stopOnError = FALSE))
+        reproducible::postProcess(sim$rstLCC_RTM, to = sim$flammableRTM, method = "near") else sim$rstLCC_RTM
+      mod$ELFlandcoverDTs <- lapply(sim$nonForestedLCCGroupsList, function(g)
+        makeLandcoverDT(rstLCC = rstLCC, flammableRTM = sim$flammableRTM,
+                        forestedLCC = P(sim)$forestedLCC, nonForestedLCCGroups = g))
+    }
+    return(lapply(seq_len(n), function(i) {
+      se <- data.table::as.data.table(sim$sppEquivs[[i]])
+      list(sppEquiv = se, nonForestedLCCGroups = sim$nonForestedLCCGroupsList[[i]],
+           missingLCCgroup = sim$missingLCCgroupList[[i]], landcoverDT = mod$ELFlandcoverDTs[[i]],
+           requiredFuelClasses = se[[fcc]])
+    }))
+  }
+  if (!data.table::is.data.table(sim$sppEquiv)) data.table::setDT(sim$sppEquiv)
+  if (is.null(mod$requiredFuelClasses))
+    mod$requiredFuelClasses <- sim$sppEquiv[[fcc]]
+  list(list(sppEquiv = sim$sppEquiv, nonForestedLCCGroups = sim$nonForestedLCCGroups,
+            missingLCCgroup = sim$missingLCCgroup, landcoverDT = sim$landcoverDT,
+            requiredFuelClasses = mod$requiredFuelClasses))
+}
+
+#' Merge the covariates of several ELFs
+#'
+#' One ELF: its covariates, untouched. Several: tables are merged by `pixelID` keeping every pixel, and
+#' rasters by adding the layers the first lacks; a column or layer the ELFs share (same name, so the same
+#' content) is taken from the first.
+#'
+#' @param tabs list, one element per ELF: `data.table`s with `pixelID`, or `SpatRaster`s.
+#' @return one `data.table` or `SpatRaster`.
+mergeCovariateTables <- function(tabs) {
+  if (length(tabs) == 1L) return(tabs[[1]])
+  if (inherits(tabs[[1]], "SpatRaster"))
+    return(Reduce(function(a, b) {
+      extra <- setdiff(names(b), names(a))
+      if (length(extra)) c(a, b[[extra]]) else a
+    }, tabs))
+  tabs <- lapply(tabs, data.table::as.data.table)
+  Reduce(function(a, b) {
+    extra <- setdiff(names(b), names(a))
+    if (!length(extra)) return(a)
+    merge(a, b[, c("pixelID", extra), with = FALSE], by = "pixelID", all = TRUE)
+  }, tabs)
+}
+
+## the ELFs' non-forest groups together (a group two ELFs share has the same name and codes); one ELF: its own
+unionLCCGroups <- function(fuelSets) {
+  if (length(fuelSets) == 1L) return(fuelSets[[1]]$nonForestedLCCGroups)
+  g <- do.call(c, lapply(fuelSets, `[[`, "nonForestedLCCGroups"))
+  g[!duplicated(names(g))]
+}
+
 #' Supply default inputs
 #'
 #' Defaults for `climateVariablesForFire`, `rstLCC_RTM`, `standAgeMap`,

tests/testthat/test-metadata.R

@@ -25,14 +25,17 @@ test_that("inputs are the expected names and classes", {
       landcoverDT              = "data.table",
       lightningMaps            = "SpatRaster",
       missingLCCgroup          = "character",
+      missingLCCgroupList      = "list",
       nonForestedLCCGroups     = "list",
+      nonForestedLCCGroupsList = "list",
       pixelGroupMap            = "SpatRaster",
       projectedClimateRasters  = "list",
       rasterToMatch            = "SpatRaster",
       rstCurrentBurn           = "SpatRaster",
       rstLCC_RTM               = "SpatRaster",
       rstLCCs                  = "list",
       sppEquiv                 = "data.table",
+      sppEquivs                = "list",
       standAgeMap              = "SpatRaster")
   )
 })

tests/testthat/test-multiELF.R

@@ -0,0 +1,57 @@
+## Several fitted ELFs in one study area (the 2-ELF Mackenzie forecast, 2026-09). Each ELF's model predicts
+## over its own pixels and a blend zone around them (fireSense_SpreadPredict), so every pixel needs every
+## ELF's covariates: each ELF's fuel classes and non-forest groups, made from its own sppEquiv. Columns two
+## ELFs share have the same name and so the same content.
+##
+## ELF A groups spruce and pine as "conifer", aspen as "decid"; non-forest groups wetland and grass.
+## ELF B keeps spruce and pine apart, aspen as "decid" (shared with A); one non-forest group "wetgrs".
+## Expectation, without hand numbers: the merged table holds the union of columns, and each column equals
+## what a one-ELF run with that ELF's objects makes.
+
+sppA <- function() data.table::data.table(LandR = c("Pice_mar", "Pinu_ban", "Popu_tre"),
+                                          FuelClass = c("conifer", "conifer", "decid"))
+sppB <- function() data.table::data.table(LandR = c("Pice_mar", "Pinu_ban", "Popu_tre"),
+                                          FuelClass = c("Pice_mar", "Pinu_ban", "decid"))
+nfA <- list(wetland = 19L, grass = 16L); nfB <- list(wetgrs = c(16L, 19L))
+
+oneELF <- function(spp, nf, missing) {
+  o <- toyObjects()
+  o$sppEquiv <- spp(); o$nonForestedLCCGroups <- nf; o$missingLCCgroup <- missing
+  o$landcoverDT <- NULL                          # made from this ELF's groups, as for each ELF below
+  toyPrepRun(o)
+}
+
+test_that("with two ELFs, the covariates are both ELFs' own, side by side", {
+  o <- toyObjects()
+  o$landcoverDT <- NULL
+  o$sppEquivs <- list(sppA(), sppB())
+  o$nonForestedLCCGroupsList <- list(nfA, nfB)
+  o$missingLCCgroupList <- list("grass", "wetgrs")
+  both <- toyPrepRun(o)
+  a <- oneELF(sppA, nfA, "grass"); b <- oneELF(sppB, nfB, "wetgrs")
+
+  sp <- covDF(both$fireSense_SpreadCovariates)
+  spA <- covDF(a$fireSense_SpreadCovariates); spB <- covDF(b$fireSense_SpreadCovariates)
+  expect_setequal(names(sp), union(names(spA), names(spB)))
+  expect_identical(sp$pixelID, spA$pixelID)
+  for (cn in setdiff(names(spA), "pixelID")) expect_equal(sp[[cn]], spA[[cn]], info = cn)
+  for (cn in setdiff(names(spB), "pixelID")) expect_equal(sp[[cn]], spB[[cn]], info = cn)
+  expect_true(all(c("conifer", "Pice_mar", "Pinu_ban", "decid", "wetland", "grass", "wetgrs") %in% names(sp)))
+
+  ## ignition covariates: the same union, each ELF's layers as in its own run
+  ig <- as.data.frame(both$fireSense_igAndEscapePred_Covariates)
+  igA <- as.data.frame(a$fireSense_igAndEscapePred_Covariates)
+  igB <- as.data.frame(b$fireSense_igAndEscapePred_Covariates)
+  expect_true(all(union(names(igA), names(igB)) %in% names(ig)))
+  key <- intersect(c("pixelID", "year"), names(ig))
+  ord <- function(d) d[do.call(order, d[key]), , drop = FALSE]
+  ig <- ord(ig); igA <- ord(igA); igB <- ord(igB)
+  for (cn in setdiff(names(igA), key)) expect_equal(ig[[cn]], igA[[cn]], info = cn)
+  for (cn in setdiff(names(igB), key)) expect_equal(ig[[cn]], igB[[cn]], info = cn)
+})
+
+test_that("mismatched per-ELF lists stop with a message", {
+  o <- toyObjects()
+  o$sppEquivs <- list(sppA(), sppB()); o$nonForestedLCCGroupsList <- list(nfA); o$missingLCCgroupList <- list("grass", "wetgrs")
+  expect_error(toyPrepRun(o), "one element per ELF")
+})