Merge pull request #1739 from infotroph/biocro-1738

Fix Biocro result summaries (#1738)

CHANGELOG.md

@@ -10,6 +10,7 @@ For more information about this file see also [Keep a Changelog](http://keepacha
 
 ### Fixes
 - `PEcAn.utils` now lazy-loads data for faster execution of functions that consult lookup tables, especially `to_ncvar`.
+- Fixed incorrect `PEcAn.BIOCRO` daily and yearly results: Was calculating every row from whole simulation instead of that day (#1738)
 
 ### Added
 

models/biocro/DESCRIPTION

@@ -20,6 +20,7 @@ Imports:
     ncdf4 (>= 1.15),
     lubridate (>= 1.6.0),
     data.table,
+    dplyr,
     XML
 Suggests:
     BioCro,

models/biocro/R/run.biocro.R

@@ -159,35 +159,58 @@ run.biocro <- function(lat, lon, metpath, soil.nc = NULL, config = config, coppi
 
   hourly.results <- do.call("rbind", hourly.results)
   hourly.results <- hourly.results[order(hourly.results$year, hourly.results$doy, hourly.results$hour),]
+
+  # Compute daily and yearly results by taking max or sum as appropriate.
+  # This notation could be more compact if we used nonstandard evaluation
+  # with bare variable names, but this way works and ensures that
+  # `R CMD check` doesn't complain about undefined variables.
+  hourly_grp <- dplyr::group_by_at(.tbl = hourly.results, .vars= c("year", "doy"))
+  daily.results <- dplyr::bind_cols(
+    dplyr::summarize_at(
+      .tbl = hourly_grp,
+      .vars = c("Stem", "Leaf", "Root", "AboveLitter", "BelowLitter",
+                "Rhizome", "Grain", "LAI", tmax="Temp"),
+      .fun = max),
+    dplyr::summarize_at(
+      .tbl = hourly_grp,
+      .vars = c("SoilEvaporation", "CanopyTrans", "precip"),
+      .fun = sum),
+    dplyr::summarize_at(
+      .tbl = hourly_grp,
+      .vars = c(tmin = "Temp"),
+      .fun = min),
+    dplyr::summarize_at(
+      .tbl = hourly_grp,
+      .vars = c(tavg = "Temp"),
+      .fun = mean))
+  # bind_cols on 4 tables leaves 3 sets of duplicate year and day columns.
+  # Let's drop these.
+  col_order <- c("year", "doy", "Stem", "Leaf", "Root",
+                 "AboveLitter", "BelowLitter", "Rhizome",
+                 "SoilEvaporation", "CanopyTrans", "Grain", "LAI",
+                 "tmax", "tmin", "tavg", "precip")
+  daily.results <- daily.results[, col_order]
 
-  daily.results <- hourly.results[, list(Stem = max(hourly.results$Stem),
-                                         Leaf = max(hourly.results$Leaf),
-                                         Root = max(hourly.results$Root),
-                                         AboveLitter = max(hourly.results$AboveLitter),
-                                         BelowLitter = max(hourly.results$BelowLitter),
-                                         Rhizome = max(hourly.results$Rhizome),
-                                         SoilEvaporation = sum(hourly.results$SoilEvaporation),
-                                         CanopyTrans = sum(hourly.results$CanopyTrans),
-                                         Grain = max(hourly.results$Grain),
-                                         LAI = max(hourly.results$LAI),
-                                         tmax = max(hourly.results$Temp),
-                                         tmin = min(hourly.results$Temp),
-                                         tavg = mean(hourly.results$Temp),
-                                         precip = sum(hourly.results$precip)),
-                                  by = "year,doy"]
-
-  annual.results <- hourly.results[, list(Stem = max(hourly.results$Stem),
-                                          Leaf = max(hourly.results$Leaf),
-                                          Root = max(hourly.results$Root),
-                                          AboveLitter = max(hourly.results$AboveLitter),
-                                          BelowLitter = max(hourly.results$BelowLitter),
-                                          Rhizome = max(hourly.results$Rhizome),
-                                          Grain = max(hourly.results$Grain),
-                                          SoilEvaporation = sum(hourly.results$SoilEvaporation),
-                                          CanopyTrans = sum(hourly.results$CanopyTrans),
-                                          map = sum(hourly.results$precip),
-                                          mat = mean(hourly.results$Temp)),
-                                   by = "year"]
+  daily_grp <- dplyr::group_by_at(.tbl = hourly.results, .vars = "year")
+  annual.results <- dplyr::bind_cols(
+      dplyr::summarize_at(
+        .tbl = daily_grp,
+        .vars = c("Stem", "Leaf", "Root", "AboveLitter", "BelowLitter",
+                  "Rhizome", "Grain"),
+        .fun = max),
+      dplyr::summarize_at(
+        .tbl = daily_grp,
+        .vars = c("SoilEvaporation", "CanopyTrans", map="precip"),
+        .fun = sum),
+      dplyr::summarize_at(
+        .tbl = daily_grp,
+        .vars = c(mat = "Temp"),
+        .fun = mean))
+  col_order <- c("year", "Stem", "Leaf", "Root", "AboveLitter", "BelowLitter",
+                 "Rhizome", "Grain", "SoilEvaporation", "CanopyTrans",
+                 "map", "mat")
+  annual.results <- annual.results[, col_order]
+
   return(list(hourly = hourly.results, 
               daily = daily.results, 
               annually = data.table::data.table(lat = lat, lon = lon, annual.results)))

models/biocro/tests/testthat/test-run.biocro.R

@@ -0,0 +1,42 @@
+context("check that BioCro output is summarized correctly")
+
+# Return precalculated BioCro 0.9 results from specified days in 2004
+# Accepts same arguments as BioCro::BioGro, ignores all but day1 and dayn
+mock_run <- function(WetDat = NULL, day1 = 1, dayn = 7, ...){
+	load("data/result.RData", envir = environment())
+	resultDT[resultDT$Year == 2004 & resultDT$DayofYear >= day1 & resultDT$DayofYear <= dayn,]
+}
+
+# Hand-calculate reference values
+ref_output <- mock_run()
+ref_met <- read.csv("data/US-Bo1.2004.csv", nrows=7*24)
+ref_leaf1 <- max(ref_output$Leaf[ref_output$DayofYear == 1])
+ref_soil5 <- sum(ref_output$SoilEvaporation[ref_output$DayofYear == 5])
+ref_mat <- mean(ref_met$Temp)
+
+# run setup
+metpath <- "data/US-Bo1"
+config <- PEcAn.settings::prepare.settings(PEcAn.settings::read.settings("data/pecan.biocro.xml"))
+config$pft$type$genus <- "Salix" # TODO should PEcAn.BIOCRO expect config$pfts$pft$type$genus instead?
+config$run$start.date <- as.POSIXct("2004-01-01")
+config$run$end.date <- as.POSIXct("2004-01-07")
+config$simulationPeriod$dateofplanting <- as.POSIXct("2004-01-01")
+config$simulationPeriod$dateofharvest <- as.POSIXct("2004-01-07")
+
+test_that("daily summarizes hourly (#1738)", {
+	skip_on_travis() # stubbing only works when BioCro package is installed
+
+	# stub out BioCro::willowGro: 
+	# calls to willowGro(...) will be replaced with calls to mock_run(...),
+	# but *only* when originating inside run.biocro.
+	mockery::stub(run.biocro, "BioCro::willowGro", mock_run) 
+
+	mock_result <- run.biocro(lat = 44, lon = -88, metpath, soil.nc = NULL, config = config, coppice.interval = 1)
+	expect_equal(nrow(mock_result$hourly), 24*7)
+	expect_equal(nrow(mock_result$daily), 7)
+	expect_equal(nrow(mock_result$annually), 1)
+	expect_gt(length(unique(mock_result$daily$tmax)), 1)
+	expect_equal(mock_result$daily$Leaf[mock_result$daily$doy == 1], ref_leaf1)
+	expect_equal(mock_result$daily$SoilEvaporation[mock_result$daily$doy == 5], ref_soil5)
+	expect_equal(mock_result$annually$mat, ref_mat)
+})