Merge pull request #4040 from ayushman1210/automated-reporting-and-visualization

feat: implement automated reporting and visualization layer

Author: dlebauer

docker/depends/pecan_package_dependencies.csv

@@ -481,6 +481,7 @@
 "purrr",">= 0.2.3","base/workflow","Imports",FALSE
 "purrr",">= 0.2.3","modules/data.atmosphere","Imports",FALSE
 "pwr","*","modules/rtm","Suggests",FALSE
+"quarto","*","modules/benchmark","Suggests",FALSE
 "R.utils","*","base/db","Imports",FALSE
 "randomForest","*","modules/assim.sequential","Suggests",FALSE
 "randtoolbox","*","base/utils","Suggests",FALSE

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/align_data.R

@@ -117,4 +117,4 @@ align_data <- function(model.calc, obvs.calc, var, align_method = "match_timeste
   }
 
   return(dat)
-} # align_data
+} # align_data

modules/benchmark/R/generate_validation_report.R

@@ -0,0 +1,64 @@
+##' 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)
+  }
+  
+  if (!requireNamespace("quarto", quietly = TRUE)) {
+    PEcAn.logger::logger.severe("The 'quarto' package is required to generate the report.")
+  }
+  
+  # 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)
+  }, finally = {
+    # Clean up the temporary template file
+    if (file.exists(temp_qmd)) {
+      file.remove(temp_qmd)
+    }
+  })
+  
+  return(invisible(output_file))
+}

modules/benchmark/R/metric_Coverage.R

@@ -1,22 +1,34 @@
 #' @name metric_Coverage
 #' @title Prediction Interval Coverage
 #' @export
-#' @param dat dataframe with columns `model_q025` and `model_q975`
+#' @param dat dataframe with columns `model_q05` and `model_q95` (or `model_q025` and `model_q975`)
 #' @param ... ignored
-#' @return A numeric value representing the fraction of observations that fall within the 95% prediction interval.
+#' @return A numeric value representing the fraction of observations that fall within the prediction interval.
 #' @details
 #' Measures the fraction of observations that fall within the model's
-#' stated 95% prediction interval.
+#' stated prediction interval (defaults to 90% interval via `model_q05`/`model_q95`, or 95% via `model_q025`/`model_q975`).
 
 metric_Coverage <- function(dat, ...) {
-  if (!"model_q025" %in% names(dat) || !"model_q975" %in% names(dat)) {
-    PEcAn.logger::logger.severe("Metric Coverage requires 'model_q025' and 'model_q975' columns in the dataset.")
+  q_low <- NULL
+  q_high <- NULL
+  
+  if (all(c("model_q05", "model_q95") %in% names(dat))) {
+    q_low <- dat$model_q05
+    q_high <- dat$model_q95
+  } else if (all(c("model_q025", "model_q975") %in% names(dat))) {
+    q_low <- dat$model_q025
+    q_high <- dat$model_q975
+  } else {
+    PEcAn.logger::logger.severe("Metric Coverage requires quantile columns ('model_q05'/'model_q95' or 'model_q025'/'model_q975') in the dataset.")
   }
 
   PEcAn.logger::logger.info("Metric: Prediction Interval Coverage")
 
-  valid <- !is.na(dat$obvs) & !is.na(dat$model_q025) & !is.na(dat$model_q975)
-  covered <- dat$obvs[valid] >= dat$model_q025[valid] & dat$obvs[valid] <= dat$model_q975[valid]
-  
+  valid <- !is.na(dat$obvs) & !is.na(q_low) & !is.na(q_high)
+  if (!any(valid)) {
+    return(NA_real_)
+  }
+
+  covered <- dat$obvs[valid] >= q_low[valid] & dat$obvs[valid] <= q_high[valid]
   return(mean(covered))
 }

modules/benchmark/R/metric_residual_plot.R

@@ -2,30 +2,105 @@
 ##'
 ##' @param metric_dat dataframe to plot, with at least columns `time`, `model`, `obvs`
 ##' @param var variable name, used as plot title
+##' @param unit measurement unit for the variable, added to the y-axis label; NULL to omit
 ##' @param filename path to save plot, or NA to not save
 ##' @param draw.plot logical: Return the plot object?
 ##' 
 ##' @author Betsy Cowdery
 ##' @export
-metric_residual_plot <- function(metric_dat, var, filename = NA, draw.plot = is.na(filename)) {
+metric_residual_plot <- function(metric_dat, var, unit = NULL, 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
+  
+  is_multi_site <- "site" %in% colnames(metric_dat) && length(unique(metric_dat$site)) > 1
+  is_multi_var <- "variable" %in% colnames(metric_dat) && length(unique(metric_dat$variable)) > 1
+  
+  if (!"site" %in% colnames(metric_dat)) metric_dat$site <- "All"
+  if (!"variable" %in% colnames(metric_dat)) metric_dat$variable <- var
+  
+  facet_groups <- split(metric_dat, list(metric_dat$site, metric_dat$variable), drop = TRUE)
+  
+  fit_results <- lapply(names(facet_groups), function(g) {
+    sub_dat <- facet_groups[[g]]
+    
+    time_num <- as.numeric(sub_dat$time)
+    if (inherits(sub_dat$time, "POSIXt")) {
+      time_num <- time_num / 86400
+    }
+    time_num <- time_num - min(time_num, na.rm = TRUE)
+    sub_dat$time_num <- time_num
+    
+    fit <- try(stats::lm(diff ~ time_num, data = sub_dat), silent = TRUE)
+    if (!inherits(fit, "try-error") && length(stats::coef(fit)) == 2 && !is.na(stats::coef(fit)[2])) {
+      intercept <- stats::coef(fit)[1]
+      slope <- stats::coef(fit)[2]
+      label_str <- sprintf("Intercept: %.3f\nSlope: %.4f / day", intercept, slope)
+      sub_dat$trend <- stats::fitted(fit)
+    } else {
+      label_str <- "Trend: N/A"
+      sub_dat$trend <- NA_real_
+    }
+    
+    annot <- data.frame(
+      site = sub_dat$site[1],
+      variable = sub_dat$variable[1],
+      label = label_str
+    )
+    
+    list(sub_dat = sub_dat, annot = annot)
+  })
+  
+  metric_dat <- do.call(rbind, lapply(fit_results, `[[`, "sub_dat"))
+  annotations <- do.call(rbind, lapply(fit_results, `[[`, "annot"))
+  
+  ylab <- if (is.null(unit)) "residual (model - obs)" else sprintf("residual (%s)", unit)
+  
+  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::geom_line(ggplot2::aes(y = .data$trend), colour = "#FF3333", linetype = "dashed", na.rm = TRUE) +
+    ggplot2::labs(title = var, x = "time", y = ylab, colour = NULL, fill = NULL) +
+    ggplot2::theme_minimal(base_size = 12)
+    
+  # Add per-panel annotations
+  p <- p + ggplot2::geom_label(
+    data = annotations,
+    ggplot2::aes(x = -Inf, y = Inf, label = .data$label),
+    hjust = -0.05, vjust = 1.1,
+    inherit.aes = FALSE,
+    alpha = 0.8
+  )
+  
+  # Add facets if applicable
+  if (is_multi_site && is_multi_var) {
+    p <- p + ggplot2::facet_wrap(~ variable + site, scales = "free_y")
+  } else if (is_multi_site) {
+    p <- p + ggplot2::facet_wrap(~ site, scales = "free_y")
+  } else if (is_multi_var) {
+    p <- p + ggplot2::facet_wrap(~ variable, scales = "free_y")
+  }
 
   if (!is.na(filename)) {
     grDevices::pdf(filename, width = 10, height = 6)
-    plot(p)
+    print(p)
     grDevices::dev.off()
   }
 
   if (draw.plot) {
     return(p)
   }
-} # metric_residual_plot
+  invisible(p)
+} # metric_residual_plot

modules/benchmark/R/metric_scatter_plot.R

@@ -1,29 +1,36 @@
 ##' 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 unit measurement unit for the variable, added to the axis labels; NULL to omit
 ##' @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)) {
+metric_scatter_plot <- function(metric_dat, var, unit = NULL, 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)
+  
+  xlab <- if (is.null(unit)) "modeled" else sprintf("modeled (%s)", unit)
+  ylab <- if (is.null(unit)) "observed" else sprintf("observed (%s)", unit)
+  
+  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 = xlab, y = ylab, colour = NULL, fill = NULL) +
+    ggplot2::theme_minimal(base_size = 12)
 
   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