Downscaling mixed systems

.Rprofile

@@ -1,8 +1,8 @@
 ## --- Custom R repositories (CRAN + R-universe) ------------------------
 options(repos = c(
-    pecanproject = "https://pecanproject.r-universe.dev",
-    ajlyons      = "https://ajlyons.r-universe.dev", # caladaptr
-    CRAN         = "https://cloud.r-project.org"
+  pecanproject = "https://pecanproject.r-universe.dev",
+  ajlyons      = "https://ajlyons.r-universe.dev", # caladaptr
+  CRAN         = "https://cloud.r-project.org"
 ))
 
 # Sys.setenv(R_LIBS_USER = file.path(
@@ -61,4 +61,4 @@ if (!interactive() && requireNamespace("knitr", quietly = TRUE)) {
     dpi = 144
   )
 }
-cat("global Rprofile loaded\n")
+cat("repository Rprofile loaded\n")

.gitignore

@@ -1,33 +1,44 @@
-# R-specific files
+########## Operating system / Editor ##########
+.DS_Store
+Thumbs.db
+.vscode/
+.idea/
+
+########## R / RStudio ##########
 .Rhistory
 .RData
 .Ruserdata
 .Rproj.user/
 *.Rproj
+Rplots.pdf
 
-# Quarto directories
-/.quarto/
+########## Quarto build outputs ##########
 /_site/
 /_book/
-*/docs/
-!/docs/**
 /_freeze/
+/.quarto/
 /_files/
 
+# Keep only top-level docs; ignore nested docs directories elsewhere
+*/docs/
+!/docs/**
 
-# VScode settings
-.vscode/
-
-# Knitr + R Markdown cache
+########## Knitr / R Markdown caches ##########
 *_cache/
 /cache/
 *_files/
-*html
-
-# System files
-.DS_Store
-Thumbs.db
+*.utf8.md
+*.knit.md
+*.html
 
-# Temporary files
-*.tmp
+########## Logs / temp ##########
 *.log
+*.tmp
+nohup.out
+
+########## Deployment / misc ##########
+rsconnect/
+
+########## Generated figures ##########
+# Do not track any rendered images; site publishing copies from _site
+figures/

000-config.R

@@ -1,78 +1,107 @@
 ### Workflow Configuration Settings ###
 
-# Check that we are in the correct working directory
-if(!basename(here::here(getwd())) == 'downscaling') {
-  PEcAn.logger::logger.error("Please run this script from the 'downscaling' directory")
-}
-
-## Global configuration settings for logging
+parser <- argparse::ArgumentParser()
 
-options(
-  # Print all tibble columns
-  tibble.width = Inf,
-  # Suppress readr::read_csv messages
-  readr.show_col_types = FALSE
+## Set development vs production mode ##
+# Dev mode speeds up workflows by subsetting data for testing and debugging
+parser$add_argument("--production",
+  type = "logical", default = TRUE,
+  help = "Set to true for production mode, false for faster development (default: TRUE)"
 )
+args <- parser$parse_args()
+PRODUCTION <- args$production
+
+# Manual override for interactive sessions
+if (rlang::is_interactive()) {
+  PRODUCTION <- TRUE
+}
 
 ## Set parallel processing options
 no_cores <- max(future::availableCores() - 1, 1)
 future::plan(future::multicore, workers = no_cores)
 
-
 # **set ccmmf_dir and pecan_outdir**
 # Define the CCMMF directory from environment variable
 ccmmf_dir <- Sys.getenv("CCMMF_DIR")
 if (ccmmf_dir == "") {
   ccmmf_dir <- "/projectnb2/dietzelab/ccmmf"
+  if (!dir.exists(ccmmf_dir)) {
+    # Fallback to repository root for local development
+    ccmmf_dir <- here::here()
+  }
 }
 pecan_outdir <- file.path(ccmmf_dir, "modelout", "ccmmf_phase_2b_mixed_pfts_20250701")
-# **Is this a test or production run?**
-# Set to FALSE during testing and development
-#
-# Global switch to toggle between fast, small scale runs for development and testing 
-# and full-scale production runs. Works by subsetting various data objects. 
-PRODUCTION <- TRUE
+
+# PEcAn model output to be analyzed
+pecan_archive_tgz <- file.path(ccmmf_dir, "lebauer_agu_2025_20251210.tgz")
 
 # **Variables to extract**
 # see docs/workflow_documentation.qmd for complete list of outputs
-outputs_to_extract <- c(
-    "TotSoilCarb",
-    "AGB"
-)
-
-if(!PRODUCTION) {
+outputs_to_extract <- c("TotSoilCarb", "AGB")
+if (!PRODUCTION) {
   # can subset for testing
   # depending on what part of the workflow you are testing
-  # outputs_to_extract <- outputs_to_extract[1]
+  outputs_to_extract <- outputs_to_extract[1]
 }
 
 ### Configuration Settings that can be set to default ###
 
 # Assume consistent directory structure for other directories
-data_dir     <- file.path(ccmmf_dir, "data")
+data_dir <- file.path(ccmmf_dir, "data")
 raw_data_dir <- file.path(ccmmf_dir, "data_raw")
 cache_dir <- file.path(ccmmf_dir, "cache")
-model_outdir  <- file.path(pecan_outdir, "out")
+model_outdir <- file.path(pecan_outdir, "out")
 
 # Misc
 set.seed(42)
-ca_albers_crs <- 3310
-
-#### Messagees ####
-PEcAn.logger::logger.info("\n\n",
-    "##### SETTINGS SUMMARY #####\n\n",
-    "Running in", ifelse(PRODUCTION, "**production**", "**development**"), "mode\n\n",
-    "### Directory Settings ###\n",
-	"- CCMMF directory:", ccmmf_dir, "\n",
-	"- data_dir       :", data_dir, "\n",
-	"- cache_dir      :", cache_dir, "\n",
-	"- raw_data_dir.  :", raw_data_dir, "\n",
-	"- pecan_outdir.  :", pecan_outdir, "\n",
-	"- model_outdir.  :", model_outdir, "\n\n",
-	"### Other Settings ###\n",
-	"- will extract variables:", paste(outputs_to_extract, collapse = ", "), "\n",
-	"- ca_albers_crs  :", ca_albers_crs, 
-       ifelse(ca_albers_crs == 3310, ", which is NAD83 / California Albers", ""), "\n",
-    wrap = FALSE
-  )
+ca_albers_crs <- "EPSG:3310"
+
+if (terra::is.lonlat(ca_albers_crs)) {
+  PEcAn.logger::logger.severe("`ca_albers_crs` must remain a projected CRS (EPSG:3310).")
+}
+
+ca_albers_info <- terra::crs(ca_albers_crs, describe = TRUE)
+if (is.null(ca_albers_info$code) || ca_albers_info$code != 3310) {
+  PEcAn.logger::logger.severe("`ca_albers_crs` is reserved for California's standard NAD83 / California Albers (EPSG:3310).")
+}
+ca_albers_name <- ca_albers_info[["name"]]
+
+#### Messages ####
+## Ensure pecan_archive_force is defined (allow override via env var)
+if (!exists("pecan_archive_force", inherits = FALSE)) {
+  pecan_archive_force <- isTRUE(as.logical(Sys.getenv("PECAN_ARCHIVE_FORCE", "FALSE")))
+}
+
+msg <- glue::glue(
+  "\n\n",
+  "##### SETTINGS SUMMARY #####\n\n",
+  "Running in {ifelse(PRODUCTION, '**PRODUCTION**', '**DEVELOPMENT**')} mode\n\n",
+  "### Directory Settings ###\n",
+  "- CCMMF directory: {ccmmf_dir}\n",
+  "- data_dir       : {data_dir}\n",
+  "- cache_dir      : {cache_dir}\n",
+  "- raw_data_dir.  : {raw_data_dir}\n",
+  "- pecan_outdir.  : {pecan_outdir}\n",
+  "- model_outdir.  : {model_outdir}\n\n",
+  "- pecan_archive_tgz : {ifelse(is.na(pecan_archive_tgz), '<unset>', pecan_archive_tgz)}\n",
+  "- pecan_archive_force : {pecan_archive_force}\n\n",
+  "### Other Settings ###\n",
+  "- will extract variables: {paste(outputs_to_extract, collapse = ', ')}\n",
+  "- ca_albers_crs : {ca_albers_crs}{if(!is.na(ca_albers_name)) paste0(' (', ca_albers_name, ')') else ''}\n"
+)
+if (!isTRUE(getOption("ccmmf.quiet_banner", FALSE))) {
+  PEcAn.logger::logger.info(msg, wrap = FALSE)
+}
+
+## Global configuration settings for logging
+
+options(
+  # Print all tibble columns
+  tibble.width = Inf,
+  # Suppress readr::read_csv messages
+  readr.show_col_types = FALSE
+)
 
+## Source all R scripts in the R/ directory
+r_scripts <- list.files(file.path(here::here(), "R"), pattern = "\\.R$", full.names = TRUE)
+lapply(r_scripts, source)

R/combine_mixed_crops.R

@@ -0,0 +1,129 @@
+#' Combine two-PFT outputs to represent mixed cropping systems
+#'
+#' Rules for combining woody and annual PFT outputs (stocks or accumulated
+#' flux totals) to represent a mixed cropping system. Supports two methods:
+#'  - "weighted": PFTs partition area (woody_cover + annual_cover = 1) and the
+#'    output is a weighted average: `woody_cover * woody_value + annual_cover * annual_value`.
+#'  - "incremental": preserves the full-area woody baseline (requires `woody_cover == 1`)
+#'    and treats annual as an increment relative to an annual initial baseline:
+#'    `annual_delta = annual_value - annual_init`; `result = woody_value + annual_cover * annual_delta`.
+#'
+#' All inputs must be vectors each of length 1 or a shared common length.
+#'
+#' This function is intended to be applied to pools or integrated fluxes (e.g., annual NPP, annual N2O flux).
+#'
+#' Validation rules (severe errors unless otherwise noted):
+#'  * No input values may be NA (including covers, pool sizes, annual_init if required)
+#'  * Covers must lie within [0,1]
+#'  * Method "incremental": `woody_cover` must be 1 (within 0.1%); if not, severe error
+#'  * Method "incremental": `annual_init` required
+#'  * Method "weighted": rows whose `woody_cover + annual_cover` differ from 1 by more than tolerance
+#'    are set to NA in the result; a single aggregated warning is emitted listing the count
+#'
+#' Inputs may be any quantity expressed per unit area (stocks such as
+#' kg/m2 or fluxes accumulated over a defined time step, e.g., kg/m2 per
+#' hour or year).
+#'
+#' @param woody_value numeric. Pool or accumulated flux for the woody PFT (kg/m2).
+#' @param annual_value numeric. Pool or accumulated flux for the annual PFT (kg/m2).
+#' @param annual_init numeric, required for method = "incremental"; the initial annual pool
+#'        size at t0 (kg C m-2). Used to compute the delta: (annual_value - annual_init).
+#' @param annual_cover numeric. Fractional cover of the annual PFT (0-1).
+#' @param woody_cover numeric. Fractional cover of the woody PFT (0-1). Must be 1 when `method` is "incremental".
+#' @param method character. One of "weighted" or "incremental".
+#'
+#' @return numeric vector of combined values.
+#'
+#' @examples
+#' # Discrete mixing (weights sum to 1)
+#' combine_mixed_crops(
+#'   woody_value = 100, annual_value = 50,
+#'   annual_cover = 0.2, woody_cover = 0.8, method = "weighted"
+#' )
+#'
+#' # Overlap: preserve woody baseline (woody_cover==1), add annual increment scaled by cover
+#' combine_mixed_crops(
+#'   woody_value = 200, annual_value = 220, annual_init = 200,
+#'   annual_cover = 0.3, woody_cover = 1.0, method = "incremental"
+#' )
+combine_mixed_crops <- function(woody_value,
+                                annual_value,
+                                annual_cover,
+                                woody_cover,
+                                annual_init = NULL,
+                                method = c("weighted", "incremental")) {
+  method <- match.arg(method)
+
+  # Internal tolerance for floating point comparisons
+  tol <- 1e-3
+
+  # Accept scalars for cover and vectors for values
+  recycled <- vctrs::vec_recycle_common(
+    woody_value = woody_value,
+    annual_value = annual_value,
+    annual_cover = annual_cover,
+    woody_cover = woody_cover,
+    annual_init = annual_init,
+    .size = vctrs::vec_size_common(
+      woody_value,
+      annual_value,
+      annual_cover,
+      woody_cover,
+      annual_init
+    )
+  )
+
+  woody_value <- recycled$woody_value
+  annual_value <- recycled$annual_value
+  annual_cover <- recycled$annual_cover
+  woody_cover <- recycled$woody_cover
+  annual_init <- recycled$annual_init
+
+  # NA checks (annual_init only if required for incremental)
+  na_checks <- list(
+    woody_value = woody_value,
+    annual_value = annual_value,
+    annual_cover = annual_cover,
+    woody_cover = woody_cover
+  )
+  if (method == "incremental") {
+    if (is.null(annual_init)) {
+      PEcAn.logger::logger.severe("incremental: annual_init is required but missing.")
+    }
+    na_checks$annual_init <- annual_init
+  }
+  na_found <- vapply(names(na_checks), function(nm) anyNA(na_checks[[nm]]), logical(1))
+  if (any(na_found)) {
+    fields <- paste(names(na_found)[na_found], collapse = ", ")
+    PEcAn.logger::logger.severe(paste0("NA values not allowed in inputs: ", fields))
+  }
+
+  # Range checks for covers
+  out_of_range_annual <- (annual_cover < -tol) | (annual_cover > 1 + tol)
+  out_of_range_woody <- (woody_cover < -tol) | (woody_cover > 1 + tol)
+  if (any(out_of_range_annual | out_of_range_woody, na.rm = TRUE)) {
+    n_bad <- sum(out_of_range_annual | out_of_range_woody, na.rm = TRUE)
+    PEcAn.logger::logger.severe("weighted: cover fractions must be in the range [0,1] (+/- tol). ", n_bad, " rows violate this constraint.")
+  }
+
+  if (method == "incremental") {
+    not_one <- abs(woody_cover - 1) > tol
+    if (any(not_one, na.rm = TRUE)) {
+      n_bad <- sum(not_one, na.rm = TRUE)
+      PEcAn.logger::logger.severe("incremental: woody_cover must be 1 (+/- ", tol, "); ", n_bad, " rows violate.")
+    }
+    res <- woody_value + annual_cover * (annual_value - annual_init)
+    return(res)
+  }
+
+  # weighted method
+  sum_cov <- woody_cover + annual_cover
+  bad_sum <- abs(sum_cov - 1) > tol | is.na(sum_cov)
+  res <- woody_cover * woody_value + annual_cover * annual_value
+  if (any(bad_sum, na.rm = TRUE)) {
+    n_bad <- sum(bad_sum, na.rm = TRUE)
+    PEcAn.logger::logger.warn(paste0("weighted: ", n_bad, " rows with cover fractions not summing to 1 (<U+00B1>tol); results set to NA for those rows."))
+    res[bad_sum] <- NA_real_
+  }
+  return(res)
+}