feat: implement automated reporting and visualization layer

modules/benchmark/DESCRIPTION

@@ -46,6 +46,7 @@ Imports:
     zoo,
     yaml
 Suggests:
+    quarto,
     PEcAn.data.land,
     testthat (>= 2.0.0)
 License: BSD_3_clause + file LICENSE

modules/benchmark/NAMESPACE

@@ -13,6 +13,7 @@ export(clean_settings_BRR)
 export(create_BRR)
 export(define_benchmark)
 export(format_wide2long)
+export(generate_validation_report)
 export(load_and_map_data)
 export(load_csv)
 export(load_data)
@@ -42,10 +43,6 @@ export(read_settings_BRR)
 export(register_metric)
 export(run_benchmark)
 importFrom(dplyr,rename)
-importFrom(ggplot2,geom_path)
-importFrom(ggplot2,geom_point)
-importFrom(ggplot2,ggplot)
-importFrom(ggplot2,labs)
 importFrom(magrittr,"%>%")
 importFrom(rlang,.data)
 importFrom(yaml,read_yaml)

modules/benchmark/R/generate_validation_report.R

@@ -0,0 +1,67 @@
+##' Generate Validation Benchmark Report
+##'
+##' @param benchmark_results A list containing `metrics` (data.frame), `aligned_data` (data.frame), and `plots` (list of ggplot objects) returned by the validation pipeline.
+##' @param output_file The path where the compiled report should be saved (e.g., "validation_report.html").
+##' @param template The path to the Quarto template. Defaults to the one provided in the package `inst/reports/Validation_report.qmd`.
+##' 
+##' @author PEcAn Project
+##' @export
+generate_validation_report <- function(benchmark_results, output_file = "Validation_report.html", template = NULL) {
+  PEcAn.logger::logger.info("Generating Validation Benchmark Report...")
+
+  if (is.null(template)) {
+    template <- system.file("reports", "Validation_report.qmd", package = "PEcAn.benchmark")
+    if (template == "") {
+      # Fallback for development mode
+      template <- file.path(getwd(), "inst", "reports", "Validation_report.qmd")
+    }
+  }
+
+  if (!file.exists(template)) {
+    PEcAn.logger::logger.severe("Template file not found:", template)
+    stop("Quarto template not found.")
+  }
+
+  if (!requireNamespace("quarto", quietly = TRUE)) {
+    PEcAn.logger::logger.severe("The 'quarto' package is required to generate the report.")
+    stop("Please install the 'quarto' R package.")
+  }
+
+  # Ensure absolute paths
+  output_file <- normalizePath(output_file, mustWork = FALSE)
+  output_dir <- dirname(output_file)
+
+  if (!dir.exists(output_dir)) {
+    dir.create(output_dir, recursive = TRUE)
+  }
+
+  # Copy template to output directory to avoid permission issues in system folders
+  temp_qmd <- file.path(output_dir, basename(template))
+  file.copy(template, temp_qmd, overwrite = TRUE)
+
+  # Quarto execute_params are converted to YAML. Complex R objects like ggplots
+  # cannot be passed via YAML. We must save them to an RDS and pass the path.
+  results_rds <- file.path(output_dir, "benchmark_results.rds")
+  saveRDS(benchmark_results, results_rds)
+
+  # Render the document
+  tryCatch({
+    quarto::quarto_render(
+      input = temp_qmd,
+      output_file = basename(output_file),
+      execute_params = list(benchmark_results = results_rds)
+    )
+
+    PEcAn.logger::logger.info("Validation report successfully generated at:", output_file)
+  }, error = function(e) {
+    PEcAn.logger::logger.severe("Failed to render validation report:", e$message)
+    stop(e)
+  }, finally = {
+    # Clean up the temporary template file
+    if (file.exists(temp_qmd)) {
+      file.remove(temp_qmd)
+    }
+  })
+
+  return(invisible(output_file))
+}

modules/benchmark/R/metric_residual_plot.R

@@ -10,22 +10,34 @@
 metric_residual_plot <- function(metric_dat, var, filename = NA, draw.plot = is.na(filename)) {
   PEcAn.logger::logger.info("Metric: Residual Plot")
 
-  metric_dat$time <- lubridate::year(as.Date(as.character(metric_dat$time), format = "%Y"))
-  metric_dat$diff <- abs(metric_dat$model - metric_dat$obvs)
-  metric_dat$zeros <- rep(0, length(metric_dat$time))
+  metric_dat <- as.data.frame(metric_dat)
 
-  p <- ggplot2::ggplot(data = metric_dat, ggplot2::aes(x = .data$time)) 
-  p <- p + ggplot2::geom_path(ggplot2::aes(y = .data$zeros), colour = "#666666", size = 2, linetype = 2, lineend = "round") 
-  p <- p + ggplot2::geom_point(ggplot2::aes(y = .data$diff), size = 4, colour = "#619CFF") 
-  p <- p + ggplot2::labs(title = var, x = "years", y = "abs(model - observation)")
+  if (!"time" %in% colnames(metric_dat)) {
+    metric_dat$time <- seq_len(nrow(metric_dat))
+  } else {
+    date.time <- try(as.Date(as.character(metric_dat$time)), silent = TRUE)
+    if (!inherits(date.time, "try-error") && !all(is.na(date.time))) {
+      metric_dat$time <- date.time
+    }
+  }
+
+  # Calculate residuals (Model - Observation)
+  metric_dat$diff <- metric_dat$model - metric_dat$obvs
+
+  p <- ggplot2::ggplot(data = metric_dat, ggplot2::aes(x = .data$time, y = .data$diff)) +
+    ggplot2::geom_hline(yintercept = 0, colour = "#666666", linewidth = 1, linetype = 2) +
+    ggplot2::geom_point(size = 2, alpha = 0.7, colour = "#619CFF") +
+    ggplot2::labs(title = var, x = "Time", y = "Residual (Model - Obs)") +
+    ggplot2::theme_minimal(base_size = 14)
 
   if (!is.na(filename)) {
     grDevices::pdf(filename, width = 10, height = 6)
-    plot(p)
+    print(p)
     grDevices::dev.off()
   }
 
   if (draw.plot) {
     return(p)
   }
+  invisible(p)
 } # metric_residual_plot

modules/benchmark/R/metric_scatter_plot.R

@@ -1,29 +1,32 @@
 ##' Scatter Plot
 ##'
 ##' @param metric_dat dataframe to plot, with at least columns `model` and `obvs`
-##' @param var ignored
+##' @param var title for the plot
 ##' @param filename path to save plot, or NA to not save
 ##' @param draw.plot logical: Return the plot object?
 ##' 
 ##' @author Betsy Cowdery
 ##' @export
-
 metric_scatter_plot <- function(metric_dat, var, filename = NA, draw.plot = is.na(filename)) {
   PEcAn.logger::logger.info("Metric: Scatter Plot")
 
-  p <- ggplot2::ggplot(data = metric_dat) 
-  p <- p + ggplot2::geom_point(ggplot2::aes(x = .data$model, y = .data$obvs), size = 4) 
-  p <- p + ggplot2::geom_abline(slope = 1, intercept = 0, colour = "#666666", 
-                       size = 2, linetype = 2)
+  metric_dat <- as.data.frame(metric_dat)
+
+  p <- ggplot2::ggplot(data = metric_dat, ggplot2::aes(x = .data$model, y = .data$obvs)) +
+    ggplot2::geom_point(size = 2, alpha = 0.7, colour = "#619CFF") +
+    ggplot2::geom_abline(slope = 1, intercept = 0, colour = "#666666", 
+                         linewidth = 1, linetype = 2) +
+    ggplot2::labs(title = var, x = "Modeled", y = "Observed") +
+    ggplot2::theme_minimal(base_size = 14)
 
   if (!is.na(filename)) {
     grDevices::pdf(filename, width = 10, height = 6)
-    plot(p)
+    print(p)
     grDevices::dev.off()
   }
 
   if (draw.plot) {
     return(p)
   }
-
+  invisible(p)
 } # metric_scatter_plot

modules/benchmark/R/metric_timeseries_plot.R

@@ -6,34 +6,61 @@
 ##' @param draw.plot logical: Return the plot object?
 ##'
 ##' @author Betsy Cowdery
-##' @importFrom ggplot2 ggplot labs geom_path geom_point
 ##' @export
-
 metric_timeseries_plot <- function(metric_dat, var, filename = NA, draw.plot = is.na(filename)) {
   PEcAn.logger::logger.info("Metric: Timeseries Plot")
 
-  # Attempt at getting around the fact that time can be annual and thus as.Date won't work
-  date.time <- try(as.Date(metric_dat$time), silent = TRUE)
-  if (inherits(date.time, "try-error")) {
-    PEcAn.logger::logger.warn("Can't coerce time column to Date format, attempting plot anyway")
-  }else{
-    metric_dat$time <- date.time
+  # Ensure metric_dat is a data.frame for ggplot2
+  metric_dat <- as.data.frame(metric_dat)
+
+  if (!"time" %in% colnames(metric_dat)) {
+    PEcAn.logger::logger.warn("Missing 'time' column in metric_dat, using row index instead.")
+    metric_dat$time <- seq_len(nrow(metric_dat))
+  } else {
+    date.time <- try(as.Date(as.character(metric_dat$time)), silent = TRUE)
+    if (!inherits(date.time, "try-error") && !all(is.na(date.time))) {
+      metric_dat$time <- date.time
+    } else {
+      PEcAn.logger::logger.warn("Can't coerce time column to Date format, using original format.")
+    }
+  }
+
+  p <- ggplot2::ggplot(data = metric_dat, ggplot2::aes(x = .data$time))
+
+  # 1. Model Ribbon (if available)
+  if (all(c("model_q05", "model_q95") %in% colnames(metric_dat))) {
+    p <- p + ggplot2::geom_ribbon(ggplot2::aes(ymin = .data$model_q05, ymax = .data$model_q95, fill = "Model 90% CI"), alpha = 0.3)
+  }
+
+  # 2. Model Mean Line
+  p <- p + ggplot2::geom_line(ggplot2::aes(y = .data$model, colour = "Model"), linewidth = 1)
+
+  # 3. Observational Points & Error Bars
+  if ("obvs_sd" %in% colnames(metric_dat)) {
+    p <- p + ggplot2::geom_pointrange(
+      ggplot2::aes(y = .data$obvs, ymin = .data$obvs - .data$obvs_sd, ymax = .data$obvs + .data$obvs_sd, colour = "Observed"), 
+      size = 0.5, alpha = 0.7
+    )
+  } else {
+    p <- p + ggplot2::geom_point(ggplot2::aes(y = .data$obvs, colour = "Observed"), size = 2, alpha = 0.7)
   }
 
-  p <- ggplot(data = metric_dat, ggplot2::aes(x = .data$time)) 
-  p <- p + labs(title = var, y = "") 
-  p <- p + geom_path(ggplot2::aes(y = .data$model, colour = "Model"), size = 2) 
-  p <- p + geom_point(ggplot2::aes(y = .data$model, colour = "Model"), size = 4) 
-  p <- p + geom_path(ggplot2::aes(y = .data$obvs, colour = "Observed"), size = 2) 
-  p <- p + geom_point(ggplot2::aes(y = .data$obvs, colour = "Observed"), size = 4)
+  # Adjust Scales
+  p <- p +
+    ggplot2::scale_colour_manual(values = c("Model" = "#619CFF", "Observed" = "#F8766D")) +
+    ggplot2::scale_fill_manual(values = c("Model 90% CI" = "#619CFF"), name = "Intervals") +
+    ggplot2::labs(title = var, x = "Time", y = "Value", color = "Source") +
+    ggplot2::theme_minimal(base_size = 14) +
+    ggplot2::theme(legend.position = "bottom")
 
   if (!is.na(filename)) {
     grDevices::pdf(filename, width = 10, height = 6)
-    plot(p)
+    print(p)
     grDevices::dev.off()
   }
 
   if (draw.plot) {
     return(p)
   }
+  invisible(p)
 } # metric_timeseries_plot

modules/benchmark/inst/reports/Validation_report.qmd

@@ -0,0 +1,80 @@
+---
+title: "PEcAn Validation Benchmark Report"
+author: "PEcAn Validation Toolkit"
+format:
+  html:
+    theme: cosmo
+    toc: true
+    toc-depth: 3
+    embed-resources: true
+params:
+  benchmark_results: null
+---
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE)
+library(ggplot2)
+```
+
+## Executive Summary
+
+The following scorecard summarizes the validation metrics computed during the benchmark run.
+
+```{r scorecard}
+# Quarto execute_params passes strings. If it's a file path, load the RDS.
+if (is.character(params$benchmark_results) && file.exists(params$benchmark_results)) {
+  results <- readRDS(params$benchmark_results)
+} else {
+  results <- params$benchmark_results
+}
+
+if (!is.null(results) && !is.null(results$metrics)) {
+  # Check if knitr is available to use kable
+  if (requireNamespace("knitr", quietly = TRUE)) {
+    knitr::kable(results$metrics, digits = 3)
+  } else {
+    print(results$metrics)
+  }
+} else {
+  cat("No metrics data provided.\n")
+}
+```
+
+## Visual Diagnostics
+
+```{r plots, results='asis', fig.width=10, fig.height=6}
+if (!is.null(results) && !is.null(results$plots)) {
+  plots <- results$plots
+
+  # if plots is a flat list of ggplot objects
+  if (inherits(plots[[1]], "ggplot")) {
+    for (p_name in names(plots)) {
+      cat(sprintf("\n### %s\n\n", p_name))
+      print(plots[[p_name]])
+      cat("\n\n")
+    }
+  } else {
+    # If plots are nested by variable: plots[[var_name]]$timeseries
+    for (var in names(plots)) {
+      cat(sprintf("\n### Variable: %s\n\n", var))
+
+      var_plots <- plots[[var]]
+
+      if (!is.null(var_plots$timeseries)) {
+        print(var_plots$timeseries)
+        cat("\n\n")
+      }
+      if (!is.null(var_plots$scatter)) {
+        print(var_plots$scatter)
+        cat("\n\n")
+      }
+      if (!is.null(var_plots$residual)) {
+        print(var_plots$residual)
+        cat("\n\n")
+      }
+    }
+  }
+} else {
+  cat("No plots provided.\n")
+}
+```