Merge branch 'develop' into automated-reporting-and-visualization

Author: ayushman1210

CHANGELOG.md

@@ -31,6 +31,7 @@ For more information about this file see also [Keep a Changelog](http://keepacha
 - PEcAn.SIPNET gains support for SIPNET v2, whose features includes management events, nitrogen cycle tracking, explicit N2O and methane fluxes, runtime setting of feature flags, and changes to the parameter set (now 73 parameters). SIPNET v1 is still fully supported, but workarounds for bugs in the legacy `sipnet.unk` version have been removed.
 - Added `PEcAn.data.land::to_co2e()` for converting SOC change, CH4, and N2O to CO2-equivalent emissions using IPCC Global Warming Potential values.
 - Added `PEcAn.data.land::event_parquet_to_json` for generating PEcAn `event.json` files from well-formatted event parquet files, with support for ensembles of events.
+- Added statewide synthetic fertilization and compost amendment event workflows for CA ag parcels. Outputs share an ensemble naming so a downstream cleaner unions them into one fertilization event type for SIPNET.
 
 ### Fixed
 - Docker GHA workflow no longer fails on pull requests opened from forks (#3618).
@@ -42,6 +43,7 @@ For more information about this file see also [Keep a Changelog](http://keepacha
 - `segment_dataframe()` now returns an empty dataframe when date filtering removes all crop-cycle segments, instead of a single row with NA columns that caused downstream segment config errors (#4007).
 
 ### Changed
+- Added `ensemble_downscale()`, a refactored version of `SDA_downscale()`.
 - `PEcAn.uncertainty::get.parameter.samples()`: replaced the `save_to_disk` flag (from #3860) with an `outdir` argument (default `settings$outdir`) controlling whether `samples.Rdata` is written; `outdir = NULL` skips the save. Existing callers are unaffected (@omkarrr2533, #4016)
 - Updated Docker architecture documentation to match current docker-compose.yml: removed portainer/minio/thredds, added rstudio/api sections, updated service lists and volumes (#3268).
 - Improved PEcAn.SIPNET documentation including README, model description, and current installation instructions (@Eshaan-byte; #3703, #3705).

base/workflow/R/runModule.run.write.configs.R

@@ -262,16 +262,15 @@ runModule.run.write.configs <- function(settings,
   }
 
   # Deprecation: internal design generation is going away. Passing input_design
-  # explicitly (the generate_joint_ensemble_design() result) will become the
-  # required path. Warn only when we are actually about to auto-generate.
-  auto_generating <-
-    (is.null(designs$ensemble) && need_ensemble) ||
-    (is.null(designs$sensitivity) && need_sa)
-  if (auto_generating) {
+  # explicitly will become the required path. Warn only about the design we are
+  # actually about to generate, so a caller who supplied one is not told to
+  # supply it again.
+  if (is.null(designs$ensemble) && need_ensemble) {
     PEcAn.logger::logger.warn(
-      "Internal input design generation is deprecated and will be removed.",
-      "Pass input_design explicitly as the list(design_matrix, samples) returned",
-      "by generate_joint_ensemble_design(); this will become required."
+      "Generating the ensemble design internally is deprecated and will be",
+      "removed. Pass input_design explicitly as the list(design_matrix,",
+      "samples) returned by generate_joint_ensemble_design(); this will",
+      "become required."
     )
   }
 

base/workflow/tests/testthat/test-runModule.run.write.configs.R

@@ -286,4 +286,24 @@ test_that(".prepare_input_designs still accepts a design supplied as X", {
 
   mockery::expect_called(gen, 0)
   expect_identical(designs$ensemble, supplied$X)
+})
+
+test_that(".prepare_input_designs does not warn about a design the caller supplied", {
+  tmp <- withr::local_tempdir()
+  settings <- make_prep_settings(tmp)
+  settings$sensitivity.analysis <- list(quantiles = c(0.025, 0.5, 0.975))
+
+  supplied <- list(design_matrix = data.frame(param = 1:3), samples = fake_bundle())
+
+  mockery::stub(.prepare_input_designs,
+                "PEcAn.uncertainty::generate_OAT_SA_design",
+                function(...) list(design_matrix = data.frame(param = 1:4)))
+
+  # the SA design is still generated internally, but there is no way to supply
+  # one, so a caller who did pass a design should not be told to pass one
+  msgs <- capture.output(
+    invisible(.prepare_input_designs(settings, input_design = supplied)),
+    type = "message"
+  )
+  expect_false(any(grepl("deprecated", msgs)))
 })

documentation/tutorials/Demo_02_Uncertainty_Analysis/uncertainty.qmd

@@ -121,13 +121,24 @@ See Demo 1 Section 6 for details on what these functions do. Briefly, they read
 
 # Write Model Configuration Files {#sec-write-configs}
 
-This step generates the model-specific configuration files that will be used to run the ecosystem model. The process involves:
+This step decides which model runs to do, then writes the configuration files for them.
+
+`generate_joint_ensemble_design()` builds the ensemble design: one row per run, holding indices that say which parameter draw each run uses. It returns those samples alongside the design, since the indices only mean anything together with the samples they point into. Because this demo also runs a sensitivity analysis, the SA design is generated internally as well; passing one in yourself isn't supported yet.
 
-2. Generating SIPNET configuration files using the `runModule.run.write.configs()` function.
 ```{r write-configs}
-settings <- PEcAn.workflow::runModule.run.write.configs(settings)
+design <- PEcAn.uncertainty::generate_joint_ensemble_design(
+  settings,
+  ensemble_size = settings$ensemble$size
+)
+
+settings <- PEcAn.workflow::runModule.run.write.configs(
+  settings,
+  input_design = design
+)
 ```
 
+Generating the design yourself lets you inspect it before any runs happen: `design$design_matrix` is the table of runs, and `design$samples` holds the parameter draws behind it, which is what the sensitivity and variance decomposition steps later report on.
+
 # Run Model Simulations
 
 This section executes the SIPNET simulations and retrieves the results.

documentation/tutorials/Demo_1_Basic_Run/pecan.xml

@@ -15,6 +15,7 @@
   </pft>
  </pfts> 
  <ensemble>
+   <size>10</size>
    <variable>NPP</variable>
   </ensemble>
  <model>

documentation/tutorials/Demo_1_Basic_Run/run_pecan.qmd

@@ -39,7 +39,7 @@ This notebook demonstrates how to:
 
 We are modeling carbon and productivity dynamics at the Niwot Ridge Forest AmeriFlux site ([US-NR1](https://ameriflux.lbl.gov/sites/siteinfo/US-NR1), a high-elevation temperate coniferous forest in Colorado. The model configuration uses the SIPNET process-based ecosystem model, parameterized with a temperate coniferous plant functional type (PFT).
 
-The simulation is run for the full year 2004 (January 1 – December 31) using AmeriFlux LBL meteorological drivers from the Niwot Ridge site. The ensemble setup specifies one model run focusing on net primary productivity (NPP) as the target output variable.
+The simulation is run for the full year 2004 (January 1 – December 31) using AmeriFlux LBL meteorological drivers from the Niwot Ridge site. The ensemble setup specifies a small ten member ensemble focusing on net primary productivity (NPP) as the target output variable, so the plots later on show the spread across members rather than a single line.
 
 This scenario is designed to be a minimal, reproducible example to demonstrate how to run SIPNET within the PEcAn workflow. In later steps, this same framework can be extended to include more ensemble members, additional PFTs, longer time periods, or alternative meteorological inputs.
 
@@ -192,12 +192,26 @@ Additional outputs include logs, a `STATUS` file that records the steps of the w
 
 # Write Model Configuration Files
 
-This step generates the model-specific configuration files and scripts that will be used to run the ecosystem model. The process involves generating SIPNET configuration files using the `runModule.run.write.configs()` function.
+This step decides which model runs to do, then writes the configuration files for them.
+
+`generate_joint_ensemble_design()` builds the design: one row per run, holding indices that say which parameter draw and which input file each run uses. It returns those samples alongside the design, since the indices only mean anything together with the samples they point into.
+
+`runModule.run.write.configs()` then takes that whole object and writes the SIPNET configuration files.
 
 ```{r write-configs}
-settings <- PEcAn.workflow::runModule.run.write.configs(settings)
+design <- PEcAn.uncertainty::generate_joint_ensemble_design(
+  settings,
+  ensemble_size = settings$ensemble$size
+)
+
+settings <- PEcAn.workflow::runModule.run.write.configs(
+  settings,
+  input_design = design
+)
 ```
 
+You can call `runModule.run.write.configs()` without a design and it will build one for itself, but that path is deprecated and will warn. Generating the design yourself also lets you inspect it before any runs happen: `design$design_matrix` is the table of runs, and `design$samples` holds the parameter draws behind it.
+
 # Run Model Simulations and Fetch Results
 
 This section executes the actual model simulations and retrieves the results. The process is managed by PEcAn's workflow system, which handles the execution of your chosen ecosystem model.
@@ -215,27 +229,68 @@ This step generates raw model outputs in model-specific format (in this case, `s
 
 After the model simulation completes, we need to extract the results and prepare them for analysis. This involves:
 
-1. Reading the run ID
+1. Reading the run IDs
 2. Setting up output paths
 3. Defining time period
 4. Loading model output
 5. Convert to a standard format
 
 ```{r get-plot-vars}
-runid <- as.character(read.table(paste(settings$outdir, "/run/", "runs.txt", sep = ""))[1, 1]) # Note: if you are using an xml from a run with multiple ensembles this line will provide only the first run id
-outdir <- file.path(settings$outdir, "/out/", runid)
+runids <- as.character(read.table(paste(settings$outdir, "/run/", "runs.txt", sep = ""))[, 1])
 start.year <- lubridate::year(settings$run$start.date)
 end.year <- lubridate::year(settings$run$end.date)
-model_output <- PEcAn.utils::read.output(
-  runid,
-  outdir,
-  start.year,
-  end.year,
-  variables = NULL,
-  dataframe = TRUE,
-  verbose = FALSE
-)
+
+# one entry per ensemble member, keyed by run id
+ensemble_output <- lapply(runids, function(runid) {
+  PEcAn.utils::read.output(
+    runid,
+    file.path(settings$outdir, "/out/", runid),
+    start.year,
+    end.year,
+    variables = NULL,
+    dataframe = TRUE,
+    verbose = FALSE
+  )
+})
+names(ensemble_output) <- runids
+
+# the members share a variable set, so the first one describes them all
+model_output <- ensemble_output[[1]]
 available_vars <- names(model_output)[!names(model_output) %in% c("posix", "time_bounds")]
+
+# sub-daily fluxes are too spiky to read as ten overlaid lines, so keep a
+# daily-averaged copy of each member for those plots
+daily_output <- lapply(ensemble_output, function(member) {
+  day <- as.Date(member$posix)
+  vars <- setdiff(names(member), c("posix", "time_bounds"))
+  averaged <- lapply(member[vars], function(x) as.numeric(tapply(x, day, mean, na.rm = TRUE)))
+  names(averaged) <- vars
+  c(list(posix = as.POSIXct(sort(unique(day)), tz = "UTC")), averaged)
+})
+
+# draw one line per ensemble member, so the spread across the ensemble is visible
+plot_ensemble <- function(vars, cols, ylab, main, labels = vars, data = ensemble_output) {
+  yrange <- range(
+    unlist(lapply(data, function(member) unlist(member[vars]))),
+    na.rm = TRUE
+  )
+  # empty plot drawn off the real time axis, so the x axis stays dates
+  plot(
+    data[[1]]$posix,
+    data[[1]][[vars[1]]],
+    type = "n",
+    ylim = yrange,
+    xlab = "Date",
+    ylab = ylab,
+    main = main
+  )
+  for (i in seq_along(vars)) {
+    for (member in data) {
+      lines(member$posix, member[[vars[i]]], col = adjustcolor(cols[i], alpha.f = 0.4))
+    }
+  }
+  legend("topright", legend = labels, col = cols, lty = 1)
+}
 ```
 
 Running this code will convert model specific output files into a standardized netCDF ([year].nc) that can be downloaded for visualization and analysis (R, Matlab, ncview, panoply, etc). This is a key step, because this standardization enables PEcAn to apply downstream analyses to outputs from different ecosystem models.
@@ -267,64 +322,67 @@ knitr::kable(vars_df, caption = "Model Output Variables and Descriptions")
 
 This section provides examples of how to create time series plots of different model variables. The examples cover various ecosystem processes including carbon fluxes, carbon pools, water variables, and structural variables like Leaf Area Index (LAI).
 
+Each plot draws one line per ensemble member. The members differ only in their parameter draws, so the spread between the lines is the parameter uncertainty carried through the model.
+
 ## Plot Carbon Fluxes
 
+Some members run negative NPP even in the growing season. That happens when the drawn respiration parameters outweigh photosynthesis, and it shows up here because this demo samples from priors rather than from a meta-analysis, so the spread is wide. Demo 3 shows what constraining these parameters with data does to that spread.
+
+GPP and NPP are sub-daily fluxes, so ten overlaid members would just fill the panel. These are averaged to daily values first, which is enough to see the members apart.
+
 ```{r plot-carbon-fluxes}
-# Plot Gross Primary Productivity (GPP) and Net Primary Productivity (NPP)
-plot(model_output$posix, model_output$GPP,
-  type = "l",
-  col = "green",
-  xlab = "Date",
+plot_ensemble(
+  vars = c("GPP", "NPP"),
+  cols = c("green", "blue"),
   ylab = "Carbon Flux (kg C m-2 s-1)",
-  main = "Carbon Fluxes Over Time"
+  main = "Carbon Fluxes Over Time (daily average)",
+  data = daily_output
 )
-lines(model_output$posix, model_output$NPP, col = "blue")
-legend("topright", legend = c("GPP", "NPP"), col = c("green", "blue"), lty = 1)
 ```
 
 ## Plot Carbon Pools
 
 ```{r plot-carbon-pools}
-# Plot Total Live Biomass and Total Soil Carbon
-plot(model_output$posix, model_output$TotLivBiom,
-  type = "l",
-  col = "darkgreen",
-  xlab = "Date",
+plot_ensemble(
+  vars = c("TotLivBiom", "TotSoilCarb"),
+  cols = c("darkgreen", "brown"),
+  labels = c("Total Live Biomass", "Total Soil Carbon"),
   ylab = "Carbon Pool (kg C m-2)",
   main = "Carbon Pools Over Time"
 )
-lines(model_output$posix, model_output$TotSoilCarb, col = "brown")
-legend("topright", legend = c("Total Live Biomass", "Total Soil Carbon"), col = c("darkgreen", "brown"), lty = 1)
 ```
 
 ## Plot Water Variables
 
+Soil moisture separates into ten lines once the growing season starts, but snow water equivalent stays a single line. That is expected: snow is driven by the met input, which is the same file for every member, while soil moisture responds to the sampled parameters.
+
 ```{r plot-water-variables}
-# Plot Soil Moisture and Snow Water Equivalent
-plot(model_output$posix, model_output$SoilMoist,
-  type = "l",
-  col = "blue",
-  xlab = "Date",
-  ylab = "Soil Moisture (kg m-2)",
-  main = "Soil Moisture Over Time"
+plot_ensemble(
+  vars = c("SoilMoist", "SWE"),
+  cols = c("blue", "lightblue"),
+  labels = c("Soil Moisture", "Snow Water Equivalent"),
+  ylab = "Water (kg m-2)",
+  main = "Soil Moisture and Snow Water Equivalent Over Time"
 )
-lines(model_output$posix, model_output$SWE, col = "lightblue")
-legend("topright", legend = c("Soil Moisture", "Snow Water Equivalent"), col = c("blue", "lightblue"), lty = 1)
 ```
 
-## Plot LAI and Biomass
+## Plot LAI and Above Ground Wood
 
 ```{r plot-lai-biomass}
-# Plot Leaf Area Index (LAI) and Above Ground Wood
-plot(model_output$posix, model_output$LAI,
-  type = "l",
-  col = "darkgreen",
-  xlab = "Date",
+plot_ensemble(
+  vars = "LAI",
+  cols = "darkgreen",
   ylab = "LAI (m2 m-2)",
   main = "Leaf Area Index Over Time"
 )
-lines(model_output$posix, model_output$AbvGrndWood, col = "brown")
-legend("topright", legend = c("LAI", "Above Ground Wood"), col = c("darkgreen", "brown"), lty = 1)
+
+plot_ensemble(
+  vars = "AbvGrndWood",
+  cols = "brown",
+  labels = "Above Ground Wood",
+  ylab = "Wood (kg C m-2)",
+  main = "Above Ground Wood Over Time"
+)
 ```
 
 # Conclusion

modules/assim.sequential/DESCRIPTION

@@ -43,6 +43,7 @@ Suggests:
     glue,
     ggpubr,
     gridExtra,
+    keras3 (>= 1.0.0),
     itertools,
     magic (>= 1.5.0),
     methods,
@@ -55,18 +56,17 @@ Suggests:
     plotrix,
     plyr (>= 1.8.4),
     randomForest,
-    keras3 (>= 1.0.0),
     raster,
     readr,
     reshape2 (>= 1.4.2),
     rlist,
     sf,
+    sp,
     stats,
     terra,
     testthat,
     tictoc,
     tidyr,
-    sp,
     utils,
     xgboost,
     XML

modules/assim.sequential/NAMESPACE

@@ -29,8 +29,10 @@ export(assessParams)
 export(block_matrix)
 export(conj_wt_wishart_sampler)
 export(construct_nimble_H)
+export(downscale_metrics)
 export(downscale_qsub_main)
 export(dwtmnorm)
+export(ensemble_downscale)
 export(get_ensemble_weights)
 export(hop_test)
 export(interactive.plotting.sda)
@@ -39,7 +41,7 @@ export(load_data_paleon_sda)
 export(matrix_network)
 export(metSplit)
 export(obs_timestep2timepoint)
-export(outlier.detector.boxplot)
+export(outlier_detector_boxplot)
 export(piecew.poly.local)
 export(post.analysis.ggplot)
 export(post.analysis.ggplot.violin)
@@ -62,6 +64,7 @@ export(sda_assemble)
 export(sda_weights_site)
 export(simple.local)
 export(stack_covariates_2_geotiff)
+export(subset_ensemble)
 export(tobit.model)
 export(tobit2space.model)
 export(tobit_model_censored)