How to handle and organize multi-site ensemble runs?

[dlebauer]

Overview

@infotroph and I were discussing how to put together multi-site ensemble forecast outputs. Current behavior is that

• PEcAn standard netcdf files are organized thus

|- out/
  - ENS-<ens number>-<site id>/
    - <YYYY>.nc

And PEcAn generates summary Rdata files including

- ensemble.<ensemble_id>.<varname>.<startYYYY-endYYYY>.Rdata
# note that ensemble_id != ens number

But for these Rdata files (and associated PDF summary plots), each site overwrites the previous site's output.

So we were looking for a way to organize these outputs by site.

CMS project approach

For the CMS project, outputs are organized thus¹:

• sda.outputN.Rdata: this series of files contains the general SDA outputs for the Nth time step. Within each file, the obs.mean contains the mean observations for each state variable of each site, and the obs.cov contains the observation covariance matrix for each state variable of each site. The forecast includes a matrix of forecast results with the dimensions of the ensemble number (rows) and the state variable number times site number (columns) (e.g., if we have 25 ensembles for each site, four state variables, and 6400 sites, then the matrix dimensions will be 25*(46400) = 2525600). Every four columns correspond to the ensemble forecasts of every state variable for each site. ... The [[6]] object contains the relative ensemble weights calculated by comparing analysis and forecast results. The param.list contains the parameters used for each ensemble member of each site.
• sda.all.forecast.analysis.Rdata: This file summarizes the forecast and analysis matrix across all runs.

sda.all.forecast.analysis.Rdata contains analysis.all and forecast.all, which are lists of matrices (data.frames, one per time step) that can be used for downscaling.

These data.frames look like:

> head(forecast.all[[1]][,1:12], n = 2)
  AbvGrndWood LAI SoilMoistFrac TotSoilCarb AbvGrndWood.1 LAI.1 SoilMoistFrac.1
1      5.0004   0          50.6    64.18524        1.6277     0            52.4
2      5.4908   0          48.4     0.98392        2.0679     0            52.4
  TotSoilCarb.1 AbvGrndWood.2 LAI.2 SoilMoistFrac.2 TotSoilCarb.2
1      30.81041        1.4220     0            52.0      11.70966
2      25.42046        1.7641     0            48.4      60.87190

Questions:

@dongchen Zhang

• What scripts creates these files?
• How do you align sites with the columns in the matrix?

@PecanProject (community):

• Generally, what would be a good way (or are there current implementations) to handle ensemble.*.Rdata and associated files from individual sites that get clobbered during multisite runs?
• Any suggested refinements to the approach above?

Footnotes

1. from geo.bu.edu:/projectnb/dietzelab/dongchen/anchorSites/NA_runs/SDA_25ens_2024_11_25/ReadMe.txt ↩

[mdietze]

• The SDA workflow generates the analysis and forecast files, they aren't generated by separate scripts that post-process the raw output
• For the alignment of sites with columns, that info is in a separate file. I've asked @DongchenZ many times to embed this info in the file itself to reduce errors/confusion. Also, that output would really be better stored as a 3-D array (site, variable, ensemble member) but SDA itself requires unwrapping a site x variable matrix into a vector, and no one ever took the extra step of reorganizing that confusing format into an array as part of the SDA output, which admittedly would take extra processing time that is arguably better spent during a postprocessing phase
• The ensemble.*.Rdata file format was never updated for multisite because no one ever needed it and no one was running the workflow step that generates it as part of their workflow. Before deciding HOW to modify that step, I think the top priority would be to discuss IF we even need that step in the first place since in the 9 years that MultiSettings has existed no one has ever needed this functionality.
• I think the original ensemble post processing functions never really served their intended purpose because they would only process a pre-specified variable as part of a monolithic workflow. I think functions that read PEcAn outputs on the fly, such as read.outputs, are much more useful.
• My strong prediction is that if you try to merge all the outputs into one file that file will be too big to work with effectively. You'll then shift to cloud-native file formats that tile the date, which is basically where we are already but with a few bells and whistles that help make the looping over files invisible. One alternative would be to tweak all of pecans outputs to be cloud native to begin with, which is a great/noble end goal but a lot of work in the short term. The other would be to have some helper functions that would allow you to more easily read ensemble outputs on the fly
• If you think it's worth pre-processsing ensemble outputs to a lower dimensional analysis-ready file (which is what Dongchen inadvertently ended up with on the SDA outputs) then (1) I'd recommend using the EFI standard and (2) you'd have to decide ahead of time what variables and what timestep you want in that summary file.

[dlebauer]

Thanks - here are two formats I've come up with. Let me know if you have any suggestions. (code is in 02_extract_sipnet_output.R; part of ccmmf/workflows#2

Note, for now I'm just taking the last time point. Once (if) we add a time dimension, the arrays will get more complex!

1. Analysis ready

> dim(ens_array)
[1] 5 2 5
> dimnames(ens_array)
$site
[1] "84389254458b3a58" "9e3b1749179a4a68" "729a6ef989c5a277" "8c23aa17d5b211ed"
[5] "57ea98e608333b76"

$variable
[1] "AGB"         "TotSoilCarb"

$ensemble
[1] "00001" "00002" "00003" "00004" "00005"

> 
> ens_array
, , ensemble = 00001

                  variable
site                    AGB TotSoilCarb
  84389254458b3a58 10.78223    10.20318
  9e3b1749179a4a68  0.81106     4.50113
  729a6ef989c5a277  0.04493     7.05939
  8c23aa17d5b211ed  0.93889     5.41103
  57ea98e608333b76  0.12848     6.57558

, , ensemble = 00002

                  variable
site                    AGB TotSoilCarb
  84389254458b3a58 12.57011    11.05740
  9e3b1749179a4a68  2.34102     4.61618
  729a6ef989c5a277  0.44770     4.18219
  8c23aa17d5b211ed  1.05084     7.69727
  57ea98e608333b76  0.52968     4.45944

, , ensemble = 00003

                  variable
site                    AGB TotSoilCarb
  84389254458b3a58 14.88861    10.36968
  9e3b1749179a4a68  1.87063     4.17231
  729a6ef989c5a277  0.64939     5.95299
  8c23aa17d5b211ed  1.62388     4.33456
  57ea98e608333b76  0.56853     4.22801

, , ensemble = 00004

                  variable
site                    AGB TotSoilCarb
  84389254458b3a58 15.38742    10.46131
  9e3b1749179a4a68  0.00540     3.25759
  729a6ef989c5a277  0.71515     3.34682
  8c23aa17d5b211ed  0.09193     8.56669
  57ea98e608333b76  0.09600     4.78530

, , ensemble = 00005

                  variable
site                    AGB TotSoilCarb
  84389254458b3a58 13.03796     8.13995
  9e3b1749179a4a68  0.42291     5.45559
  729a6ef989c5a277  0.11499     3.96766
  8c23aa17d5b211ed  0.00019     5.79587
  57ea98e608333b76  0.44818     4.63800

> 

2. EFI Std

> efi_std
# A tibble: 1,960 × 7
   time                ensemble     X     Y TotSoilCarb     AGB obs_flag
   <dttm>              <chr>    <dbl> <dbl>       <dbl>   <dbl>    <dbl>
 1 2024-01-01 23:00:00 00001    -121.  37.3        5.96 0.396          0
 2 2024-01-01 23:00:00 00001    -122.  38.2        5.96 4.58           0
 3 2024-01-01 23:00:00 00001    -122.  37.7        5.71 0.163          0

[mdietze]

So both options are close to the EFI standard (which allows both long and array options) but neither is exactly compliant. First, I'll note that I was confused that the first example didn't have time but the second did -- I'm going to assume there IS a time dimension, but if not just drop from what I'm recommending. If you go with the array option my recommendation would be [time, site, ensemble, variable] or separate [time, site, ensemble] arrays for each variable (the latter being more normal for netCDF, and probably more efficient on memory as you need maps and timeseries of a single variables much more ofter than you need to calculate across-variable covariances). If you go with the long option, my recommendation would be [time, site, ensemble, variable, prediction]. In both cases I think site id is a better choice than lon/lat since it eliminates a potentially error prone and computationally expensive distance-based matching step. I'd also recommend making ensemble numeric instead of character.

As for the choice between array and long formats, there's pros and cons of both. On the one hand, Array is more compact, and netCDF is elegantly self-documenting, but on the other hand netCDF is a pain in butt to actually work with in practice and a real barrier to inexperienced users. Long tabular (e.g. saved as CSV) is easier for most users but less compact. It might also be easier to link to the vector polygon file, as shapefiles are tabular not array-based (though may require a pivot-wider as I don't know if shapefiles allow you to repeat the primary key -- which could be ghastly if you're pivoting time, ensemble, and variable)

[dlebauer]

Got it in the prev. examples I was only using one time point and was just focusing on making the table EFI compliant. I decided to keep both array and table format for now to see what's easier to work with. My inclination is that tables will be easier.

And then I thought - surely it can't be harder to cram it into a NetCDF file than to write an EML file 😆?!

Here are updated arrays and table:

File sizes (2 variables, 8 time points, 20 ensemble members, and 98 sites)

• CSV: 2.0M
• list of arrays in RData file: 146K
• NetCDF 125K

CSV

> head(efi_long)                                                              
# A tibble: 6 × 5
  datetime            site             ensemble variable    prediction
  <dttm>              <chr>                <dbl> <chr>            <dbl>
1 2022-12-31 22:59:24 57ea98e608333b76         1 AGB            0.0720 
2 2022-12-31 22:59:24 57ea98e608333b76         1 TotSoilCarb    6.63   
3 2023-12-31 22:59:24 57ea98e608333b76         1 AGB            0.128  
4 2023-12-31 22:59:24 57ea98e608333b76         1 TotSoilCarb    6.58   
5 2022-12-31 22:59:24 729a6ef989c5a277         1 AGB            0.0010

Arrays

> head(ens_arrays)
$AGB
, , ensemble = 1

                     site
datetime              84389254458b3a58 9e3b1749179a4a68 729a6ef989c5a277
  2022-12-31 22:59:24         11.42190          0.74507          0.00103
  2023-12-31 22:59:24         10.78343          0.81111          0.04493
                     site
datetime              8c23aa17d5b211ed 57ea98e608333b76
  2022-12-31 22:59:24          1.02317          0.07197
  2023-12-31 22:59:24          0.93900          0.12845

, , ensemble = 2

                     site
datetime              84389254458b3a58 9e3b1749179a4a68 729a6ef989c5a277
  2022-12-31 22:59:24         13.05292          2.24812          0.30700
  2023-12-31 22:59:24         12.57112          2.34125          0.44772
                     site
datetime              8c23aa17d5b211ed 57ea98e608333b76
  2022-12-31 22:59:24          1.12851          0.47633
  2023-12-31 22:59:24          1.05095          0.52972

NetCDF file:

ncdump efi_forecast.nc 
netcdf efi_forecast {
dimensions:
        ntime = 8 ;
        site = 98 ;
        ensemble = 20 ;
variables:
        double ntime(ntime) ;
                ntime:units = "days since 1970-01-01 00:00:00" ;
                ntime:long_name = "Time middle averaging period" ;
                ntime:calendar = "standard" ;
        float time(ntime) ;
                time:units = "days since 1970-01-01 00:00:00" ;
                time:long_name = "Time dimension" ;
        int site(site) ;
                site:long_name = "Site ID" ;
        float lat(site) ;
                lat:units = "degrees_north" ;
                lat:long_name = "Latitude" ;
        float lon(site) ;
                lon:units = "degrees_east" ;
                lon:long_name = "Longitude" ;
        int ensemble(ensemble) ;
                ensemble:long_name = "ensemble member" ;
        float AGB(ensemble, site, ntime) ;
                AGB:units = "kg C m-2" ;
                AGB:long_name = "Total aboveground biomass" ;
        float TotSoilCarb(ensemble, site, ntime) ;
                TotSoilCarb:units = "kg C m-2" ;
                TotSoilCarb:long_name = "Total Soil Carbon" ;

// global attributes:
                :model_name = "SIPNET" ;
                :model_version = "v1.3" ;
                :forecast_iteration_id = "4812-4294967294" ;
                :forecast_time = "2025-03-12 07:48:52" ;
                :obs_flag = 0. ;
                :creation_date = "2025-03-13" ;
}

[mdietze]

EML: If you want to be a purist on EFI standard compliance, both netCDF and CSV end up with a EML because netCDF is only self documenting on file format and is not complete metadata. In the short term I wasn't worrying about generating metadata, thought it's a good point that it weould be much appreciated by end users

CSV vs netCDF file size: 100% agree that "raw" CSV is 10x bigger, and that difference will probably increase as the size of the output increases, but Carl B. would then ask, what's the size of the tabular if you then zip the CSV or output the tablular as parquet (e.g. arrow::write_parquet)? It will still be bigger, but I'm going to predict the difference won't be as dramatic. Parquet output can also be both compressed and tiled, with the latter making reading more efficient than a raw CSV (at least for large file sizes) though probably not any more efficient than netCDF also lot requiring you to load the entire object into local memory.

[dlebauer]

As an update - the code is in https://github.com/ccmmf/workflows/pull/2/files?diff=split#diff-481ce7e714b84da2b312b6fce6bd15c391e48a368531d74fcc065ffb10a43388

Plan is to move it to PEcAn.utils functions, perhaps a single pecan2efi function with the option of output type (and perhaps a non-exported helper function for each output type)

[dlebauer]

Here is a proposal for new functions to support multi-site ensemble runs. (Originally started as a GH issue, but figured there are enough unknowns and open questions that it would be helpful to get feedback here first.

Here is the draft:

Description

Create new functions that

1. read ensemble outputs from multi-site runs
2. convert these to EFI-standard output formats:

• long tidy CSV
• array RDS
• array NetCDF

Proposed Solution

A workflow might look like

ensemble_df <- read_multisite_ensemble(
  outdir = settings$outdir,
  site_ids = c("SITE1", "SITE2"),
  ens_ids = 1:10,
  variables = c("AGB", "TotSoilCarb"),
  start_year = 2021,
  end_year = 2025
)

ensemble_array <- ensemble2array(ensemble_df)
saveRDS(ensemble_array, "ensemble_output.rds")

ensemble_to_csv(ensemble_df, "ensemble_output.csv")

This is mostly based on https://github.com/ccmmf/workflows/blob/fb6aeb3b678ecabc89bc7af706aeafa05b799426/downscale/02_extract_sipnet_output.R#L85

New functions (please suggest changes to names and args)

• read_multisite_ensemble(outdir, site_ids, ens_ids, variables, start_year, end_year)
  • returns ensemble_df, a long format tidy data w/ cols datetime, site_id, ensemble, variable, prediction
• ensemble2array(ensemble_df)
• ensemble_nc2efi_csv(ensemble_df, outfile, metadata = FALSE) writes CSV
• ensemble_nc2efi_array(ensemble_array, outfile, ...) writes RDS
• ensemble_nc2efi_nc(ensemble_array, outfile, ...) writes netcdf
• generate_efi_metadata(ensemble_df, ...) if metadata = TRUE, this helper function could write out EFI+EML - compliant xml or json.

I've written code that produces the arrays and netcdf files, but have been using the tabular CSV files in my workflows.

Considerations

• If we can parse ensemble and site ids from the directory structure, outdir would be the only required argument to read_multisite_ensemble.
• Should we consider refactoring existing read.ensemble and read.ensemble.ts functions? The outputs are a list (still easy to use), but neither of the ensemble x site implementations currently use these.
• How to handle time slices? What should we use for time slicing? Current read.ensemble functions have start and end date. We could also allow 'all', 'last', etc, or have a wrapper function that calls multiple times for a vector of dates (each call would have start_date == end_date).

Additional Context

• prototypes developed for @ccmmf project:
  • the script I wrote for reading and writing ensemble outputs https://github.com/ccmmf/workflows/blob/fb6aeb3b678ecabc89bc7af706aeafa05b799426/downscale/02_extract_sipnet_output.R#L85
  • read_ncdir written by @infotroph https://github.com/ccmmf/workflows/blob/22e34f246ebc8f83fea1e541772272e9fe9243fc/phase_1b_ca_woody/irri_compare_20250319.Rmd#L28
• existing functions
  • read.ensemble https://github.com/PecanProject/pecan/blob/develop/modules/uncertainty/R/ensemble.R
  • read.ensemble.ts https://github.com/PecanProject/pecan/blob/develop/modules/uncertainty/R/run.ensemble.analysis.R