Fix climateYear shadowing, stale climate rasters, categorical resampling and rstLCC setup

fireSense_dataPrepPredict.Rmd

@@ -45,18 +45,20 @@ download.file(url = "https://img.shields.io/badge/Made%20with-Markdown-1f425f.pn
 
 ### Module summary
 
-<!-- TODO -->
-fireSense [@Marchal:2017a; @Marchal:2017b; @Marchal:2019]
+Prepares, each year, the covariate tables that the fireSense [@Marchal:2017a; @Marchal:2017b; @Marchal:2019] predict modules use:
 
-Provide a brief summary of what the module does / how to use the module.
+- `fireSense_igAndEscapePred_Covariates` for *fireSense_IgnitionPredict* and *fireSense_EscapePredict*: fuel classes, non-forest landcover, `youngAge`, ignition climate and lightning days, aggregated by `igAggFactor`.
+- `fireSense_SpreadCovariates` for *fireSense_SpreadPredict*: the same fuel, landcover and `youngAge` columns plus spread climate, at the resolution of `flammableRTM`.
 
-Module documentation should be written so that others can use your module.
-This is a template for module documentation, and should be changed to reflect your module.
+Fuel classes come from `cohortData` and `pixelGroupMap`, grouped by the `fuelClassCol` column of `sppEquiv`.
+The covariates are built by the same *fireSenseUtils* functions that *fireSense_dataPrepFit* uses, so they match the fitted models.
 
 ### Module inputs and parameters
 
-Describe input data required by the module and how to obtain it (e.g., directly from online sources or supplied by other modules)
-If `sourceURL` is specified, `downloadData("fireSense_dataPrepPredict", "..")` may be sufficient.
+`cohortData`, `pixelGroupMap` and `rstCurrentBurn` come from the vegetation and fire modules during the simulation.
+`sppEquiv`, `nonForestedLCCGroups`, `missingLCCgroup`, `lightningMaps`, `flammableRTM` and `landcoverDT` should be the ones used by *fireSense_dataPrepFit*.
+Climate is either `currentClimateRasters`, or `projectedClimateRasters` with layers named `year<year>`.
+`igAggFactor` is overwritten in `init` with the value set in the other modules.
 
 Table \@ref(tab:moduleInputs-fireSense-dataPrepPredict) shows the full list of module inputs.
 
@@ -78,18 +80,16 @@ knitr::kable(df_params, caption = "List of (ref:fireSense-dataPrepPredict) param
 
 ### Events
 
-<!-- TODO -->
-Describe what happens for each event type.
+All events after `init`, except `save`, repeat every `fireTimeStep` years.
 
-### Plotting
+- `init`: aligns `standAgeMap` and `rstLCC_RTM` to `rasterToMatch`; builds `landcoverDT` if absent; builds `nonForest_timeSinceDisturbance` if absent, from the fire polygons of the `cutoffForYoungAge` years up to `dataYear`.
+- `getClimateRasters` (from `.runInitialTime`): a supplied `currentClimateRasters` (e.g. from the `climateYear` module) is left alone. If it is absent, or this module built it for another year, takes layer `year<Y>` of each element of `projectedClimateRasters`, where `Y` is `climateYear` if supplied, else `time(sim)`. Stops if it does not match `pixelGroupMap`.
+- `prepIgAndEscPredictData` (from `.runInitialTime`): builds `fireSense_igAndEscapePred_Covariates`. Scheduled if `whichModulesToPrepare` has `fireSense_IgnitionPredict` or `fireSense_EscapePredict`.
+- `prepSpreadPredictData` (from `.runInitialTime`): builds `fireSense_SpreadCovariates`. Scheduled if `whichModulesToPrepare` has `fireSense_SpreadPredict`.
+- `ageNonForest` (from `time(sim) + 1`): adds 1 to `nonForest_timeSinceDisturbance` and resets pixels burned in `rstCurrentBurn` to 0.
+- `save`: does nothing except emit a message. The module never schedules it.
 
-<!-- TODO -->
-Write what is plotted.
-
-### Saving
-
-<!-- TODO -->
-Write what is saved.
+The module does not plot or save anything.
 
 ### Module outputs
 
@@ -103,8 +103,8 @@ knitr::kable(df_outputs, caption = "List of (ref:fireSense-dataPrepPredict) outp
 
 ### Links to other modules
 
-<!-- TODO: link to other fireSense modules -->
-Describe any anticipated linkages to other modules, such as modules that supply input data or do post-hoc analysis.
+Runs after *Biomass_borealDataPrep*, *fireSense_dataPrepFit*, *fireSense_IgnitionFit* and *fireSense_SpreadFit*, and supplies *fireSense_IgnitionPredict*, *fireSense_EscapePredict* and *fireSense_SpreadPredict*.
+It is normally run as part of the [fireSense](https://github.com/PredictiveEcology/fireSense) module group.
 
 ### Getting help
 

tests/testthat/helper-toyPredict.R

@@ -0,0 +1,89 @@
+## A tiny in-memory landscape for the event-level tests. Nothing is downloaded: every input the
+## module's `.inputObjects()` would otherwise build is supplied, so the guards there short-circuit.
+##
+## 4 x 4 pixels of 250 m. Landcover is categorical (210 forest, 20 water, 50 herb), which is the
+## point of the `method = "near"` test. Climate is one MDC stack whose `year<yyyy>` layer is
+## constant at (yyyy - 1900), so a wrong year is visible in the values.
+
+toyCRSp <- "EPSG:3978"
+
+## setup.R's `moduleName`/`testPaths` are not in scope for helpers, so resolve them here.
+toyModuleName <- "fireSense_dataPrepPredict"
+
+toyRastP <- function(vals, nam = "layer", n = 4) {
+  r <- terra::rast(nrows = n, ncols = n, xmin = 0, xmax = 1000, ymin = 0, ymax = 1000,
+                   crs = toyCRSp)
+  r <- terra::setValues(r, vals)
+  names(r) <- nam
+  r
+}
+
+toyLCCvalsP <- function() rep(c(210L, 210L, 50L, 20L), times = 4)
+
+## MDC stack: layer "year2001" is 101 everywhere, "year2002" is 102, ...
+toyClimateP <- function(years = 2001:2005) {
+  lyrs <- lapply(years, function(y) toyRastP(rep(y - 1900, 16), paste0("year", y)))
+  list(MDC = terra::rast(lyrs))
+}
+
+toyInputsP <- function(...) {
+  rtm <- toyRastP(rep(1L, 16), "rtm")
+  lcc <- toyRastP(toyLCCvalsP(), "lcc")
+  out <- list(
+    climateVariablesForFire = list(ignition = "MDC", spread = "MDC"),
+    projectedClimateRasters = toyClimateP(),
+    pixelGroupMap = toyRastP(rep(1L, 16), "pixelGroup"),
+    rasterToMatch = rtm,
+    rstLCC_RTM = lcc,
+    ## supplying rstLCCs takes the `suppliedElsewhere("rstLCCs", ...)` branch of
+    ## .inputObjects(), so nothing is downloaded
+    rstLCCs = list(year2001 = lcc),
+    flammableRTM = toyRastP(as.integer(toyLCCvalsP() != 20L), "flammable"),
+    standAgeMap = toyRastP(rep(80L, 16), "standAge"),
+    nonForest_timeSinceDisturbance = toyRastP(rep(30L, 16), "TSD"),
+    nonForestedLCCGroups = list(herb = 50L),
+    sppEquiv = data.table::data.table(LandR = "Pice_mar", FuelClass = "BlkSprc"),
+    cohortData = data.table::data.table(pixelGroup = 1L,
+                                        speciesCode = factor("Pice_mar"),
+                                        age = 80L, B = 3000L)
+  )
+  utils::modifyList(out, list(...))
+}
+
+toyParamsP <- function(...) {
+  utils::modifyList(list(whichModulesToPrepare = character(0),
+                         dataYear = 2001,
+                         .useCache = FALSE), list(...))
+}
+
+toyPathsP <- function() {
+  root <- withr::local_tempdir(.local_envir = parent.frame())
+  mp <- dirname(normalizePath(file.path("..", ".."), winslash = "/", mustWork = TRUE))
+  list(cachePath = file.path(root, "cache"), inputPath = file.path(root, "inputs"),
+       modulePath = mp, outputPath = file.path(root, "outputs"))
+}
+
+toySimInitP <- function(objects = toyInputsP(), params = toyParamsP(),
+                        start = 2001, end = 2003) {
+  withr::local_options(spades.useRequire = FALSE, spades.moduleCodeChecks = FALSE,
+                       reproducible.verbose = -2)
+  suppressMessages(SpaDES.core::simInit(
+    times = list(start = start, end = end), modules = toyModuleName,
+    params = stats::setNames(list(params), toyModuleName),
+    objects = objects, paths = toyPathsP()))
+}
+
+## this module's mod$ objects in a simList
+toyModP <- function(sim) sim[[".modObjs"]][[toyModuleName]]
+
+## The module's own functions. Under `convertToPackage()` they live in the module's namespace;
+## reaching them through the simList works either way.
+toyFunP <- function(sim, nm) get(nm, envir = sim[[".mods"]][[toyModuleName]])
+
+## Run named events of this module through spades(), which is what gives the module's functions
+## their `mod` and `P()` context. Messages are muffled.
+toyRunEventsP <- function(sim, events) {
+  withr::local_options(spades.useRequire = FALSE, reproducible.verbose = -2)
+  suppressMessages(suppressWarnings(
+    SpaDES.core::spades(sim, events = stats::setNames(list(events), toyModuleName), debug = FALSE)))
+}

tests/testthat/setup-toyLandscape.R

@@ -0,0 +1,114 @@
+## A 4 x 4 toy landscape on which every covariate can be worked out by hand.
+## A setup file rather than a helper, so that it shares an environment with `moduleName`
+## and `testPaths` from setup.R.
+##
+## Tests run inside the namespace of the package rendition; tables are handled with base
+## subsetting and data.table::set() only, so nothing depends on data.table-awareness there.
+##
+## The module calls postProcess() and Cache() unqualified but does not list `reproducible`
+## in `reqdPkgs`; in a project another module attaches it. It is always installed
+## (SpaDES.core imports it), so attach it here.
+suppressPackageStartupMessages(library(reproducible))
+
+## Cell numbers (row-major) and what is in them:
+##
+##    1  2  3  4      PG1 PG1 PG2 PG2      pixel groups 1-4 are forest with cohorts
+##    5  6  7  8      PG3 PG3 PG4 PG4
+##    9 10 11 12      wet wet grs grs      non-forest: wetland, grass
+##   13 14 15 16      --- --- wet  X       13, 14: forest land cover but no cohorts
+##                                         15: wetland burned 5 y ago; 16: not flammable
+##
+## PG1: Pice_mar 80 y B 2000 + Pinu_ban 40 y B 1000 -> class1 B = 3000
+## PG2: Popu_tre 50 y B 500                         -> class2 B = 500
+## PG3: Pice_mar 10 y B 300 (max age 10 <= cutoff 15) -> youngAge
+## PG4: Pice_mar 60 y B 100 + Popu_tre 60 y B 50    -> class1 B = 100, class2 B = 50
+toyRast <- function(vals) {
+  terra::rast(nrows = 4, ncols = 4, xmin = 0, xmax = 4, ymin = 0, ymax = 4, vals = vals,
+              crs = "EPSG:3005")
+}
+
+toyLCC <- function() toyRast(c(rep(210, 8), 19, 19, 16, 16, 210, 210, 19, 20))
+
+toyCohortData <- function() {
+  data.table::data.table(
+    pixelGroup  = c(1L, 1L, 2L, 3L, 4L, 4L),
+    speciesCode = factor(c("Pice_mar", "Pinu_ban", "Popu_tre", "Pice_mar", "Pice_mar", "Popu_tre")),
+    age         = c(80L, 40L, 50L, 10L, 60L, 60L),
+    B           = c(2000L, 1000L, 500L, 300L, 100L, 50L)
+  )
+}
+
+toySppEquiv <- function() {
+  data.table::data.table(LandR = c("Pice_mar", "Pinu_ban", "Popu_tre"),
+                         FuelClass = c("class1", "class1", "class2"))
+}
+
+toyLandcoverDT <- function() {
+  data.table::data.table(pixelID = 1:15,
+                         wetland = as.integer(1:15 %in% c(9, 10, 15)),
+                         grass   = as.integer(1:15 %in% c(11, 12)))
+}
+
+## climate: layer `year<Y>` of MDC is (Y - 2000) * 100 + cell number, so both the year that
+## was taken and the cell it came from can be read off any value; `Tmax` is its negative
+toyClimate <- function(years = 2001:2003) {
+  mk <- function(sign) {
+    r <- terra::rast(lapply(years, function(y) toyRast(sign * ((y - 2000) * 100 + 1:16))))
+    names(r) <- paste0("year", years)
+    r
+  }
+  list(MDC = mk(1), Tmax = mk(-1))
+}
+
+toyObjects <- function() {
+  list(
+    rasterToMatch = toyRast(1),
+    studyArea = terra::as.polygons(terra::ext(toyRast(1)), crs = "EPSG:3005"),
+    flammableRTM = toyRast(c(rep(1, 15), 0)),
+    pixelGroupMap = toyRast(c(1, 1, 2, 2, 3, 3, 4, 4, rep(NA, 8))),
+    cohortData = toyCohortData(),
+    sppEquiv = toySppEquiv(),
+    landcoverDT = toyLandcoverDT(),
+    nonForestedLCCGroups = list(wetland = 19L, grass = 16L),
+    missingLCCgroup = "grass",
+    ## years since fire: 30 everywhere except the wetland cell 15 and the forest cell 1
+    nonForest_timeSinceDisturbance = toyRast(c(5, rep(30, 13), 5, 30)),
+    standAgeMap = toyRast(50),
+    ## supplying `rstLCCs` takes the branch of `.inputObjects` that does not download NTEMS
+    ## landcover; `rstLCC_RTM` is then its last element
+    rstLCCs = list(toyLCC()),
+    projectedClimateRasters = toyClimate(),
+    climateVariablesForFire = list(ignition = "MDC", spread = "MDC"),
+    lightningMaps = {
+      r <- c(toyRast(1:16), toyRast(1001:1016))
+      names(r) <- c("lightningDays", "lightningDensity")
+      r
+    }
+  )
+}
+
+toyPrepSim <- function(objects = toyObjects(), params = list(),
+                       times = list(start = 2001, end = 2001)) {
+  params <- utils::modifyList(list(igAggFactor = 2), params)
+  sim <- SpaDES.core::simInit(
+    times = c(times, timeunit = "year"),
+    modules = moduleName,
+    params = stats::setNames(list(params), moduleName),
+    objects = objects,
+    paths = testPaths
+  )
+  sim
+}
+
+toyPrepRun <- function(...) SpaDES.core::spades(toyPrepSim(...), debug = FALSE)
+
+## a covariate table as a data.frame ordered by pixelID
+covDF <- function(dt) {
+  df <- as.data.frame(dt)
+  df[order(df$pixelID), , drop = FALSE]
+}
+
+evOf <- function(dt, type) {
+  df <- as.data.frame(dt)
+  df[df$moduleName == "fireSense_dataPrepPredict" & df$eventType == type, , drop = FALSE]
+}

tests/testthat/test-ageNonForest.R

@@ -0,0 +1,39 @@
+## `ageNonForest()` is the one helper in this module with no dependencies at all: it takes
+## two rasters and returns one. Every value below is worked out by hand from the toy TSD.
+
+test_that("ageNonForest adds one year everywhere when nothing burned", {
+  TSD <- toyRast(c(5, rep(30, 13), 5, 30))
+  out <- ageNonForest(TSD = TSD, rstCurrentBurn = NULL, timeStep = 1)
+  expect_s4_class(out, "SpatRaster")
+  ## hand-computed: every cell + 1
+  expect_identical(as.vector(terra::values(out)), c(6, rep(31, 13), 6, 31))
+  ## geometry must survive: the values are written back with setValues()
+  expect_true(terra::compareGeom(out, TSD, stopOnError = FALSE))
+})
+
+test_that("ageNonForest resets burned pixels to zero and ages the rest", {
+  TSD <- toyRast(c(5, rep(30, 13), 5, 30))
+  ## cell 2 burned (1); cell 3 is NA (no burn data, i.e. did not burn); cell 16 burned
+  burn <- toyRast(c(0, 1, NA, rep(0, 12), 1))
+  out <- ageNonForest(TSD = TSD, rstCurrentBurn = burn, timeStep = 1)
+  ## hand-computed: 5+1; burned -> 0; NA is unburned so 30+1; cell 15 (TSD 5, unburned) -> 6;
+  ## cell 16 burned -> 0
+  expect_identical(as.vector(terra::values(out)), c(6, 0, rep(31, 12), 6, 0))
+})
+
+test_that("ageNonForest treats NA and 0 in rstCurrentBurn identically", {
+  TSD <- toyRast(rep(10, 16))
+  allNA <- toyRast(rep(NA_real_, 16))
+  allZero <- toyRast(rep(0, 16))
+  expect_identical(as.vector(terra::values(ageNonForest(TSD, allNA, 1))), rep(11, 16))
+  expect_identical(as.vector(terra::values(ageNonForest(TSD, allNA, 1))),
+                   as.vector(terra::values(ageNonForest(TSD, allZero, 1))))
+})
+
+test_that("ageNonForest ignores timeStep, as documented", {
+  ## `timeStep` is accepted but the raster always advances by exactly one year. Asserted so
+  ## that honouring it becomes a deliberate, visible change rather than a silent one.
+  TSD <- toyRast(rep(10, 16))
+  expect_identical(as.vector(terra::values(ageNonForest(TSD, NULL, timeStep = 1))),
+                   as.vector(terra::values(ageNonForest(TSD, NULL, timeStep = 10))))
+})