Merge pull request #88 from PredictiveEcology/suz-calcRootC

calcRootC allows input pools to be split by SW/HW

Author: cboisvenue

DESCRIPTION

@@ -5,7 +5,7 @@ Description: Implementation of several components of the Carbon Budget Model of
     Canadian Forest Service (v3).
 URL:
     https://github.com/PredictiveEcology/CBMutils
-Version: 2.5.4.9002
+Version: 2.5.5.9000
 Authors@R: c(
     person("Céline",    "Boisvenue", email = "celine.boisvenue@nrcan-rncan.gc.ca", role = c("aut", "cre")),
     person("Alex M",    "Chubaty",   email = "achubaty@for-cast.ca",               role = c("aut"), comment = c(ORCID = "0000-0001-7146-8135")),

R/CBM-tools_calcRootC.R

@@ -1,14 +1,16 @@
 
 utils::globalVariables(c(
-  "..aboveGroundColumns"
+  "sw", "AG", "SoftwoodAG", "HardwoodAG", "SoftwoodAGB", "HardwoodAGB",
+  "SoftwoodRootB", "HardwoodRootB", "SoftwoodRootProp", "HardwoodRootProp"
 ))
 
 #' `calcRootC`
 #'
 #' `calcRootC` calculates the mass of carbon in roots pools from above ground pools
 #'
-#' @param aboveGroundC data.table with mass of carbon (tonnes/ha) in the Merch, Foliage and Other pools
-#' @param sw_hw a boolean vector indicating if the cohort is softwood (0) or hardwood (1)
+#' @param aboveGroundC data.table of above ground biomass with the mass of carbon (tonnes/ha) in each pool.
+#' Columns can either be `sw` (Softwood = TRUE, hardwood = FALSE), `Merch`, `Foliage` and `Other`,
+#' or `SoftwoodMerch`, `HardwoodMerch`, `SoftwoodFoliage`, `HardwoodFoliage`, `SoftwoodOther`, and `HardwoodOther`.
 #' @param a_sw "a" value for softwood root biomass
 #' @param b_sw "b" value for softwood root biomass
 #' @param a_hw "a" value for hardwood root biomass
@@ -27,47 +29,65 @@ utils::globalVariables(c(
 #' @returns data.table with mass of carbon (tonnes/ha) in coarseRoots and fineRoots pools.
 #' @export
 #'
-calcRootC <- function(aboveGroundC, sw_hw,
+calcRootC <- function(aboveGroundC,
                       a_sw = 0.222, b_sw = 1,
                       a_hw = 1.576, b_hw = 0.615,
                       a_frp = 0.072, b_frp = 0.354, c_frp = -0.060212,
                       biomassToCarbonRate = 0.5){
 
-  aboveGroundColumns <- c("Merch", "Foliage", "Other")
+  # Choose column set
+  AGcols <- list(
+    c("SoftwoodMerch", "HardwoodMerch", "SoftwoodFoliage", "HardwoodFoliage", "SoftwoodOther", "HardwoodOther"),
+    c("sw", "Merch", "Foliage", "Other")
+  )
 
-  if(!all(aboveGroundColumns %in% names(aboveGroundC))) {
-    stop("aboveGroundC needs the columns: ", paste(aboveGroundColumns, collapse = ", "))
-  }
-  aboveGroundC <- as.data.table(aboveGroundC)
+  whichCols <- which(sapply(AGcols, function(cols) all(tolower(cols) %in% tolower(names(aboveGroundC)))))
+  if (length(whichCols) == 0) stop(
+    "aboveGroundC needs one of these column sets:\n- ",
+    paste(sapply(AGcols, function(AGcol) paste(shQuote(AGcol), collapse = ", ")),
+          collapse = "\n- "))
+
+  AGcols <- AGcols[[whichCols[[1]]]]
+  aboveGroundC <- data.table::as.data.table(aboveGroundC)
+  data.table::setnames(aboveGroundC, tolower(AGcols), AGcols, skip_absent = TRUE)
+  aboveGroundC <- aboveGroundC[, .SD, .SDcols = AGcols]
 
   # Calculate the total above ground mass of carbon
-  totAGC <- rowSums(aboveGroundC[, ..aboveGroundColumns])
+  if (whichCols == 1){
 
-  # Convert Mg/ha of Carbon to Mg/ha of biomass
-  totAGB <- totAGC / biomassToCarbonRate
+    aboveGroundC[, SoftwoodAG := rowSums(aboveGroundC[, .(SoftwoodMerch, SoftwoodFoliage, SoftwoodOther)])]
+    aboveGroundC[, HardwoodAG := rowSums(aboveGroundC[, .(HardwoodMerch, HardwoodFoliage, HardwoodOther)])]
 
-  # Calculate root total biomass
-  if(!all(sw_hw %in% c(1,0))) {
-    stop("sw_hw needs to be a boolean vector")
+  }else{
+
+    if (!is.logical(aboveGroundC$sw)) stop("aboveGroundC 'sw' column must be logical")
+
+    aboveGroundC[, AG := rowSums(aboveGroundC[, .(Merch, Foliage, Other)])]
+    aboveGroundC[, SoftwoodAG := data.table::fifelse( sw, AG, 0)]
+    aboveGroundC[, HardwoodAG := data.table::fifelse(!sw, AG, 0)]
+    aboveGroundC[, AG := NULL]
   }
 
-  rootTotBiom <- ifelse(sw_hw == 0,
-                        a_sw * totAGB^b_sw,
-                        a_hw * totAGB^b_hw)
+  # Convert Mg/ha of Carbon to Mg/ha of biomass
+  aboveGroundC[, SoftwoodAGB := SoftwoodAG / biomassToCarbonRate]
+  aboveGroundC[, HardwoodAGB := HardwoodAG / biomassToCarbonRate]
+
+  # Calculate root total biomass
+  aboveGroundC[, SoftwoodRootB := a_sw * SoftwoodAGB^b_sw]
+  aboveGroundC[, HardwoodRootB := a_hw * HardwoodAGB^b_hw]
 
   # Calculate the proportion of fine roots
-  fineRootProp <- a_frp + b_frp * exp(c_frp * rootTotBiom)
+  aboveGroundC[, SoftwoodRootProp := a_frp + b_frp * exp(c_frp * SoftwoodRootB)]
+  aboveGroundC[, HardwoodRootProp := a_frp + b_frp * exp(c_frp * HardwoodRootB)]
 
+  # Calculate tonnes/ha of carbon
+  aboveGroundC[, SoftwoodCoarseRoots := biomassToCarbonRate * SoftwoodRootB * (1 - SoftwoodRootProp)]
+  aboveGroundC[, SoftwoodFineRoots   := biomassToCarbonRate * SoftwoodRootB * SoftwoodRootProp]
+  aboveGroundC[, HardwoodCoarseRoots := biomassToCarbonRate * HardwoodRootB * (1 - HardwoodRootProp)]
+  aboveGroundC[, HardwoodFineRoots   := biomassToCarbonRate * HardwoodRootB * HardwoodRootProp]
 
-  # Calculate the proportion of fine roots
-  rootBiom <- data.table(
-    coarseRoots = rootTotBiom * (1 - fineRootProp),
-    fineRoots = rootTotBiom * fineRootProp
-  )
+  return(aboveGroundC[, .(SoftwoodCoarseRoots, HardwoodCoarseRoots, SoftwoodFineRoots, HardwoodFineRoots)])
+}
 
-  # Reconvert into tonnes/ha of carbon
-  rootC <- rootBiom * biomassToCarbonRate
 
-  return(rootC)
 
-}

man/calcRootC.Rd

@@ -2,31 +2,52 @@
 if (!testthat::is_testing()) source(testthat::test_path("setup.R"))
 
 test_that("calcRootC", {
-  ABC <- data.table(Merch = c(0,1),
+
+  # Columns type 1
+  ABC <- data.table(sw      = c(TRUE, FALSE),
+                    Merch   = c(0,1),
                     Foliage = c(3,1),
-                    Other = c(1,1))
-  sw_hw = c(0,1)
-  rootC <- calcRootC(aboveGroundC = ABC, sw_hw = sw_hw)
-
-  # Works correctly
-  expect_identical(
-    round(rootC,3),
-    data.table(coarseRoots = c(0.542, 1.570),
-               fineRoots = c(0.346, 0.802))
+                    Other   = c(1,1))
+
+  rootC <- calcRootC(aboveGroundC = ABC)
+
+  expect_equal(
+    rootC[, .(CoarseRoots = SoftwoodCoarseRoots + HardwoodCoarseRoots,
+              FineRoots   = SoftwoodFineRoots   + HardwoodFineRoots)],
+    data.table(CoarseRoots = c(0.542, 1.570),
+               FineRoots   = c(0.346, 0.802)),
+    tolerance = 0.001, scale = 1
   )
-  expect_named(rootC, c("coarseRoots", "fineRoots"))
 
-  # Error with missspecified inputs
-  expect_error(
-    calcRootC(aboveGroundC = ABC[,.(Merch, Foliage, other = Other)], sw_hw = sw_hw)
+  # Columns type 2
+  ABC <- data.table(SoftwoodMerch   = c(0,0),
+                    SoftwoodFoliage = c(3,0),
+                    SoftwoodOther   = c(1,0),
+                    HardwoodMerch   = c(0,1),
+                    HardwoodFoliage = c(0,1),
+                    HardwoodOther   = c(0,1))
+
+  rootC <- calcRootC(aboveGroundC = ABC)
+
+  expect_equal(
+    rootC,
+    data.table(SoftwoodCoarseRoots = c(0.542, 0.000),
+               HardwoodCoarseRoots = c(0.000, 1.570),
+               SoftwoodFineRoots   = c(0.346, 0.000),
+               HardwoodFineRoots   = c(0.000, 0.802)),
+    tolerance = 0.001, scale = 1
   )
 
+  # Error with miss specified inputs
   expect_error(
-    calcRootC(aboveGroundC = ABC[,.(Merch, Foliage)], sw_hw = sw_hw)
+    calcRootC(aboveGroundC = ABC[,.(sw, Merch, Foliage)])
   )
 
   expect_error(
-    calcRootC(aboveGroundC = ABC, sw_hw = c("sw", "hw"))
+    calcRootC(aboveGroundC = ABC[,.(sw = c(0, 1), Merch, Foliage, Other)])
+  )
+  expect_error(
+    calcRootC(aboveGroundC = ABC[,.(sw = c("sw", "hw"), Merch, Foliage, Other)])
   )
-
 })
+