Merge pull request #4039 from Tejas7007/GH-4038-ilamb-multimodel

Add multi-model CMIP6/TRENDY benchmarking to inst/ilamb (#4038)

Author: mdietze

CHANGELOG.md

@@ -13,6 +13,7 @@ For more information about this file see also [Keep a Changelog](http://keepacha
     with units attached as an attribute.
 - New package `PEcAn.RothC` runs the RothC soil carbon model.
 - Added `inst/ilamb/` pipeline in PEcAn.benchmark to convert downscaled SDA reanalysis GeoTIFFs into ILAMB-compatible CF netCDF for carbon-cycle benchmarking (#4019).
+- Added multi-model benchmarking scripts to `inst/ilamb/` in PEcAn.benchmark: build the CMIP6 and TRENDY ensembles, score the PEcAn reanalysis and its individual SDA members against observational benchmarks over two evaluation windows, and summarize ensemble spread (#4038).
 - Added PEcAn.PEPRMT model, including a demo run with example data
 - Add `format_try_for_ma()` and `try_trait_mapping()` to `PEcAn.data.remote` to convert trait data from the external TRY database into the tabular format required by the PEcAn meta-analysis module (#3717).
 - Add function `qsub_sda()` for submitting SDA batch jobs by splitting a large number of sites into multiple small groups of sites (#3634).

modules/benchmark/inst/ilamb/README.md

@@ -103,3 +103,196 @@ loading. Set `ILAMB_OUTPUT_DIR` to point the tests at your output directory.
 - A known structural discontinuity exists in the underlying LandTrendr input
   around 2017-2018 (see the ORNL DAAC documentation); it is preserved as-is in
   the converted output rather than adjusted here.
+
+---
+
+# Multi-Model Benchmarking (CMIP6 and TRENDY)
+
+Building on the conversion pipeline above, these scripts benchmark the PEcAn
+reanalysis against two major model intercomparison ensembles, CMIP6 and TRENDY
+(Global Carbon Budget), using ILAMB, and score the individual PEcAn ensemble
+members so that PEcAn's skill spread can be compared directly with the model
+ensembles' spread.
+
+Comparison is over North America for three state variables with established
+observational benchmarks: vegetation carbon (`cVeg`), soil carbon (`cSoil`),
+and leaf area index (`lai`). Soil moisture is supported by the conversion step
+but is not yet included in the multi-model scoring (see Scope below).
+
+## Pipeline
+
+Run in order; each step writes inputs for the next. The first script is the
+conversion tool documented above.
+
+| Step | Script | Purpose |
+|------|--------|---------|
+| 1 | `convert_geotiff_to_ilamb.py` | PEcAn GeoTIFF ensemble maps to CF netCDF |
+| 2 | `build_cmip6_ensemble.py` | Download + regrid CMIP6 historical fields |
+| 3 | `build_cmip6_ssp.py` | Splice CMIP6 historical + ssp245 for a longer record |
+| 4 | `build_trendy_ensemble.py` | Download + regrid TRENDY (GCB) fields |
+| 5 | `build_window_ensembles.py` | Slice all models to an evaluation window; build ensemble means |
+| 6 | `build_pecan_members.py` | Per-member PEcAn fields for spread analysis |
+| 7 | `make_spread_figures.py` | Spread figure from ILAMB scores |
+
+All model fields are regridded onto the same 0.5 degree North American grid the
+conversion step produces, so PEcAn, the observational benchmarks, and every
+model share one grid. LAI is reduced to its July value throughout, to match the
+July snapshot of the PEcAn product. The TRENDY ensemble mean is built only from
+members that pass a per-variable physical-plausibility screen; the members used
+are recorded in `trendy_ensemble_manifest.json`.
+
+## Evaluation windows
+
+The PEcAn reanalysis and the models are compared over two windows, which answer
+different questions:
+
+- **2012-2014**: a mean-state snapshot using the full set of CMIP6 historical
+  models (25), the broadest representative model sample.
+- **2015-2023**: a longer record for interannual variability and trends, using
+  the CMIP6 models that provide a continuous historical + ssp245 land-carbon
+  record (14).
+
+Reporting both is deliberate: the longer window supports variability and trend
+analysis that three years cannot, while the 2012-2014 window retains the full
+model sample. The two windows also differ in CMIP6 composition, and that
+difference is itself informative; see the note on soil carbon below.
+
+## ILAMB configuration
+
+`pecan_ilamb.cfg` defines three confrontations:
+
+| Variable | Benchmark | Notes |
+|----------|-----------|-------|
+| Biomass (cVeg) | Xu & Saatchi 2021 | RMSE skipped; mass-weighted |
+| Leaf Area Index | GIMMS LAI4g | seasonal-cycle and RMSE skipped |
+| Soil Carbon | HWSD2 | RMSE skipped; mass-weighted |
+
+The configuration applies no time or region restriction itself; the evaluation
+window is set entirely by the windowed input files (step 5), and the North
+American extent is set by the grid.
+
+## Results
+
+ILAMB overall scores (0-1, higher is a closer match to the benchmark). PEcAn is
+the reanalysis ensemble mean; CMIP6 and TRENDY are ensemble means.
+
+**2012-2014 (25 CMIP6 models)**
+
+| Variable | PEcAn | CMIP6 | TRENDY |
+|----------|-------|-------|--------|
+| Biomass | 0.483 | 0.477 | 0.410 |
+| Leaf Area Index | 0.510 | 0.448 | 0.453 |
+| Soil Carbon | 0.625 | 0.634 | 0.478 |
+
+**2015-2023 (14 CMIP6 models)**
+
+| Variable | PEcAn | CMIP6 | TRENDY |
+|----------|-------|-------|--------|
+| Biomass | 0.464 | 0.473 | 0.409 |
+| Leaf Area Index | 0.513 | 0.448 | 0.454 |
+| Soil Carbon | 0.648 | 0.715 | 0.478 |
+
+PEcAn scores above both model ensembles on leaf area index in every window, and
+above the TRENDY ensemble on all three variables. On biomass and soil carbon it
+is comparable to the CMIP6 ensemble.
+
+The higher CMIP6 soil-carbon score in the 2015-2023 window (0.715) is a
+composition effect, not a change in the benchmark: the 14 models with a
+continuous ssp245 record happen to be stronger soil-carbon performers. The full
+25-model window (0.634) is the representative figure, where PEcAn and CMIP6 are
+close.
+
+## Ensemble spread
+
+Scoring the 100 individual PEcAn members alongside the individual CMIP6 and
+TRENDY models shows that the PEcAn members vary far less in skill than the
+models do. Standard deviation of member scores:
+
+| Variable | PEcAn (100) | CMIP6 | TRENDY |
+|----------|-------------|-------|--------|
+| Biomass (2012-2014) | 0.003 | 0.057 | 0.086 |
+| LAI (2012-2014) | 0.004 | 0.082 | 0.145 |
+| Soil Carbon (2012-2014) | 0.003 | 0.166 | 0.137 |
+| Biomass (2015-2023) | 0.014 | 0.040 | 0.086 |
+| LAI (2015-2023) | 0.002 | 0.085 | 0.147 |
+| Soil Carbon (2015-2023) | 0.001 | 0.146 | 0.137 |
+
+The narrow PEcAn spread reflects skill, not identical members: the members
+differ in their spatial carbon fields (for soil carbon, on the order of one
+percent of the field, member to member) yet match the benchmark about equally
+well. Whether that spread is appropriately sized relative to the error, that is,
+whether the ensemble is well-calibrated, is a separate question, addressed by
+probabilistic scoring in a later contribution.
+
+`fig_spread_clouds.png` shows the per-member scores for all three ensembles
+across both windows.
+
+## Scope and caveats
+
+- **North America only.** The comparison uses the North American PEcAn grid.
+- **Biomass pool mismatch.** PEcAn `cVeg` is above-ground wood carbon, while the
+  Xu & Saatchi benchmark is total live biomass (including roots). PEcAn is
+  therefore expected to read somewhat low on biomass by construction.
+- **Soil-carbon depth mismatch.** PEcAn soil carbon integrates 0-200 cm, deeper
+  than the benchmark, so PEcAn reads higher in deep-carbon regions.
+- **Soil moisture deferred.** The conversion step handles soil moisture, but
+  multi-model soil-moisture scoring is not yet included.
+- **Probabilistic scoring is future work.** The spread results above motivate
+  ensemble calibration metrics (rank histograms, reliability), which are a
+  planned follow-on rather than part of this contribution.
+
+## Reproducing the figure
+
+`make_spread_figures.py` reads ILAMB `scores.csv` output for the two windows.
+Its paths refer to the analysis working tree where the scoring runs were
+produced; point them at your own ILAMB build directories to regenerate the
+figure.
+
+## Data sources
+
+- CMIP6: ESGF, via `intake-esgf` (historical and ssp245, r1i1p1f1).
+- TRENDY: Global Carbon Budget 2024, public download index.
+- Benchmarks: Xu & Saatchi 2021 (biomass), GIMMS LAI4g (leaf area index),
+  HWSD2 (soil carbon).
+
+## Customizing the benchmark
+
+The benchmark is designed to stay usable as models and datasets evolve. Three
+common changes: updating or adding models, adding benchmark datasets, and
+adjusting the scoring weights. The per-field meaning of every configuration
+option is documented inline in `pecan_ilamb.cfg`.
+
+### Updating or adding a model
+
+Each model is just a directory under the model root (for example
+`ilamb_models_2012_2014/`), containing `cVeg.nc`, `cSoil.nc`, and `lai.nc` on
+the analysis grid. ILAMB scores every model directory it finds there, so:
+
+- To add a model, place its directory (with the three variable files, regridded
+  to the analysis grid) under the model root and re-run `ilamb-run`.
+- To update a model to a newer version (for example a new CMIP6 or TRENDY
+  release), regenerate its files with the relevant build script and overwrite
+  that model's directory. The build scripts take the source and output paths as
+  arguments, so pointing them at a new data release does not require code edits.
+
+The ensemble means (CMIP6, TRENDY, PEcAn) are rebuilt from the individual
+members by `build_window_ensembles.py`, so re-run it after changing the set of
+member models so the means reflect the new set.
+
+### Adding a benchmark dataset
+
+Benchmarks are defined in `pecan_ilamb.cfg`. To add one for an existing
+variable, copy the dataset block under the relevant `[h2: ...]` variable, rename
+it, and point its `source` at the new benchmark netCDF (placed under
+`ILAMB_ROOT/DATA/...`). When more than one benchmark is listed under a variable,
+their `weight` values set how they combine. To score a new variable entirely,
+add a new `[h2: ...]` block and make sure the model files contain that variable.
+
+### Adjusting weights
+
+Scoring weights live in `pecan_ilamb.cfg` at two levels: the `weight` on each
+`[h2: ...]` variable sets how much that variable counts toward the overall
+score relative to the other variables, and the `weight` on each dataset block
+sets how much that dataset counts within its variable. Both are relative
+weights, not percentages. See the inline comments in the config for the full
+per-field reference.

modules/benchmark/inst/ilamb/build_cmip6_ensemble.py

@@ -0,0 +1,85 @@
+"""
+Download and regrid CMIP6 historical land-carbon fields for benchmarking.
+
+Retrieves cVeg, cSoil, and lai for a set of CMIP6 models (one realization each,
+r1i1p1f1) from the ESGF archive via intake-esgf, then regrids every field onto
+the PEcAn analysis grid so all models, observations, and the reanalysis share a
+common grid for ILAMB scoring. Downloads are cached, so re-runs resume rather
+than refetch.
+
+LAI is reduced to the July value to match the growing-season convention used on
+the reanalysis side; carbon pools use the full monthly series. Longitudes are
+converted to the [-180, 180] convention where needed before interpolation.
+
+Output: one directory per model, each containing cVeg.nc / cSoil.nc / lai.nc.
+The full historical span is written here; restriction to an evaluation window
+is handled downstream by build_window_ensembles.py.
+
+Paths default to the values used for this analysis; override with the flags
+below to run elsewhere.
+"""
+import argparse
+import os
+import glob
+import xarray as xr
+import intake_esgf
+from intake_esgf import ESGFCatalog
+
+MODELS = ["ACCESS-ESM1-5","AWI-ESM-1-1-LR","BCC-CSM2-MR","BCC-ESM1","CESM2",
+          "CESM2-FV2","CESM2-WACCM","CESM2-WACCM-FV2","CMCC-CM2-SR5","CMCC-ESM2",
+          "CanESM5","CanESM5-1","EC-Earth3-CC","EC-Earth3-Veg","EC-Earth3-Veg-LR",
+          "GFDL-ESM4","IPSL-CM5A2-INCA","IPSL-CM6A-LR","IPSL-CM6A-LR-INCA","KIOST-ESM",
+          "MPI-ESM-1-2-HAM","MPI-ESM1-2-LR","NorESM2-LM","NorESM2-MM","SAM0-UNICON","TaiESM1"]
+
+ap = argparse.ArgumentParser(description=__doc__)
+ap.add_argument("--cache",   default="/projectnb/dietzelab/tdahiya/cmip6_cache")
+ap.add_argument("--outroot", default="/projectnb/dietzelab/tdahiya/ilamb_models_multi")
+ap.add_argument("--pecan-grid", dest="pec",
+                default="/projectnb/dietzelab/tdahiya/ilamb_models/PEcAn/cVeg.nc",
+                help="A PEcAn netCDF whose lat/lon define the target grid.")
+args = ap.parse_args()
+CACHE, OUTROOT, PEC = args.cache, args.outroot, args.pec
+
+intake_esgf.conf.set(local_cache=[CACHE])
+
+print("=== Download (resumes from cache) ===", flush=True)
+cat = ESGFCatalog()
+cat.search(project="CMIP6", experiment_id="historical", source_id=MODELS,
+           variable_id=["cVeg","cSoil","lai"], table_id=["Lmon","Emon"], variant_label="r1i1p1f1")
+cat.to_path_dict()
+print("Download done.\n=== Processing ===", flush=True)
+
+pec = xr.open_dataset(PEC, decode_times=False)
+tgt_lat, tgt_lon = pec.lat.values, pec.lon.values
+
+def proc(model, var):
+    files = sorted(glob.glob(f"{CACHE}/**/{model}/**/{var}_*.nc", recursive=True))
+    if not files: return None
+    da = xr.open_mfdataset(files, decode_times=True, combine="by_coords")[var]
+    units = da.attrs.get("units","")
+    da = da.sel(time=slice("2000-01-01","2014-12-31"))
+    if var == "lai":
+        da = da.sel(time=da.time.dt.month == 7)
+    if float(da.lon.max()) > 180 and tgt_lon.min() < 0:
+        da = da.assign_coords(lon=(((da.lon+180)%360)-180)).sortby("lon")
+    da = da.interp(lat=tgt_lat, lon=tgt_lon, method="linear").load()
+    da.attrs["units"]=units; da.name=var
+    da["time"].attrs.pop("bounds",None); da["time"].encoding.pop("bounds",None)
+    da["lat"].attrs.update({"units":"degrees_north","standard_name":"latitude"})
+    da["lon"].attrs.update({"units":"degrees_east","standard_name":"longitude"})
+    return da
+
+ok=[]
+for m in MODELS:
+    try:
+        od=f"{OUTROOT}/{m}"; os.makedirs(od, exist_ok=True); n=0
+        for var in ["cVeg","cSoil","lai"]:
+            da=proc(m,var)
+            if da is not None:
+                xr.Dataset({var:da}).to_netcdf(f"{od}/{var}.nc"); n+=1
+        print(f"  {m}: {'OK' if n==3 else 'PARTIAL'} ({n}/3)", flush=True)
+        if n==3: ok.append(m)
+    except Exception as e:
+        print(f"  {m}: FAILED - {str(e)[:90]}", flush=True)
+
+print(f"\nDone. {len(ok)}/{len(MODELS)} fully processed:", ok)

modules/benchmark/inst/ilamb/build_cmip6_ssp.py

@@ -0,0 +1,89 @@
+"""
+Build a continuous CMIP6 record by splicing historical and ssp245 runs.
+
+Same models, variables, and regridding as the historical builder, but joins the
+end of each model's historical run to the start of its ssp245 scenario so the
+record extends past 2014. Only models that provide both pieces with land-carbon
+output are kept; the rest are reported as having no usable scenario run.
+
+This spliced record is the source for the longer evaluation window, where
+interannual variability and trends can be assessed. LAI is again reduced to July.
+
+Output: one directory per model with cVeg.nc / cSoil.nc / lai.nc.
+
+Paths default to the values used for this analysis; override with the flags
+below to run elsewhere.
+"""
+import argparse
+import os
+import glob
+import xarray as xr
+import intake_esgf
+from intake_esgf import ESGFCatalog
+
+MODELS = ["ACCESS-ESM1-5","AWI-ESM-1-1-LR","BCC-CSM2-MR","BCC-ESM1","CESM2",
+          "CESM2-FV2","CESM2-WACCM","CESM2-WACCM-FV2","CMCC-CM2-SR5","CMCC-ESM2",
+          "CanESM5","CanESM5-1","EC-Earth3-CC","EC-Earth3-Veg","EC-Earth3-Veg-LR",
+          "GFDL-ESM4","IPSL-CM5A2-INCA","IPSL-CM6A-LR","IPSL-CM6A-LR-INCA","KIOST-ESM",
+          "MPI-ESM-1-2-HAM","MPI-ESM1-2-LR","NorESM2-LM","NorESM2-MM","SAM0-UNICON","TaiESM1"]
+
+ap = argparse.ArgumentParser(description=__doc__)
+ap.add_argument("--cache",   default="/projectnb/dietzelab/tdahiya/cmip6_cache")
+ap.add_argument("--outroot", default="/projectnb/dietzelab/tdahiya/ilamb_models_2012_2023")
+ap.add_argument("--pecan-grid", dest="pec",
+                default="/projectnb/dietzelab/tdahiya/ilamb_models/PEcAn/cVeg.nc",
+                help="A PEcAn netCDF whose lat/lon define the target grid.")
+args = ap.parse_args()
+CACHE, OUTROOT, PEC = args.cache, args.outroot, args.pec
+
+intake_esgf.conf.set(local_cache=[CACHE])
+
+print("=== Download ssp245 (historical already cached) ===", flush=True)
+cat = ESGFCatalog()
+cat.search(project="CMIP6", experiment_id="ssp245", source_id=MODELS,
+           variable_id=["cVeg","cSoil","lai"], table_id=["Lmon","Emon"], variant_label="r1i1p1f1")
+cat.to_path_dict()
+print("Download done.\n=== Splice hist(2012-2014)+ssp245(2015-2023) ===", flush=True)
+
+pec = xr.open_dataset(PEC, decode_times=False)
+tgt_lat, tgt_lon = pec.lat.values, pec.lon.values
+
+def load_exp(model, var, exp, t0, t1):
+    files = sorted(glob.glob(f"{CACHE}/**/{model}/**/{var}_*.nc", recursive=True))
+    files = [f for f in files if exp in os.path.basename(f)]
+    if not files: return None
+    da = xr.open_mfdataset(files, decode_times=True, combine="by_coords")[var]
+    return da.sel(time=slice(t0, t1))
+
+def proc(model, var):
+    hist = load_exp(model, var, "historical", "2012-01-01", "2014-12-31")
+    ssp  = load_exp(model, var, "ssp245",     "2015-01-01", "2023-12-31")
+    if hist is None or ssp is None: return None
+    units = hist.attrs.get("units","")
+    da = xr.concat([hist, ssp], dim="time").sortby("time")
+    if var == "lai":
+        da = da.sel(time=da.time.dt.month == 7)
+    if float(da.lon.max()) > 180 and tgt_lon.min() < 0:
+        da = da.assign_coords(lon=(((da.lon+180)%360)-180)).sortby("lon")
+    da = da.interp(lat=tgt_lat, lon=tgt_lon, method="linear").load()
+    da.attrs["units"]=units; da.name=var
+    da["time"].attrs.pop("bounds",None); da["time"].encoding.pop("bounds",None)
+    da["lat"].attrs.update({"units":"degrees_north","standard_name":"latitude"})
+    da["lon"].attrs.update({"units":"degrees_east","standard_name":"longitude"})
+    return da
+
+ok=[]; no_ssp=[]
+for m in MODELS:
+    try:
+        od=f"{OUTROOT}/{m}"; os.makedirs(od, exist_ok=True); n=0
+        for var in ["cVeg","cSoil","lai"]:
+            da=proc(m,var)
+            if da is not None:
+                xr.Dataset({var:da}).to_netcdf(f"{od}/{var}.nc"); n+=1
+        print(f"  {m}: {'OK' if n==3 else 'PARTIAL'} ({n}/3)", flush=True)
+        if n==3: ok.append(m)
+        elif n==0: no_ssp.append(m)
+    except Exception as e:
+        print(f"  {m}: FAILED - {str(e)[:90]}", flush=True)
+print(f"\nDone. {len(ok)}/{len(MODELS)} full splice: {ok}", flush=True)
+print(f"No usable ssp245 (likely no scenario run): {no_ssp}", flush=True)