From 9429756b6ae07e80f0b8486437aa5b2c144bf6d4 Mon Sep 17 00:00:00 2001 From: mvertens Date: Fri, 2 Oct 2026 13:36:47 +0200 Subject: [PATCH 01/16] Skip empty history files before GenTS metadata pull --- scripts/gen_timeseries.py | 75 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 74 insertions(+), 1 deletion(-) diff --git a/scripts/gen_timeseries.py b/scripts/gen_timeseries.py index 3fef067..252137d 100755 --- a/scripts/gen_timeseries.py +++ b/scripts/gen_timeseries.py @@ -14,6 +14,7 @@ import glob import argparse import logging +from concurrent.futures import ProcessPoolExecutor # Determine local directory path: _LOCAL_PATH = os.path.dirname(os.path.abspath(__file__)) @@ -21,11 +22,73 @@ from pathlib import Path # Time series generation imports +import netCDF4 from gents.hfcollection import HFCollection from gents.timeseries import TSCollection from cmip7_prep.include_patterns import all_include_patterns +# ++++++++++++++++++++++++++++++ +# Unusable history file detection +# ++++++++++++++++++++++++++++++ + + +def time_record_count(path): + """ + Returns the number of time records in a history file: 0 if it holds none, + None if it has no time dimension, or -1 if it will not open. + """ + try: + with netCDF4.Dataset(path) as ds: + time_dim = ds.dimensions.get("time") + if time_dim is None: + return None + return time_dim.size + except OSError: + return -1 + + +def find_unusable_files(inputdir, include_patterns, workers, logger): + """ + Returns the history files GenTS cannot derive time bounds from: those with + zero time records, and those that will not open. + + GenTS takes a minimum over each file's time coordinate to order the files + into a series, so a file with no records raises a zero-size reduction in + gents.meta. Files like that are written when a run segment opens a history + stream but ends before any sample reaches it. They hold no data, so they + are dropped here instead of being left for GenTS to warn about one logged + traceback at a time. + """ + candidates = sorted( + { + path + for pattern in include_patterns + for path in glob.glob(os.path.join(inputdir, pattern)) + } + ) + if not candidates: + return [] + + unusable = [] + with ProcessPoolExecutor(max_workers=min(workers, len(candidates))) as executor: + for path, count in zip(candidates, executor.map(time_record_count, candidates)): + if count == 0: + logger.warning("Skipping %s: zero time records", path) + unusable.append(path) + elif count == -1: + logger.warning("Skipping %s: cannot be opened", path) + unusable.append(path) + + if unusable: + logger.warning( + "Excluding %d of %d history file(s) from time series generation", + len(unusable), + len(candidates), + ) + return unusable + + # ++++++++++++++++++++++++++++++ # Input argument parser function # ++++++++++++++++++++++++++++++ @@ -208,11 +271,17 @@ def main(): sys.exit(0) logger.info(f"include patterns are {include_patterns}") +<<<<<<< HEAD varlist = ( [variable.strip() for variable in args.varlist.split(",") if variable.strip()] if args.varlist else None ) +======= + # Drop files GenTS cannot read time bounds from before they reach it + unusable_files = find_unusable_files(inputdir, include_patterns, workers, logger) + +>>>>>>> 92bd470 (Skip empty history files before GenTS metadata pull) # Determine how time series will be created if not args.years_spec: @@ -220,6 +289,8 @@ def main(): logger.info("Starting hf_collection") hf_collection = HFCollection(inputdir, num_processes=workers) hf_collection = hf_collection.include(include_patterns) + if unusable_files: + hf_collection = hf_collection.exclude(unusable_files) logger.info("Finished hf_collection") # Create base TSCollection @@ -244,12 +315,14 @@ def main(): logger.info("Year increment for time series generation is %s", nyears) hf_collection = HFCollection(inputdir, num_processes=workers) + if unusable_files: + hf_collection = hf_collection.exclude(unusable_files) for include_pattern in include_patterns: logger.info("Processing files with pattern: %s", include_pattern) for year in range(year_first, year_last + 1, nyears): logger.info(f"Processing from year {year} to year {year+nyears-1}") - hfp_collection = hf_collection.include_patterns([include_pattern]) + hfp_collection = hf_collection.include([include_pattern]) hfp_collection = hfp_collection.include_years(year, year + nyears - 1) logger.info(f"files to process for year {year} are") From 105cb0608b9d381192c3b71e75c1a11d751b8b87 Mon Sep 17 00:00:00 2001 From: mvertens Date: Sat, 3 Oct 2026 10:46:04 +0200 Subject: [PATCH 02/16] Fix license_id, numpy-dependent doctest, and test collection --- .pylintrc | 4 +- INSTALL.md | 59 ++++++++++++++++++++++++++++ data/cmor_dataset.json | 2 +- pytest.ini | 2 +- scripts/gen_timeseries.py | 75 +----------------------------------- src/cmip7_prep/cmor_utils.py | 4 +- 6 files changed, 66 insertions(+), 80 deletions(-) diff --git a/.pylintrc b/.pylintrc index 4948095..93a61a9 100644 --- a/.pylintrc +++ b/.pylintrc @@ -1,6 +1,6 @@ [MASTER] ignore-paths=build,dist,.venv,.eggs,scripts -extension-pkg-allow-list=cmor +extension-pkg-allow-list=cmor,netCDF4 init-hook='import sys, os; src = os.path.join(os.getcwd(), "src"); sys.path.insert(0, src) if src not in sys.path else None' [MESSAGES CONTROL] @@ -17,5 +17,5 @@ disable= max-line-length=100 [TYPECHECK] -ignored-modules=xarray,xesmf,cmor,numpy,geocat.comp,yaml,cftime,pytest,dask +ignored-modules=xarray,xesmf,cmor,numpy,geocat.comp,yaml,cftime,pytest,dask,netCDF4 ignored-classes=cmor diff --git a/INSTALL.md b/INSTALL.md index f2301d6..7bc32e1 100644 --- a/INSTALL.md +++ b/INSTALL.md @@ -49,6 +49,40 @@ pip install gents dulwich cmip7-data-request-api These three are not on conda-forge, so they must come from pip. +Note that `pip install ` does nothing when the package is already present. +To move to a newer release, ask for it explicitly: + +```bash +pip install -U gents dulwich cmip7-data-request-api +``` + +## Local development environment (macOS) + +The steps above create the shared analysis environment on NIRD. For running the +test suite and editing the code on a laptop, create a private environment +instead. Everything the tests need, including CMOR and ESMF, has an +`osx-arm64` build on conda-forge, so no part of the stack has to be skipped. + +```bash +conda create -n cmip7-dev -c conda-forge python=3.12 cmor=3.15 \ + numpy xarray netcdf4 cftime pyyaml pytest pytest-cov \ + xesmf=0.8.7 esmpy=8.9.0 geocat-comp=2025.10.01 +``` + +```bash +conda activate cmip7-dev +pip install -e . --no-deps +pip install gents +``` + +```bash +pytest -q +``` + +Use **Miniforge** rather than Anaconda for this (see the BLAS entry under +Troubleshooting). `dulwich` and `cmip7-data-request-api` are only needed for +the data-request and CV tooling, so a test-only environment can leave them out. + ## CMIP7 CMOR tables This recipe only provides the CMOR **library** (3.15). The matching **tables** @@ -67,3 +101,28 @@ labels). Keep that checkout up to date, or CMOR will reject otherwise-valid runs - **`license_id "..." could not be found` / `grid_label "..." is invalid`:** the `cmip7-cmor-tables` checkout is out of date. Update it to the version matching CMOR 3.15 / the CMIP7 controlled vocabulary. +- **`zero-size array to reduction operation minimum which has no identity`, + logged with a traceback while pulling metadata:** one or more history files + hold no time records, which happens when a run segment opens a history stream + but ends before a sample is written. GenTS takes a minimum over each file's + time coordinate, so an empty file fails. This is **not fatal** -- + `HFCollection.pull_metadata` defaults to `raise_errors=False`, logs the + failure with `exc_info=True` and drops the file, so the run continues. + `gen_timeseries.py` screens such files out up front so they appear as one + warning each instead of a traceback. +- **`Error importing numpy: you should not try to import numpy from its source + directory`:** misleading message; numpy prints it whenever its compiled + extension fails to load for any reason. Get the real cause with + ```bash + python -c "import traceback + try: import numpy + except ImportError as e: traceback.print_exception(e.__cause__)" + ``` + A mixed-channel Anaconda install can leave `libcblas`/`libblas` symlinked to + an OpenBLAS build that was never installed, in which case the cause names a + missing `.dylib`. A conda-forge-only base (Miniforge) avoids this; repairing + a broken one needs `conda install --force-reinstall openblas blas`. +- **`environment.yml` disagrees with these instructions:** it is an export of + the older `cmordev_env`, not the `cmordev_env_312` built above, so its pins + (including `cmor` and `gents`) are not the versions this recipe installs. + Treat the steps here as authoritative. diff --git a/data/cmor_dataset.json b/data/cmor_dataset.json index 449b1cc..fcd2666 100644 --- a/data/cmor_dataset.json +++ b/data/cmor_dataset.json @@ -31,7 +31,7 @@ "sub_experiment": "none", "archive_id": "WCRP", "area_label" : "u", - "license_id": "CC-BY-4-0", + "license_id": "CC-BY-4.0", "realization_index": "r1", "initialization_index": "i1", "physics_index": "p1", diff --git a/pytest.ini b/pytest.ini index e58dcb8..d73e499 100644 --- a/pytest.ini +++ b/pytest.ini @@ -1,4 +1,4 @@ [pytest] addopts = --doctest-modules -norecursedirs = scripts +norecursedirs = scripts cmip7-cmor-tables pythonpath = src . diff --git a/scripts/gen_timeseries.py b/scripts/gen_timeseries.py index 252137d..3fef067 100755 --- a/scripts/gen_timeseries.py +++ b/scripts/gen_timeseries.py @@ -14,7 +14,6 @@ import glob import argparse import logging -from concurrent.futures import ProcessPoolExecutor # Determine local directory path: _LOCAL_PATH = os.path.dirname(os.path.abspath(__file__)) @@ -22,73 +21,11 @@ from pathlib import Path # Time series generation imports -import netCDF4 from gents.hfcollection import HFCollection from gents.timeseries import TSCollection from cmip7_prep.include_patterns import all_include_patterns -# ++++++++++++++++++++++++++++++ -# Unusable history file detection -# ++++++++++++++++++++++++++++++ - - -def time_record_count(path): - """ - Returns the number of time records in a history file: 0 if it holds none, - None if it has no time dimension, or -1 if it will not open. - """ - try: - with netCDF4.Dataset(path) as ds: - time_dim = ds.dimensions.get("time") - if time_dim is None: - return None - return time_dim.size - except OSError: - return -1 - - -def find_unusable_files(inputdir, include_patterns, workers, logger): - """ - Returns the history files GenTS cannot derive time bounds from: those with - zero time records, and those that will not open. - - GenTS takes a minimum over each file's time coordinate to order the files - into a series, so a file with no records raises a zero-size reduction in - gents.meta. Files like that are written when a run segment opens a history - stream but ends before any sample reaches it. They hold no data, so they - are dropped here instead of being left for GenTS to warn about one logged - traceback at a time. - """ - candidates = sorted( - { - path - for pattern in include_patterns - for path in glob.glob(os.path.join(inputdir, pattern)) - } - ) - if not candidates: - return [] - - unusable = [] - with ProcessPoolExecutor(max_workers=min(workers, len(candidates))) as executor: - for path, count in zip(candidates, executor.map(time_record_count, candidates)): - if count == 0: - logger.warning("Skipping %s: zero time records", path) - unusable.append(path) - elif count == -1: - logger.warning("Skipping %s: cannot be opened", path) - unusable.append(path) - - if unusable: - logger.warning( - "Excluding %d of %d history file(s) from time series generation", - len(unusable), - len(candidates), - ) - return unusable - - # ++++++++++++++++++++++++++++++ # Input argument parser function # ++++++++++++++++++++++++++++++ @@ -271,17 +208,11 @@ def main(): sys.exit(0) logger.info(f"include patterns are {include_patterns}") -<<<<<<< HEAD varlist = ( [variable.strip() for variable in args.varlist.split(",") if variable.strip()] if args.varlist else None ) -======= - # Drop files GenTS cannot read time bounds from before they reach it - unusable_files = find_unusable_files(inputdir, include_patterns, workers, logger) - ->>>>>>> 92bd470 (Skip empty history files before GenTS metadata pull) # Determine how time series will be created if not args.years_spec: @@ -289,8 +220,6 @@ def main(): logger.info("Starting hf_collection") hf_collection = HFCollection(inputdir, num_processes=workers) hf_collection = hf_collection.include(include_patterns) - if unusable_files: - hf_collection = hf_collection.exclude(unusable_files) logger.info("Finished hf_collection") # Create base TSCollection @@ -315,14 +244,12 @@ def main(): logger.info("Year increment for time series generation is %s", nyears) hf_collection = HFCollection(inputdir, num_processes=workers) - if unusable_files: - hf_collection = hf_collection.exclude(unusable_files) for include_pattern in include_patterns: logger.info("Processing files with pattern: %s", include_pattern) for year in range(year_first, year_last + 1, nyears): logger.info(f"Processing from year {year} to year {year+nyears-1}") - hfp_collection = hf_collection.include([include_pattern]) + hfp_collection = hf_collection.include_patterns([include_pattern]) hfp_collection = hfp_collection.include_years(year, year + nyears - 1) logger.info(f"files to process for year {year} are") diff --git a/src/cmip7_prep/cmor_utils.py b/src/cmip7_prep/cmor_utils.py index dd99c03..958cca6 100644 --- a/src/cmip7_prep/cmor_utils.py +++ b/src/cmip7_prep/cmor_utils.py @@ -204,7 +204,7 @@ def roll_for_monotonic_with_bounds(lon, lon_bnds): >>> lon array([ 0., 90., 180., 270.]) >>> shift - np.int64(-2) + -2 >>> lon2 = np.array([0., 90., 180., 270.]) >>> lon2_bnds = np.array([[-45., 45.], [45., 135.], [135., 225.], [225., 315.]]) >>> _, _, shift2 = roll_for_monotonic_with_bounds(lon2, lon2_bnds) @@ -217,7 +217,7 @@ def roll_for_monotonic_with_bounds(lon, lon_bnds): shift = k[0] + 1 lon = np.roll(lon, -shift) lon_bnds = np.roll(lon_bnds, -shift, axis=0) - return lon, lon_bnds, -shift + return lon, lon_bnds, int(-shift) return lon, lon_bnds, 0 From cf3ef296964597f15e47a3688fd3c736746c5abf Mon Sep 17 00:00:00 2001 From: mvertens Date: Sat, 3 Oct 2026 21:21:26 +0200 Subject: [PATCH 03/16] Add synthetic timeseries tests for every realm and frequency --- .pylintrc | 4 +- tests/test_timeseries_synthetic.py | 218 +++++++++++++++++++++++++++++ 2 files changed, 220 insertions(+), 2 deletions(-) create mode 100644 tests/test_timeseries_synthetic.py diff --git a/.pylintrc b/.pylintrc index 93a61a9..965e631 100644 --- a/.pylintrc +++ b/.pylintrc @@ -1,6 +1,6 @@ [MASTER] ignore-paths=build,dist,.venv,.eggs,scripts -extension-pkg-allow-list=cmor,netCDF4 +extension-pkg-allow-list=cmor init-hook='import sys, os; src = os.path.join(os.getcwd(), "src"); sys.path.insert(0, src) if src not in sys.path else None' [MESSAGES CONTROL] @@ -17,5 +17,5 @@ disable= max-line-length=100 [TYPECHECK] -ignored-modules=xarray,xesmf,cmor,numpy,geocat.comp,yaml,cftime,pytest,dask,netCDF4 +ignored-modules=xarray,xesmf,cmor,numpy,geocat.comp,yaml,cftime,pytest,dask,gents ignored-classes=cmor diff --git a/tests/test_timeseries_synthetic.py b/tests/test_timeseries_synthetic.py new file mode 100644 index 0000000..838a05c --- /dev/null +++ b/tests/test_timeseries_synthetic.py @@ -0,0 +1,218 @@ +"""Synthetic end-to-end checks of time series generation, by realm and frequency. + +These tests build throwaway history files for every realm and frequency the +include-pattern tables declare, then run the real GenTS collection and time +series steps over them. Nothing here touches model output on disk, so the whole +sweep runs on a laptop in seconds. + +What this is meant to catch: + +* an include pattern that selects nothing, which the driver reports only as + "No input files to process" followed by a clean exit; +* a realm or frequency quietly dropped from a table; +* a GenTS API change, of the kind that deprecated ``include_patterns`` in 1.4.0; +* a stream naming convention that moved on without the tables following. + +What it deliberately does not cover: regridding (the ESMF weights live on NIRD), +CMOR table and controlled-vocabulary compliance, and anything about whether the +numbers are physically meaningful. Those belong to the on-cluster tier. +""" + +import glob +import os +import warnings + +import numpy as np +import pytest +import xarray as xr +from gents.hfcollection import HFCollection +from gents.timeseries import TSCollection + +from cmip7_prep.include_patterns import all_include_patterns, load_include_patterns + +MODELS = ["noresm", "cesm"] + +# Days between samples, per frequency, so a stream's time axis is at least +# plausible for the frequency it claims to be. +STEP_DAYS = { + "1hr": 1.0 / 24.0, + "3hr": 3.0 / 24.0, + "6hr": 6.0 / 24.0, + "day": 1.0, + "mon": 30.0, + "yr": 365.0, +} + +# Neutral variable names: these tests are about file selection and series +# assembly, not about any realm's actual variables. +VARIABLES = ["FIELD_A", "FIELD_B"] + +RECORDS_PER_FILE = 2 +FILES_PER_STREAM = 2 + + +def _cases(): + """Return every (model, realm, frequency) the tables declare.""" + cases = [] + for model in MODELS: + for realm, by_frequency in load_include_patterns(model).items(): + for frequency in sorted(by_frequency): + cases.append((model, realm, frequency)) + return cases + + +# pytest labels each case from the values themselves, e.g. "noresm-seaIce-day", +# so no explicit ids are needed. +CASES = _cases() + + +def _patterns(model, realm, frequency): + """Return the include patterns for one case, as the driver asks for them.""" + ice_sheet = "gris" if realm == "landIce" else None + return all_include_patterns(model, realm, ice_sheet, None, [frequency]) + + +def _filename(fragment, index): + """Return a history file name embedding an include-pattern fragment. + + Patterns are stream fragments such as ``cam.h0a``, ``cice.h1.`` or + ``cism.gris.h``; stripping the dots lets one rule serve all of them. + """ + return f"case.{fragment.strip('.')}.{1526 + index:04d}-01-01-00000.nc" + + +def _write_history_file(path, frequency, first_record): + """Write a minimal history file with RECORDS_PER_FILE samples. + + Time is left as raw offsets with CF attributes set by hand, so the file + looks like model output rather than something xarray round-tripped: no + fill value on the coordinates, and an unlimited time dimension. + """ + step = STEP_DAYS[frequency] + offsets = np.array( + [(first_record + n) * step for n in range(RECORDS_PER_FILE)], dtype="f8" + ) + dataset = xr.Dataset( + { + name: ( + ("time", "lat", "lon"), + np.zeros((RECORDS_PER_FILE, 2, 2), dtype="f4"), + {"units": "1"}, + ) + for name in VARIABLES + }, + coords={ + "time": ("time", offsets), + "lat": ("lat", np.array([-45.0, 45.0])), + "lon": ("lon", np.array([0.0, 180.0])), + }, + ) + dataset["time_bnds"] = ( + ("time", "nbnd"), + np.array([[o, o + step] for o in offsets]), + ) + dataset["time"].attrs.update( + units="days since 1850-01-01", calendar="noleap", bounds="time_bnds" + ) + no_fill = {name: {"_FillValue": None} for name in ("time", "lat", "lon")} + dataset.to_netcdf(path, unlimited_dims=["time"], encoding=no_fill) + + +def _build_stream(directory, fragment, frequency): + """Write FILES_PER_STREAM consecutive history files for one stream.""" + written = [] + for index in range(FILES_PER_STREAM): + path = os.path.join(directory, _filename(fragment, index)) + _write_history_file(path, frequency, index * RECORDS_PER_FILE) + written.append(path) + return written + + +# ------------------------------------------------------------------ patterns + + +@pytest.mark.parametrize("model,realm,frequency", CASES) +def test_patterns_are_declared(model, realm, frequency): + """Every declared realm and frequency yields at least one pattern.""" + patterns = _patterns(model, realm, frequency) + assert patterns, f"{model}/{realm}/{frequency} declares no patterns" + assert all(isinstance(p, str) and p for p in patterns) + assert not any("{" in p for p in patterns), "unsubstituted placeholder" + + +@pytest.mark.parametrize("model,realm,frequency", CASES) +def test_patterns_select_their_own_files(tmp_path, model, realm, frequency): + """Files named after a pattern are found by the glob the driver builds. + + The driver wraps each pattern as ``**`` and globs the input + directory (see gen_timeseries.py), so that exact form is what is tested. + """ + patterns = _patterns(model, realm, frequency) + for fragment in patterns: + _build_stream(tmp_path, fragment, frequency) + + for fragment in patterns: + found = glob.glob(os.path.join(tmp_path, f"*{fragment}*")) + assert found, f"pattern {fragment!r} matched nothing it should have" + + +# ------------------------------------------------------------- series output + + +@pytest.mark.parametrize("model,realm,frequency", CASES) +def test_timeseries_generated_per_variable(tmp_path, model, realm, frequency): + """GenTS turns each synthetic stream into one time series per variable.""" + inputdir = tmp_path / "hist" + outputdir = tmp_path / "timeseries" + inputdir.mkdir() + outputdir.mkdir() + + patterns = _patterns(model, realm, frequency) + fragment = patterns[0] + written = _build_stream(inputdir, fragment, frequency) + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + collection = HFCollection(str(inputdir), num_processes=1) + collection = collection.include([f"*{fragment}*"]) + assert len(collection) == len(written) + collection.pull_metadata(show_progress=False) + + series = TSCollection(collection, str(outputdir), num_processes=1) + series = series.apply_overwrite("*") + series.execute() + + produced = [ + os.path.basename(p) + for p in glob.glob(str(outputdir / "**" / "*.nc"), recursive=True) + ] + assert produced, f"{model}/{realm}/{frequency}: no time series written" + for variable in VARIABLES: + matching = [name for name in produced if variable in name] + assert len(matching) == 1, ( + f"{model}/{realm}/{frequency}: expected one series for {variable}, " + f"got {matching}" + ) + + +# -------------------------------------------------------------- API contract + + +def test_unknown_realm_is_rejected(): + """A realm no table defines raises rather than returning nothing.""" + with pytest.raises(ValueError, match="No include_patterns"): + all_include_patterns("noresm", "no_such_realm") + + +def test_unknown_frequency_is_rejected(): + """A frequency a realm does not define raises and lists what is available.""" + with pytest.raises(ValueError, match="No include_patterns"): + all_include_patterns("noresm", "seaIce", None, None, ["1hr"]) + + +def test_ice_sheet_placeholder_is_substituted(): + """The landIce patterns carry the ice sheet given, not a placeholder.""" + for sheet in ("gris", "ais"): + patterns = all_include_patterns("noresm", "landIce", sheet, None, ["yr"]) + assert patterns + assert all(sheet in p for p in patterns), patterns From ad3bc4789a388b8f02a0a60baa1dc9948f894995 Mon Sep 17 00:00:00 2001 From: mvertens Date: Sun, 4 Oct 2026 14:29:18 +0200 Subject: [PATCH 04/16] addition of a one-command reference-case runner with CV (controlled vocabulary) validation --- scripts/run_reference_case.py | 315 +++++++++++++++++++++++++++ src/cmip7_prep/cv_lookup.py | 89 ++++++++ src/cmip7_prep/reference_run.py | 365 ++++++++++++++++++++++++++++++++ tests/test_cv_lookup.py | 125 +++++++++++ tests/test_reference_run.py | 320 ++++++++++++++++++++++++++++ 5 files changed, 1214 insertions(+) create mode 100644 scripts/run_reference_case.py create mode 100644 src/cmip7_prep/cv_lookup.py create mode 100644 src/cmip7_prep/reference_run.py create mode 100644 tests/test_cv_lookup.py create mode 100644 tests/test_reference_run.py diff --git a/scripts/run_reference_case.py b/scripts/run_reference_case.py new file mode 100644 index 0000000..f063ab1 --- /dev/null +++ b/scripts/run_reference_case.py @@ -0,0 +1,315 @@ +#!/usr/bin/env python3 + +"""Run the whole chain over one archived case, for every realm, in one command. + + python scripts/run_reference_case.py --case-dir --outdir + +For each realm the model declares, this runs gen_timeseries.py, then +cmor_driver.py and validate_cmor_output.py for each of that realm's +frequencies, and prints one table at the end saying what succeeded, what +failed, and how much was produced. + +It is meant to be run on the machine holding the data, interactively -- there +is no batch submission here. A full run takes a while, so start it under +``tmux`` or ``screen`` if the connection might drop. Use ``--dry-run`` first to +see exactly which commands would be issued. + +Validation is what decides whether a run was good: validate_cmor_output.py is +invoked with --strict, so a missing variable or a CMOR log error fails the run. +""" + +import argparse +import logging +import subprocess +import sys +import time +from pathlib import Path + +_LOCAL_PATH = Path(__file__).resolve().parent +sys.path.insert(0, str(_LOCAL_PATH.parent / "src")) + +# pylint: disable=wrong-import-position +from cmip7_prep.cv_lookup import validate as validate_cv +from cmip7_prep.reference_run import ( + ATM_RESOLUTIONS, + STAGES, + Plan, + Step, + build_plan, +) + +logger = logging.getLogger("run_reference_case") + + +def parse_arguments(): + """Parse and return the command-line arguments.""" + parser = argparse.ArgumentParser( + description=( + "Run time series generation, CMORization and validation over one " + "archived case, for every realm, and report the result." + ) + ) + + # Grouped so --help states plainly which arguments must be given. All five + # describe the case itself: a wrong value produces output that looks valid, + # so none of them may be defaulted. + required = parser.add_argument_group("required arguments") + required.add_argument( + "--case-dir", + required=True, + help="Archived case root, holding atm/, lnd/, ice/, glc/ and ocn/", + ) + required.add_argument( + "--outdir", + required=True, + help="Directory to write timeseries/ and cmor/ output into", + ) + required.add_argument( + "--model", + choices=["noresm", "cesm"], + required=True, + help=( + "Model whose include patterns and variable mappings to use. " + "The wrong one produces plausible but wrong output." + ), + ) + required.add_argument( + "--atmos-resolution", + choices=list(ATM_RESOLUTIONS), + required=True, + help=( + "Grid the atmosphere and land were run on. Ocean and sea ice are " + "always on the model's own tripolar grid, and land ice is written " + "on its native projected grid, so those are derived from this." + ), + ) + required.add_argument( + "--experiment", + required=True, + help=( + "CMIP7 experiment_id of the case, e.g. piControl. Checked against " + "the controlled vocabulary before any work starts." + ), + ) + + selection = parser.add_argument_group("selecting what to run") + selection.add_argument( + "--realms", + nargs="+", + default=None, + metavar="REALM", + help="Realms to process (default: every realm the model declares)", + ) + selection.add_argument( + "--frequencies", + nargs="+", + default=None, + metavar="FREQ", + help="Frequencies to process (default: every frequency each realm declares)", + ) + selection.add_argument( + "--years", + default=None, + help=( + "Years to process, as gen_timeseries.py's --years-spec: " + "first:last:increment, e.g. 1526:1527:2 for two years in one file" + ), + ) + selection.add_argument( + "--stages", + nargs="+", + choices=list(STAGES), + default=list(STAGES), + help="Stages to run (default: all three)", + ) + selection.add_argument( + "--ice-sheet", + choices=["gris", "ais"], + default=None, + help="Ice sheet for the landIce realm (default: gris)", + ) + + behaviour = parser.add_argument_group("how to run") + behaviour.add_argument( + "--workers", + type=int, + default=4, + help=( + "Worker processes per stage (default: 4, chosen to be polite on a " + "shared interactive node rather than to be fast)" + ), + ) + behaviour.add_argument( + "--dry-run", + action="store_true", + help="Print the commands that would run, and exit", + ) + behaviour.add_argument( + "--keep-going", + action="store_true", + help="Continue after a failing step instead of stopping", + ) + behaviour.add_argument( + "--tables-root", + default=None, + help=( + "cmip7-cmor-tables checkout, used to check --experiment against the " + "controlled vocabulary (default: the one beside this repo)" + ), + ) + behaviour.add_argument("--log-level", default="INFO", help="Default: INFO") + return parser.parse_args() + + +def run_step(step: Step, log_dir: Path) -> dict: + """Run one step, tee its output to a log file, and return the outcome.""" + log_path = log_dir / step.log_name + logger.info("[%s] %s", step.key, " ".join(step.command)) + started = time.time() + with open(log_path, "w", encoding="utf-8") as log_file: + log_file.write(" ".join(step.command) + "\n\n") + log_file.flush() + completed = subprocess.run( + step.command, + stdout=log_file, + stderr=subprocess.STDOUT, + check=False, + ) + elapsed = time.time() - started + outcome = { + "key": step.key, + "stage": step.stage, + "realm": step.realm, + "frequency": step.frequency, + "returncode": completed.returncode, + "seconds": elapsed, + "log": str(log_path), + } + if completed.returncode == 0: + logger.info("[%s] ok in %.1fs", step.key, elapsed) + else: + logger.error( + "[%s] FAILED (exit %d) after %.1fs -- see %s", + step.key, + completed.returncode, + elapsed, + log_path, + ) + return outcome + + +def print_plan(plan: Plan) -> None: + """Print the commands a run would issue, without running them.""" + for stage in STAGES: + steps = plan.for_stage(stage) + if not steps: + continue + print(f"\n=== {stage} ({len(steps)} step(s)) ===") + for step in steps: + print(f" {step.key}") + print(f" {' '.join(step.command)}") + if plan.skipped: + print("\n=== skipped ===") + for reason in plan.skipped: + print(f" {reason}") + + +def print_report(outcomes: list[dict], plan: Plan) -> None: + """Print the end-of-run summary table.""" + print("\n" + "=" * 78) + print("REFERENCE CASE RUN") + print("=" * 78) + + width = max((len(o["key"]) for o in outcomes), default=10) + for outcome in outcomes: + status = ( + "ok" if outcome["returncode"] == 0 else f"FAILED ({outcome['returncode']})" + ) + print(f" {outcome['key']:<{width}} {outcome['seconds']:7.1f}s {status}") + + failed = [o for o in outcomes if o["returncode"] != 0] + print("-" * 78) + print(f" {len(outcomes) - len(failed)} ok, {len(failed)} failed") + + if plan.skipped: + print("\n skipped:") + for reason in plan.skipped: + print(f" {reason}") + + if failed: + print("\n logs for the failures:") + for outcome in failed: + print(f" {outcome['key']}: {outcome['log']}") + + +def main(): + """Build the plan, run it, and report.""" + args = parse_arguments() + logging.basicConfig( + level=getattr(logging, args.log_level.upper(), logging.INFO), + format="%(asctime)s %(levelname)s %(name)s: %(message)s", + ) + + case_dir = Path(args.case_dir) + if not case_dir.is_dir(): + logger.error("No such case directory: %s", case_dir) + sys.exit(2) + + # Fail on a bad experiment now rather than after the time series stage: CMOR + # would reject it at write time, by which point the slow part has run. + tables_root = Path(args.tables_root or _LOCAL_PATH.parent / "cmip7-cmor-tables") + try: + validate_cv(tables_root, "experiment_id", args.experiment) + except (ValueError, KeyError, FileNotFoundError) as exc: + logger.error("%s", exc) + sys.exit(2) + + plan = build_plan( + case_dir=case_dir, + outdir=args.outdir, + model=args.model, + realms=args.realms, + frequencies=args.frequencies, + years=args.years, + stages=args.stages, + atmos_resolution=args.atmos_resolution, + experiment=args.experiment, + workers=args.workers, + ice_sheet=args.ice_sheet, + scripts_dir=_LOCAL_PATH, + ) + + if not plan.steps: + logger.error("Nothing to do. Skipped: %s", plan.skipped or "(nothing)") + sys.exit(2) + + logger.info( + "%d step(s) over realms %s", + len(plan.steps), + sorted({step.realm for step in plan.steps}), + ) + + if args.dry_run: + print_plan(plan) + return + + outdir = Path(args.outdir) + log_dir = outdir / "logs" + log_dir.mkdir(parents=True, exist_ok=True) + + outcomes = [] + for step in plan.steps: + outcome = run_step(step, log_dir) + outcomes.append(outcome) + if outcome["returncode"] != 0 and not args.keep_going: + logger.error("Stopping at the first failure; pass --keep-going to continue") + break + + print_report(outcomes, plan) + + failed = any(outcome["returncode"] != 0 for outcome in outcomes) + sys.exit(1 if failed else 0) + + +if __name__ == "__main__": + main() diff --git a/src/cmip7_prep/cv_lookup.py b/src/cmip7_prep/cv_lookup.py new file mode 100644 index 0000000..0416f44 --- /dev/null +++ b/src/cmip7_prep/cv_lookup.py @@ -0,0 +1,89 @@ +"""Reading the controlled vocabulary, to reject bad values before a run starts. + +CMOR validates global attributes against the controlled vocabulary when it +writes, and refuses the write if one is not registered. That check is the right +one, but it arrives late: by then the time series stage has already run, which +is the expensive part. The helpers here read the same file CMOR reads, so a +misspelled experiment can be rejected in the first second instead of the last. + +The vocabulary lives in the CMOR tables checkout, not in this package, so every +function takes the tables root. It is read on demand and cached, since a run +asks about a handful of values at most. +""" + +from __future__ import annotations + +import difflib +import json +import os +from functools import lru_cache +from pathlib import Path + +# Where the vocabulary sits relative to the tables checkout. +CV_RELATIVE_PATH = Path("tables-cvs") / "cmor-cvs.json" + +# The one key the file wraps everything in. +CV_ROOT_KEY = "CV" + + +def cv_path(tables_root: os.PathLike | str) -> Path: + """Return the path of the controlled vocabulary within a tables checkout.""" + return Path(tables_root) / CV_RELATIVE_PATH + + +@lru_cache(maxsize=4) +def _load(tables_root: str) -> dict: + """Return the vocabulary's contents, cached per tables checkout.""" + path = cv_path(tables_root) + if not path.is_file(): + raise FileNotFoundError( + f"No controlled vocabulary at {path}. Point --tables-root at a " + "cmip7-cmor-tables checkout, or update the one you have." + ) + document = json.loads(path.read_text(encoding="utf-8")) + return document.get(CV_ROOT_KEY, document) + + +def allowed_values(tables_root: os.PathLike | str, key: str) -> list[str]: + """Return the registered values of one controlled attribute, sorted. + + ``key`` is a vocabulary key such as ``experiment_id`` or ``source_id``. + Some attributes nest -- the licences live at ``license.license_id`` -- so a + dotted path walks down to them. An unknown key raises, listing what is + controlled at that level, so a typo in the key is as loud as one in a value. + """ + entry = _load(str(tables_root)) + walked = [] + for part in key.split("."): + if not isinstance(entry, dict) or part not in entry: + where = ".".join(walked) or "the vocabulary" + raise KeyError( + f"{key!r} is not a controlled attribute in " + f"{cv_path(tables_root)}; {where} controls: " + f"{sorted(entry) if isinstance(entry, dict) else entry}" + ) + entry = entry[part] + walked.append(part) + if isinstance(entry, dict): + return sorted(entry) + if isinstance(entry, list): + return sorted(str(value) for value in entry) + return [str(entry)] + + +def validate(tables_root: os.PathLike | str, key: str, value: str) -> None: + """Raise ValueError unless ``value`` is registered for ``key``. + + The message names the nearest registered values, since the usual cause is a + small misspelling and 75 experiment ids are too many to print. + """ + allowed = allowed_values(tables_root, key) + if value in allowed: + return + suggestions = difflib.get_close_matches(value, allowed, n=5, cutoff=0.5) + hint = ( + f" Did you mean: {', '.join(suggestions)}?" + if suggestions + else f" It controls {len(allowed)} value(s); see {cv_path(tables_root)}." + ) + raise ValueError(f"{key} {value!r} is not in the controlled vocabulary.{hint}") diff --git a/src/cmip7_prep/reference_run.py b/src/cmip7_prep/reference_run.py new file mode 100644 index 0000000..cf325ce --- /dev/null +++ b/src/cmip7_prep/reference_run.py @@ -0,0 +1,365 @@ +"""Plan a full-chain run over one archived case, for every realm. + +This is the planning half of ``scripts/run_reference_case.py``: it works out +which commands a complete run needs, in which order, without running anything. +Keeping it here rather than in the script means the plan can be tested, and +printed for inspection before a long run is started. + +The archive layout assumed is CIME's short-term archive: + + //hist/ history files, one component per realm + //proc/... where processed output conventionally goes + +though the output root is given separately, so a run never has to write back +into an archive it may not own. +""" + +from __future__ import annotations + +import os +from dataclasses import dataclass, field +from pathlib import Path +from typing import Sequence + +from .include_patterns import load_include_patterns + +# Which component directory holds each realm's history files. Mirrors +# REALM_COMPONENT_MAP in scripts/cmor_driver.py. +REALM_COMPONENT = { + "atmos": "atm", + "aerosol": "atm", + "atmosChem": "atm", + "land": "lnd", + "ocean": "ocn", + "ocnBgchem": "ocn", + "seaIce": "ice", + "landIce": "glc", +} + +# The ice sheet to process when a realm is per-ice-sheet and none was named. +DEFAULT_ICE_SHEET = "gris" + +# Grids the atmosphere and land may be run on. 'regular' means the output +# already carries lat/lon and is not regridded. +ATM_RESOLUTIONS = ("ne30", "ne16", "regular") + +# Realms whose input grid is the atmosphere/land grid the case was run on. +ATM_GRID_REALMS = frozenset({"atmos", "atmosChem", "aerosol", "land"}) + +# Realms written on their native grid, with no ESMF regridding. CISM land-ice +# output carries projected x/y coordinates which cmor_writer georeferences with +# the ice sheet's own map projection (see cism_grid.project_xy_to_latlon), so no +# weight files are involved and the resolution would be ignored. The value +# below is a pass-through: it names the entry that builds a regridder and +# discards it, which is what the unregridded path expects. +NATIVE_GRID_REALMS = frozenset({"landIce"}) +NATIVE_GRID = "regular" + +# The ocean/sea-ice grid each model writes on. These realms never share the +# atmosphere's grid, so passing one an atmosphere resolution would regrid +# through the wrong weights and produce plausible but wrong output. +OCEAN_GRID = {"noresm": "tnx1v4", "cesm": "tx2_3v2"} + +OCEAN_GRID_REALMS = frozenset({"ocean", "ocnBgchem", "seaIce"}) + +STAGES = ("timeseries", "cmor", "validate") + + +@dataclass +class Step: + """One command the run will execute. + + ``key`` identifies the step in the final report; ``realm`` and ``frequency`` + are carried so a failure can be attributed without parsing the command. + """ + + key: str + stage: str + realm: str + command: list[str] + frequency: str | None = None + log_name: str = "" + + def __post_init__(self) -> None: + if not self.log_name: + self.log_name = f"{self.key}.log" + + +@dataclass +class Plan: + """A complete run: its steps, and the directories they use.""" + + steps: list[Step] = field(default_factory=list) + timeseries_root: Path | None = None + cmor_root: Path | None = None + skipped: list[str] = field(default_factory=list) + + def for_stage(self, stage: str) -> list[Step]: + """Return the steps belonging to one stage.""" + return [step for step in self.steps if step.stage == stage] + + +def realms_for(model: str) -> list[str]: + """Return the realms a model's include-pattern table declares.""" + return list(load_include_patterns(model)) + + +def frequencies_for(model: str, realm: str) -> list[str]: + """Return the frequencies declared for one model and realm.""" + return sorted(load_include_patterns(model)[realm]) + + +def resolution_for(model: str, realm: str, atm_resolution: str) -> str: + """Return the input grid name for one realm. + + Every realm's grid is derived, so no caller has to know the model's grid + layout. ``atm_resolution`` is the grid the atmosphere and land were run on + and applies only to those realms; ocean and sea ice are always on the + model's tripolar grid; land ice is written on its native projected grid and + is not regridded at all. + """ + if atm_resolution not in ATM_RESOLUTIONS: + raise ValueError( + f"Unknown atmosphere resolution {atm_resolution!r}; " + f"choose from {list(ATM_RESOLUTIONS)}" + ) + if realm in ATM_GRID_REALMS: + return atm_resolution + if realm in OCEAN_GRID_REALMS: + try: + return OCEAN_GRID[model] + except KeyError: + raise ValueError( + f"No ocean grid known for model={model!r}; " + f"known models: {sorted(OCEAN_GRID)}" + ) from None + if realm in NATIVE_GRID_REALMS: + return NATIVE_GRID + raise ValueError(f"No input grid known for realm {realm!r}") + + +def history_dir(case_dir: os.PathLike | str, realm: str) -> Path: + """Return the history directory holding one realm's output.""" + component = REALM_COMPONENT[realm] + return Path(case_dir) / component / "hist" + + +def build_plan( + case_dir: os.PathLike | str, + outdir: os.PathLike | str, + *, + model: str, + realms: Sequence[str] | None = None, + frequencies: Sequence[str] | None = None, + years: str | None = None, + stages: Sequence[str] = STAGES, + atmos_resolution: str, + experiment: str, + workers: int = 4, + ice_sheet: str | None = None, + scripts_dir: os.PathLike | str = "scripts", +) -> Plan: + """Return the full plan for one reference-case run. + + ``realms`` of None takes every realm the model declares; a realm whose + history directory is absent is recorded in ``Plan.skipped`` rather than + failing the run, since an archived case need not hold every component. + + ``atmos_resolution`` is the grid the atmosphere and land were run on; every + other realm's grid is derived from it or from the model. + + ``years`` is passed through to gen_timeseries.py as ``--years-spec`` and so + uses its format, ``first:last:increment``. + """ + unknown = [stage for stage in stages if stage not in STAGES] + if unknown: + raise ValueError(f"Unknown stage(s) {unknown}; choose from {list(STAGES)}") + + case_dir = Path(case_dir) + outdir = Path(outdir) + scripts_dir = Path(scripts_dir) + timeseries_root = outdir / "timeseries" + cmor_root = outdir / "cmor" + + plan = Plan(timeseries_root=timeseries_root, cmor_root=cmor_root) + wanted_realms = list(realms) if realms else realms_for(model) + + for realm in wanted_realms: + if realm not in REALM_COMPONENT: + raise ValueError( + f"Realm {realm!r} has no component directory; " + f"known realms: {sorted(REALM_COMPONENT)}" + ) + inputdir = history_dir(case_dir, realm) + if not inputdir.is_dir(): + plan.skipped.append(f"{realm}: no history directory at {inputdir}") + continue + + realm_frequencies = _realm_frequencies(model, realm, frequencies) + if not realm_frequencies: + plan.skipped.append( + f"{realm}: none of the requested frequencies are declared" + ) + continue + + sheet = ice_sheet or DEFAULT_ICE_SHEET if realm == "landIce" else None + ts_dir = timeseries_root / realm + + realm_resolution = resolution_for(model, realm, atmos_resolution) + + if "timeseries" in stages: + plan.steps.append( + _timeseries_step( + scripts_dir=scripts_dir, + realm=realm, + inputdir=inputdir, + ts_dir=ts_dir, + model=model, + frequencies=realm_frequencies, + years=years, + workers=workers, + sheet=sheet, + ) + ) + + for frequency in realm_frequencies: + if "cmor" in stages: + plan.steps.append( + _cmor_step( + scripts_dir=scripts_dir, + realm=realm, + frequency=frequency, + ts_dir=ts_dir, + cmor_root=cmor_root, + model=model, + resolution=realm_resolution, + experiment=experiment, + workers=workers, + sheet=sheet, + ) + ) + if "validate" in stages: + plan.steps.append( + _validate_step( + scripts_dir=scripts_dir, + realm=realm, + frequency=frequency, + cmor_root=cmor_root, + model=model, + experiment=experiment, + ) + ) + + return plan + + +def _realm_frequencies( + model: str, realm: str, requested: Sequence[str] | None +) -> list[str]: + """Return the declared frequencies of a realm, narrowed to those requested.""" + declared = frequencies_for(model, realm) + if not requested: + return declared + return [frequency for frequency in declared if frequency in requested] + + +def _timeseries_step( + *, scripts_dir, realm, inputdir, ts_dir, model, frequencies, years, workers, sheet +) -> Step: + """Return the gen_timeseries.py step for one realm.""" + command = [ + "python", + str(Path(scripts_dir) / "gen_timeseries.py"), + "--inputdir", + str(inputdir), + "--outputdir", + str(ts_dir), + "--realm", + realm, + "--model", + model, + "--workers", + str(workers), + "--frequency", + *frequencies, + ] + if years: + command += ["--years-spec", years] + if sheet: + command += ["--ice-sheet", sheet] + return Step( + key=f"timeseries-{realm}", stage="timeseries", realm=realm, command=command + ) + + +def _cmor_step( + *, + scripts_dir, + realm, + frequency, + ts_dir, + cmor_root, + model, + resolution, + experiment, + workers, + sheet, +) -> Step: + """Return the cmor_driver.py step for one realm and frequency.""" + command = [ + "python", + str(Path(scripts_dir) / "cmor_driver.py"), + "--realm", + realm, + "--frequency", + frequency, + "--tsdir", + str(ts_dir), + "--outdir", + str(cmor_root), + "--model", + model, + "--resolution", + resolution, + "--experiment", + experiment, + "--workers", + str(workers), + ] + if sheet: + command += ["--ice-sheet", sheet] + return Step( + key=f"cmor-{realm}-{frequency}", + stage="cmor", + realm=realm, + frequency=frequency, + command=command, + ) + + +def _validate_step( + *, scripts_dir, realm, frequency, cmor_root, model, experiment +) -> Step: + """Return the validate_cmor_output.py step for one realm and frequency.""" + command = [ + "python", + str(Path(scripts_dir) / "validate_cmor_output.py"), + "--model", + model, + "--realm", + realm, + "--frequency", + frequency, + "--experiment", + experiment, + "--root-output-path", + str(cmor_root), + "--strict", + ] + return Step( + key=f"validate-{realm}-{frequency}", + stage="validate", + realm=realm, + frequency=frequency, + command=command, + ) diff --git a/tests/test_cv_lookup.py b/tests/test_cv_lookup.py new file mode 100644 index 0000000..f4fab3b --- /dev/null +++ b/tests/test_cv_lookup.py @@ -0,0 +1,125 @@ +"""Tests for reading the controlled vocabulary to reject bad values early.""" + +import json + +import pytest + +from cmip7_prep import cv_lookup + +TABLES_ROOT = "cmip7-cmor-tables" + + +@pytest.fixture(name="fake_tables") +def fake_tables_fixture(tmp_path): + """Return a tables root holding a small stand-in vocabulary.""" + cv_dir = tmp_path / "tables-cvs" + cv_dir.mkdir() + (cv_dir / "cmor-cvs.json").write_text( + json.dumps( + { + "CV": { + "experiment_id": { + "piControl": {"description": "control"}, + "historical": {"description": "historical"}, + "abrupt-4xCO2": {"description": "idealised"}, + }, + "license": ["CC-BY-4.0"], + "required_global_attributes": ["source_id"], + } + } + ), + encoding="utf-8", + ) + cv_lookup._load.cache_clear() # pylint: disable=protected-access + return tmp_path + + +# ---------------------------------------------------------------- the real CV + + +class TestAgainstTheRepoVocabulary: + """Checks against the vocabulary actually checked out in this repo.""" + + def test_experiment_ids_are_readable(self): + """The repo's vocabulary yields a non-trivial list of experiments.""" + experiments = cv_lookup.allowed_values(TABLES_ROOT, "experiment_id") + assert len(experiments) > 20 + assert "piControl" in experiments + + def test_a_valid_experiment_passes(self): + """A registered experiment validates silently.""" + cv_lookup.validate(TABLES_ROOT, "experiment_id", "piControl") + + def test_the_license_we_fixed_is_registered(self): + """CC-BY-4.0 is the registered spelling, not CC-BY-4-0. + + A typo here is what stopped CMOR writing anything at all, so it is + worth asserting the vocabulary still agrees with data/cmor_dataset.json. + """ + licenses = cv_lookup.allowed_values(TABLES_ROOT, "license.license_id") + assert "CC-BY-4.0" in licenses + assert "CC-BY-4-0" not in licenses + + +# --------------------------------------------------------- the stand-in tables + + +class TestLookup: + """Tests for allowed_values and validate.""" + + def test_values_are_sorted(self, fake_tables): + """Registered values come back sorted, for stable messages.""" + assert cv_lookup.allowed_values(fake_tables, "experiment_id") == [ + "abrupt-4xCO2", + "historical", + "piControl", + ] + + def test_a_list_valued_key_works(self, fake_tables): + """A key whose value is a list is read the same way as a mapping.""" + assert cv_lookup.allowed_values(fake_tables, "license") == ["CC-BY-4.0"] + + def test_a_dotted_path_walks_into_nested_keys(self): + """The licences nest under license.license_id in the real vocabulary.""" + assert "CC-BY-4.0" in cv_lookup.allowed_values( + TABLES_ROOT, "license.license_id" + ) + + def test_a_dotted_path_that_goes_nowhere_raises(self): + """A wrong path says what is controlled at the level it reached.""" + with pytest.raises(KeyError, match="not a controlled attribute"): + cv_lookup.allowed_values(TABLES_ROOT, "license.no_such_level") + + def test_unknown_key_lists_what_is_controlled(self, fake_tables): + """A typo in the attribute name is as loud as one in the value.""" + with pytest.raises(KeyError, match="not a controlled attribute"): + cv_lookup.allowed_values(fake_tables, "experimnet_id") + + def test_missing_vocabulary_says_where_it_looked(self, tmp_path): + """A tables root without a vocabulary names the path it wanted.""" + cv_lookup._load.cache_clear() # pylint: disable=protected-access + with pytest.raises(FileNotFoundError, match="cmor-cvs.json"): + cv_lookup.allowed_values(tmp_path, "experiment_id") + + +class TestValidate: + """Tests for the rejection messages.""" + + def test_valid_value_is_accepted(self, fake_tables): + """A registered value raises nothing.""" + cv_lookup.validate(fake_tables, "experiment_id", "historical") + + def test_near_miss_suggests_the_real_thing(self, fake_tables): + """A misspelling is answered with the nearest registered values.""" + with pytest.raises(ValueError, match="Did you mean: piControl"): + cv_lookup.validate(fake_tables, "experiment_id", "piControll") + + def test_wild_value_points_at_the_file(self, fake_tables): + """With nothing close, the message says where the list lives.""" + with pytest.raises(ValueError, match="cmor-cvs.json"): + cv_lookup.validate(fake_tables, "experiment_id", "banana") + + def test_case_matters(self, fake_tables): + """Vocabulary values are case sensitive, as CMOR treats them.""" + with pytest.raises(ValueError, match="not in the controlled vocabulary"): + cv_lookup.validate(fake_tables, "experiment_id", "picontrol") diff --git a/tests/test_reference_run.py b/tests/test_reference_run.py new file mode 100644 index 0000000..424467e --- /dev/null +++ b/tests/test_reference_run.py @@ -0,0 +1,320 @@ +"""Tests for planning a full-chain reference-case run.""" + +import pytest + +from cmip7_prep.reference_run import ( + ATM_RESOLUTIONS, + NATIVE_GRID, + OCEAN_GRID, + REALM_COMPONENT, + STAGES, + build_plan, + frequencies_for, + history_dir, + realms_for, + resolution_for, +) + + +@pytest.fixture(name="case_dir") +def case_dir_fixture(tmp_path): + """Return a case directory with every component's history directory.""" + for component in sorted(set(REALM_COMPONENT.values())): + (tmp_path / component / "hist").mkdir(parents=True) + return tmp_path + + +def _plan(case_dir, outdir, **kwargs): + """Call build_plan with a model, which it requires. + + The model has no default on purpose: the wrong one yields plausible but + wrong output, so a caller must say which it means. + """ + kwargs.setdefault("model", "noresm") + kwargs.setdefault("atmos_resolution", "ne16") + kwargs.setdefault("experiment", "piControl") + return build_plan(case_dir, outdir, **kwargs) + + +def _command_of(plan, key): + """Return the command of the step with one key.""" + for step in plan.steps: + if step.key == key: + return step.command + raise AssertionError(f"no step {key!r} in {[s.key for s in plan.steps]}") + + +# ------------------------------------------------------------------- lookups + + +def test_realms_come_from_the_tables(): + """The realms planned for are the ones the include-pattern table declares.""" + assert realms_for("noresm") == [ + "atmos", + "atmosChem", + "aerosol", + "land", + "seaIce", + "landIce", + ] + + +def test_every_realm_has_a_component_directory(): + """No declared realm is missing from the component map. + + A realm present in the tables but absent here would be silently skipped, + so the two must stay in step. + """ + for model in ("noresm", "cesm"): + for realm in realms_for(model): + assert realm in REALM_COMPONENT, f"{model}/{realm} has no component" + + +def test_history_dir_follows_the_archive_layout(tmp_path): + """History lives under //hist.""" + assert history_dir(tmp_path, "seaIce") == tmp_path / "ice" / "hist" + assert history_dir(tmp_path, "atmosChem") == tmp_path / "atm" / "hist" + + +# ---------------------------------------------------------------- plan shape + + +class TestPlanShape: + """Tests for which steps a plan contains.""" + + def test_covers_every_realm_by_default(self, case_dir, tmp_path): + """One timeseries step per realm, when no realms are named.""" + plan = _plan(case_dir, tmp_path / "out") + timeseries = plan.for_stage("timeseries") + assert {step.realm for step in timeseries} == set(realms_for("noresm")) + assert len(timeseries) == 6 + + def test_cmor_and_validate_are_per_frequency(self, case_dir, tmp_path): + """CMOR and validation run once per realm and frequency.""" + plan = _plan(case_dir, tmp_path / "out", realms=["seaIce"]) + assert [step.key for step in plan.for_stage("cmor")] == [ + "cmor-seaIce-day", + "cmor-seaIce-mon", + ] + assert [step.key for step in plan.for_stage("validate")] == [ + "validate-seaIce-day", + "validate-seaIce-mon", + ] + + def test_realms_can_be_narrowed(self, case_dir, tmp_path): + """Naming realms restricts the plan to them.""" + plan = _plan(case_dir, tmp_path / "out", realms=["land", "seaIce"]) + assert {step.realm for step in plan.steps} == {"land", "seaIce"} + + def test_frequencies_can_be_narrowed(self, case_dir, tmp_path): + """Naming frequencies restricts each realm to the ones it declares.""" + plan = _plan( + case_dir, tmp_path / "out", realms=["atmos"], frequencies=["mon", "day"] + ) + assert [step.frequency for step in plan.for_stage("cmor")] == ["day", "mon"] + + def test_frequency_a_realm_does_not_declare_is_skipped(self, case_dir, tmp_path): + """Asking seaIce for 1hr skips it rather than planning a doomed step.""" + plan = _plan(case_dir, tmp_path / "out", realms=["seaIce"], frequencies=["1hr"]) + assert not plan.steps + assert any("seaIce" in reason for reason in plan.skipped) + + def test_stages_can_be_narrowed(self, case_dir, tmp_path): + """Only the requested stages are planned.""" + plan = _plan( + case_dir, tmp_path / "out", realms=["seaIce"], stages=["timeseries"] + ) + assert {step.stage for step in plan.steps} == {"timeseries"} + + def test_unknown_stage_is_rejected(self, case_dir, tmp_path): + """A misspelled stage fails immediately, not mid-run.""" + with pytest.raises(ValueError, match="Unknown stage"): + _plan(case_dir, tmp_path / "out", stages=["timeseries", "cmorise"]) + + def test_missing_component_is_skipped_not_fatal(self, tmp_path): + """A case without every component plans what it can and says so. + + An archived case need not hold every component, so a missing directory + is reported rather than ending the run. + """ + case = tmp_path / "case" + (case / "ice" / "hist").mkdir(parents=True) + plan = _plan(case, tmp_path / "out") + assert {step.realm for step in plan.steps} == {"seaIce"} + assert len(plan.skipped) == 5 + assert all("no history directory" in reason for reason in plan.skipped) + + def test_all_stages_present_by_default(self, case_dir, tmp_path): + """The default plan runs all three stages.""" + plan = _plan(case_dir, tmp_path / "out", realms=["seaIce"]) + assert {step.stage for step in plan.steps} == set(STAGES) + + +# -------------------------------------------------------------- command form + + +class TestCommands: + """Tests for the commands a plan builds.""" + + def test_timeseries_command_points_at_the_right_history_dir( + self, case_dir, tmp_path + ): + """Each realm reads from its own component's history directory.""" + plan = _plan(case_dir, tmp_path / "out", realms=["land"]) + command = _command_of(plan, "timeseries-land") + assert "--inputdir" in command + assert command[command.index("--inputdir") + 1] == str( + case_dir / "lnd" / "hist" + ) + + def test_timeseries_passes_every_declared_frequency(self, case_dir, tmp_path): + """One timeseries run covers all of a realm's frequencies.""" + plan = _plan(case_dir, tmp_path / "out", realms=["land"]) + command = _command_of(plan, "timeseries-land") + at = command.index("--frequency") + assert command[at + 1 : at + 5] == frequencies_for("noresm", "land") + + def test_years_are_passed_through(self, case_dir, tmp_path): + """--years reaches gen_timeseries.py as --years-spec.""" + plan = _plan(case_dir, tmp_path / "out", realms=["seaIce"], years="1526:1527:2") + command = _command_of(plan, "timeseries-seaIce") + assert command[command.index("--years-spec") + 1] == "1526:1527:2" + + def test_no_years_spec_when_not_asked_for(self, case_dir, tmp_path): + """Without --years the flag is absent, so the script's default applies.""" + plan = _plan(case_dir, tmp_path / "out", realms=["seaIce"]) + assert "--years-spec" not in _command_of(plan, "timeseries-seaIce") + + def test_ice_sheet_only_for_landice(self, case_dir, tmp_path): + """landIce gets --ice-sheet; other realms do not.""" + plan = _plan(case_dir, tmp_path / "out", realms=["landIce", "seaIce"]) + assert "--ice-sheet" in _command_of(plan, "timeseries-landIce") + assert "--ice-sheet" not in _command_of(plan, "timeseries-seaIce") + + def test_ice_sheet_can_be_chosen(self, case_dir, tmp_path): + """The ice sheet given is the one passed on.""" + plan = _plan(case_dir, tmp_path / "out", realms=["landIce"], ice_sheet="ais") + command = _command_of(plan, "timeseries-landIce") + assert command[command.index("--ice-sheet") + 1] == "ais" + + def test_cmor_reads_the_timeseries_the_first_stage_writes(self, case_dir, tmp_path): + """The CMOR step's --tsdir is the timeseries step's --outputdir.""" + out = tmp_path / "out" + plan = _plan(case_dir, out, realms=["seaIce"]) + wrote = _command_of(plan, "timeseries-seaIce") + reads = _command_of(plan, "cmor-seaIce-mon") + assert ( + wrote[wrote.index("--outputdir") + 1] == reads[reads.index("--tsdir") + 1] + ) + + def test_validate_reads_where_cmor_writes(self, case_dir, tmp_path): + """Validation's --root-output-path is the CMOR step's --outdir.""" + plan = _plan(case_dir, tmp_path / "out", realms=["seaIce"]) + wrote = _command_of(plan, "cmor-seaIce-mon") + reads = _command_of(plan, "validate-seaIce-mon") + assert ( + wrote[wrote.index("--outdir") + 1] + == reads[reads.index("--root-output-path") + 1] + ) + + def test_validate_is_strict(self, case_dir, tmp_path): + """Validation fails the run on missing variables or log errors.""" + plan = _plan(case_dir, tmp_path / "out", realms=["seaIce"]) + assert "--strict" in _command_of(plan, "validate-seaIce-mon") + + def test_resolution_and_experiment_reach_cmor(self, case_dir, tmp_path): + """Resolution and experiment are passed to the CMOR step. + + Resolution is checked on an atmosphere realm, since that is the only + kind that takes the value given; see TestResolutionPerRealm. + """ + plan = _plan( + case_dir, + tmp_path / "out", + realms=["atmos"], + atmos_resolution="ne30", + experiment="historical", + ) + command = _command_of(plan, "cmor-atmos-mon") + assert command[command.index("--resolution") + 1] == "ne30" + assert command[command.index("--experiment") + 1] == "historical" + + def test_every_step_has_a_distinct_log_name(self, case_dir, tmp_path): + """Logs cannot overwrite one another.""" + plan = _plan(case_dir, tmp_path / "out") + names = [step.log_name for step in plan.steps] + assert len(names) == len(set(names)) + + +@pytest.mark.parametrize("omitted", ["model", "atmos_resolution", "experiment"]) +def test_case_properties_are_required(omitted): + """build_plan refuses to guess a property of the case. + + Model, resolution and experiment each change what the output means while + leaving it looking valid, so none of them may be defaulted. + """ + supplied = { + "model": "noresm", + "atmos_resolution": "ne16", + "experiment": "piControl", + } + del supplied[omitted] + with pytest.raises(TypeError, match=omitted): + build_plan("case", "out", **supplied) + + +# ----------------------------------------------------------- input grid names + + +class TestResolutionPerRealm: + """Tests for choosing each realm's input grid. + + The grid is not a free choice: sea ice is on the model's tripolar grid + whatever the atmosphere was run on, so passing one realm another's grid + would regrid through the wrong weights and yield plausible wrong output. + """ + + @pytest.mark.parametrize("realm", ["atmos", "atmosChem", "aerosol", "land"]) + @pytest.mark.parametrize("atm", ATM_RESOLUTIONS) + def test_atmosphere_realms_take_the_given_grid(self, realm, atm): + """Atmosphere and land use the resolution the case was run at.""" + assert resolution_for("noresm", realm, atm) == atm + + @pytest.mark.parametrize("realm", ["seaIce", "ocean", "ocnBgchem"]) + @pytest.mark.parametrize("model", ["noresm", "cesm"]) + def test_ocean_realms_take_the_model_grid(self, realm, model): + """Ocean and sea ice ignore the atmosphere resolution entirely.""" + assert resolution_for(model, realm, "ne16") == OCEAN_GRID[model] + + def test_unknown_atmosphere_resolution_is_rejected(self): + """A resolution with no weight files fails before any run starts.""" + with pytest.raises(ValueError, match="Unknown atmosphere resolution"): + resolution_for("noresm", "atmos", "ne120") + + @pytest.mark.parametrize("model", ["noresm", "cesm"]) + @pytest.mark.parametrize("atm", ATM_RESOLUTIONS) + def test_landice_is_never_regridded(self, model, atm): + """land ice takes the pass-through grid whatever else was asked for. + + CISM output is written on its native projected grid, georeferenced from + its x/y coordinates, so no weight files apply and the atmosphere's + resolution is irrelevant to it. + """ + assert resolution_for(model, "landIce", atm) == NATIVE_GRID + + def test_cmor_step_carries_the_realm_grid(self, case_dir, tmp_path): + """Each CMOR step is told its own realm's grid, not the atmosphere's.""" + plan = _plan(case_dir, tmp_path / "out", realms=["seaIce", "atmos"]) + sea = _command_of(plan, "cmor-seaIce-mon") + atm = _command_of(plan, "cmor-atmos-mon") + assert sea[sea.index("--resolution") + 1] == "tnx1v4" + assert atm[atm.index("--resolution") + 1] == "ne16" + + def test_landice_runs_every_stage_without_extra_input(self, case_dir, tmp_path): + """land ice needs nothing beyond the atmosphere resolution.""" + plan = _plan(case_dir, tmp_path / "out", realms=["landIce"]) + assert {step.stage for step in plan.steps} == set(STAGES) + assert not plan.skipped + command = _command_of(plan, "cmor-landIce-yr") + assert command[command.index("--resolution") + 1] == NATIVE_GRID + assert command[command.index("--ice-sheet") + 1] == "gris" From 89582134f71d1eb11023387d446baaf0c44f6556 Mon Sep 17 00:00:00 2001 From: mvertens Date: Sun, 4 Oct 2026 22:49:08 +0200 Subject: [PATCH 05/16] now require the driver inputs that were silently defaulted and clarify the help messages for gen_timeseries and cmor_driver to make the required arguments explicit --- scripts/cmor_driver.py | 311 +++++++++++++++----------------- scripts/gen_timeseries.py | 142 +++++++++------ scripts/run_reference_case.py | 61 +++---- src/cmip7_prep/grids.py | 84 +++++++++ src/cmip7_prep/reference_run.py | 80 ++------ tests/test_grids.py | 61 +++++++ tests/test_reference_run.py | 128 +++++++------ 7 files changed, 489 insertions(+), 378 deletions(-) create mode 100644 src/cmip7_prep/grids.py create mode 100644 tests/test_grids.py diff --git a/scripts/cmor_driver.py b/scripts/cmor_driver.py index ae734c5..ea12f71 100755 --- a/scripts/cmor_driver.py +++ b/scripts/cmor_driver.py @@ -40,6 +40,7 @@ patterns_for_variable, ) from cmip7_prep.mapping_compat import Mapping +from cmip7_prep.grids import ATM_RESOLUTIONS, resolution_for from cmip7_prep.regrid import zonal_mean_on_pressure_grid, regrid_to_latlon_ds from cmip7_prep.pipeline import ( realize_regrid_prepare, @@ -84,46 +85,23 @@ }, } + # If CESM archives time series under a component directory, a realm has to be # mapped to the component that wrote it. Several realms share a component. -REALM_COMPONENT_MAP = { - "atmos": "atm", - "aerosol": "atm", - "atmosChem": "atm", - "land": "lnd", - "ocean": "ocn", - "ocnBgchem": "ocn", - "seaIce": "ice", - "landIce": "glc", -} - - def parse_args(): + """Parse the command line, grouped so --help reads as four decisions.""" parser = argparse.ArgumentParser( - description="CMIP7 monthly processing for atm/lnd realms" + description="CMIP7 processing of model time series into CMORized output" ) - parser.add_argument( - "--version", - action="store_true", - help="Show program version and exit", - ) - parser.add_argument( - "--cmip-vars", - nargs="*", - help="List of CMIP variable names to process directly (bypasses variable search)", - ) - parser.add_argument( - "--workers", - type=int, - default=1, - help="Number of Dask workers (default: set to 1 for serial execution)", - ) - parser.add_argument( - "--overwrite", - action="store_true", - help="Overwrite existing timeseries outputs (default: False)", + + required = parser.add_argument_group("required arguments") + required.add_argument( + "--model", + choices=["cesm", "noresm"], + required=True, + help="Model whose variable mappings to use", ) - parser.add_argument( + required.add_argument( "--realm", choices=[ "atmos", @@ -135,125 +113,138 @@ def parse_args(): "seaIce", "landIce", ], - default="atmos", - help="Realm to process. (Default: atmos)", + required=True, + help="Realm to process", ) - parser.add_argument( - "--ice-sheet", - choices=["gris", "ais"], - default=None, + required.add_argument( + "--frequency", + type=str, + required=True, + choices=["mon", "day", "6hr", "3hr", "1hr", "yr"], + help="Frequency of the data being translated", + ) + required.add_argument( + "--tsdir", + type=str, + required=True, help=( - "Ice sheet to process for the landIce realm: 'gris' (Greenland) or " - "'ais' (Antarctica). Required when --realm landIce; ignored otherwise." + "Directory of time series files to CMORize, as produced by " + "gen_timeseries.py or another time series utility" ), ) - parser.add_argument( - "--resolution", + required.add_argument( + "--experiment", type=str, - choices=[ - "ne16", - "ne30", - "tx2_3v2", - "tnx1v4", - "regular", - ], - default="ne30", - help="input_grid name (Default: ne30)", + required=True, + help="CMIP7 experiment_id of the case", ) - parser.add_argument( - "--ocn-static-file", + + case = parser.add_argument_group("other properties of the case") + + case.add_argument( + "--atmos-res", type=str, + choices=list(ATM_RESOLUTIONS), default=None, - help="Path to static file for CESM/MOM variables (optional)", + help=( + "Grid the atmosphere and land were run on. Required when --realm " + "is atmos, atmosChem, aerosol or land; ignored otherwise, since " + "ocean and sea ice are on the model's own tripolar grid and land " + "ice is written on its native projected grid." + ), ) - parser.add_argument( - "--tsdir", + case.add_argument( + "--ocn-static-file", type=str, - help="Time series directory (optional)." - "If not specified, will use a preset time series test directory", + default=None, + help=( + "MOM6 static file supplying the ocean fx fields deptho, areacello " + "and sftof, which come from nowhere else. Used for CESM with " + "--realm ocean; ignored for other models and realms." + ), ) - parser.add_argument( # Move to a wrapper script - "--caseroot", type=str, help="Case root directory" + case.add_argument( + "--ice-sheet", + choices=["gris", "ais"], + default=None, + help=( + "Ice sheet to process for the landIce realm: 'gris' (Greenland) or " + "'ais' (Antarctica). Required when --realm landIce; ignored otherwise." + ), ) - parser.add_argument( # Move to a wrapper script - "--cimeroot", type=str, help="CIME root directory" + case.add_argument( + "--realization-initialization-physics-forcing", + type=str, + default="r1i1p1f1", + help="Variant indices in the format rXiYpZfW (default: r1i1p1f1)", ) - parser.add_argument( - "--test", action="store_true", help="Run in test mode with default paths" + + selection = parser.add_argument_group("which variables to produce") + selection.add_argument( + "--cmip-vars", + nargs="*", + help="CMIP variable names to process directly (bypasses variable search)", ) - parser.add_argument( - "--frequency", - type=str, - default="mon", - choices=["mon", "day", "6hr", "3hr", "1hr", "yr"], + selection.add_argument( + "--run-all-from-yaml", + action="store_true", help=( - "Frequency of data to be translated " - "(mon, day, 6hr, 3hr, 1hr, yr), (Default: mon)" + "Ignore the CMIP7 data request variable list and run every variable " + "defined in the YAML mapping file" ), ) - parser.add_argument( + selection.add_argument( + "--custom-yaml", + default=False, + help="Custom YAML mapping file, overriding the packaged one (optional)", + ) + + paths = parser.add_argument_group("input and output") + + paths.add_argument( "--outdir", type=str, default=".", - help="Output directory for CMORized files. (Default .)", + help="Output directory for CMORized files (default: the current directory)", ) - parser.add_argument( - "--experiment", + paths.add_argument( + "--tables-root", type=str, - default="piControl", - help="Experiment name for data request. (Default piControl)", + default=None, + help="cmip7-cmor-tables checkout, overriding the model-specific default", ) - parser.add_argument( - "--model", - choices=["cesm", "noresm"], - default="cesm", - help="Model to use, default: cesm", + paths.add_argument( + "--overwrite", + action="store_true", + help="Overwrite existing outputs (default: False)", ) - parser.add_argument( + + behaviour = parser.add_argument_group("how to run") + behaviour.add_argument( + "--workers", + type=int, + default=1, + help="Number of Dask workers (default: 1, serial execution)", + ) + behaviour.add_argument( + "--test", action="store_true", help="Run in test mode with default paths" + ) + behaviour.add_argument( "--debug", action="store_true", help="Enable debug logging output", ) - parser.add_argument( + behaviour.add_argument( "--log-level", default="INFO", choices=["DEBUG", "INFO", "WARNING", "ERROR"], - help="log output level", - ) - parser.add_argument( - "--custom-yaml", - default=False, - help="Path to custom YAML mapping file (optional, overrides default packaged YAML)", - ) - parser.add_argument( - "--tables-root", - type=str, - default=None, - help="Path to cmip7-cmor-tables directory (optional, overrides model-specific default)", - ) - parser.add_argument( - "--run-all-from-yaml", - action="store_true", - help="Override CMIP7 data request variable list and run all variables defined in the YAML mapping file", - ) - parser.add_argument( - "--realization-initialization-physics-forcing", - type=str, - default="r1i1p1f1", - help="Realization, initialization, physics, and forcing indices in the format rXiYpZfW (default: r1i1p1f1)", + help="Log output level (default: INFO)", ) args = parser.parse_args() return args -def get_version(): - # Use dynamic version from cmip7_prep - from cmip7_prep import __version__ - - return __version__ - - def _priority_for_logging(data_request, cmip_var) -> str: """Return a stable priority label for logs, including synthetic variables.""" if getattr(cmip_var, "is_synthetic", False): @@ -715,26 +706,49 @@ def main(): # Set variables used below OUTDIR = args.outdir - resolution = args.resolution + # The grid to regrid from follows from the model and realm; see grids.py. + try: + resolution = resolution_for(args.model, args.realm, args.atmos_res) + except ValueError as exc: + logger.error("%s", exc) + sys.exit(2) + logger.info( + "input grid for realm %s is %s (atmosphere/land ran at %s)", + args.realm, + resolution, + args.atmos_res, + ) model = args.model frequency = args.frequency realm = args.realm logger.debug("Realm is %s", realm) ripf_index = args.realization_initialization_physics_forcing - # Ocean fx fields (areacello, deptho, sftof) are read from the MOM6 static - # file and merged into the native data and the CMOR output. CESM only. + # Ocean only: ocean fx fields (deptho, areacello, sftof) are read + # from the MOM6 static file and merged into the native data and + # the CMOR output. ocn_fx_fields = None - if model == "cesm": - if realm in ["ocean", "seaIce"]: - if args.ocn_static_file: - ocn_fx_fields = ocean_fx_fields(args.ocn_static_file) - logger.info( - f"Loaded ocean fx fields from {args.ocn_static_file}: {list(ocn_fx_fields.keys())}" - ) + if model == "cesm" and realm == "ocean": + if args.ocn_static_file: + ocn_fx_fields = ocean_fx_fields(args.ocn_static_file) + logger.info( + "Loaded ocean fx fields from %s: %s", + args.ocn_static_file, + list(ocn_fx_fields), + ) + elif args.ocn_static_file: + logger.warning( + "Ignoring --ocn-static-file: it supplies ocean fx fields, which " + "model=%s realm=%s does not use", + model, + realm, + ) # Determine time series directory (TSDIR) - TSDIR = _resolve_tsdir(args, model, realm) + TSDIR = Path(args.tsdir) + if not TSDIR.is_dir(): + logger.error("Time series directory %s does not exist", TSDIR) + sys.exit(1) # Make output directory if it does not exist OUTDIR = Path(args.outdir) @@ -1026,48 +1040,5 @@ def main(): ) -def _resolve_tsdir(args, model, realm): - """Return the time series directory (TSDIR) for this run.""" - TSDIR = None - if args.tsdir: - TSDIR = Path(args.tsdir) - if not TSDIR.exists(): - logger.error(f"Time series directory {TSDIR} does not exist") - sys.exit(1) - else: - if model == "noresm": - logger.error("must specify --tsdir as an input argument for noresm model") - sys.exit(1) - elif model == "cesm": - if args.caseroot and args.cimeroot: - caseroot = args.caseroot - cimeroot = args.cimeroot - sys.path.append(cimeroot) - _LIBDIR = os.path.join(cimeroot, "CIME", "Tools") - sys.path.append(_LIBDIR) - try: - from CIME.case import Case - except ImportError as e: - logger.error(f"Error importing CIME modules: {e}") - sys.exit(1) - with Case(caseroot, read_only=True) as case: - inputroot = case.get_value("DOUT_S_ROOT") - component = REALM_COMPONENT_MAP.get(realm) - if component is None: - logger.error(f"no time series directory exists for realm {realm}") - sys.exit(1) - TSDIR = ( - Path(inputroot) / component / "proc" / "tseries" / args.frequency - ) - else: - logger.error("no time series directory found for cesm model") - sys.exit(1) - return TSDIR - - if __name__ == "__main__": - args = parse_args() - if getattr(args, "version", False): - print(f"cmor_driver.py version: {get_version()}") - sys.exit(0) main() diff --git a/scripts/gen_timeseries.py b/scripts/gen_timeseries.py index d35902c..a155458 100755 --- a/scripts/gen_timeseries.py +++ b/scripts/gen_timeseries.py @@ -26,7 +26,15 @@ from gents.hfcollection import HFCollection from gents.timeseries import TSCollection -from cmip7_prep.include_patterns import all_include_patterns +from cmip7_prep.include_patterns import all_include_patterns, load_include_patterns + +MODELS = ["cesm", "noresm"] + +# Realms come from the include-pattern tables rather than a list kept here, so +# a realm added to a table is immediately runnable. +DECLARED_REALMS = sorted( + {realm for model in MODELS for realm in load_include_patterns(model)} +) # ++++++++++++++++++++++++++++++ # Unusable history file detection @@ -52,6 +60,19 @@ def time_record_count(path): return -1 +def _time_record_counts(paths, workers): + """Return each path's time record count, in order. + + A pool is only started when there is more than one worker and more than one + file to look at: asking for one worker should not cost a pool, and some + environments refuse to create one at all. + """ + if workers > 1 and len(paths) > 1: + with ProcessPoolExecutor(max_workers=min(workers, len(paths))) as executor: + return list(executor.map(time_record_count, paths)) + return [time_record_count(path) for path in paths] + + def find_unusable_files(inputdir, include_patterns, workers, logger): """ Returns the history files GenTS cannot derive time bounds from: those with @@ -75,14 +96,13 @@ def find_unusable_files(inputdir, include_patterns, workers, logger): return [] unusable = [] - with ProcessPoolExecutor(max_workers=min(workers, len(candidates))) as executor: - for path, count in zip(candidates, executor.map(time_record_count, candidates)): - if count == 0: - logger.warning("Skipping %s: zero time records", path) - unusable.append(path) - elif count == -1: - logger.warning("Skipping %s: cannot be opened", path) - unusable.append(path) + for path, count in zip(candidates, _time_record_counts(candidates, workers)): + if count == 0: + logger.warning("Skipping %s: zero time records", path) + unusable.append(path) + elif count == -1: + logger.warning("Skipping %s: cannot be opened", path) + unusable.append(path) if unusable: logger.warning( @@ -104,43 +124,47 @@ def parse_arguments(): python module and outputs the final argument object. """ - # Create parser object: parser = argparse.ArgumentParser( description="Utility to create time series for all time slice files in a directory" ) - parser.add_argument( - "--debug", action="store_true", help="Turn on debug output (False by default)." - ) - - parser.add_argument( + required = parser.add_argument_group("required arguments") + required.add_argument( "--inputdir", type=str, - help="Comma separated full pathnames of directories containing input spectral element data files (required)", required=True, + help="Full pathname of the directory containing the input history files", ) - parser.add_argument( - "--varlist", - type=str, - help="Comma separated list of variables to process", - default=None, + required.add_argument( + "--model", + choices=MODELS, + required=True, + help=( + "Model whose include patterns to use. The wrong one selects the " + "wrong history streams, so there is no default." + ), ) - parser.add_argument( + required.add_argument( "--realm", - choices=["atmos", "land", "seaIce", "landIce"], - help="Realm to process - sets include patterns for time series (required)", + choices=DECLARED_REALMS, required=True, + help="Realm to process; sets the include patterns for the time series", ) - parser.add_argument( - "--ice-sheet", - choices=["gris", "ais"], + + selection = parser.add_argument_group("selecting what to process") + selection.add_argument( + "--frequency", + nargs="+", default=None, + metavar="FREQ", help=( - "Ice sheet for the landIce realm: 'gris' (Greenland) or 'ais' " - "(Antarctica). Required when --realm landIce; ignored otherwise." + "Only generate time series for these frequencies, e.g. '--frequency " + "6hr' or '--frequency mon day'. Frequencies are those defined for the " + "realm in _include_patterns.yaml. " + "(Default: every frequency the realm defines.)" ), ) - parser.add_argument( + selection.add_argument( "--sampling", choices=["tavg", "tpt"], default=None, @@ -150,54 +174,54 @@ def parse_arguments(): "CMIP7 variable being produced later." ), ) - parser.add_argument( - "--frequency", - nargs="+", + selection.add_argument( + "--varlist", + type=str, + default=None, + help="Comma separated list of variables to process (default: all of them)", + ) + selection.add_argument( + "--ice-sheet", + choices=["gris", "ais"], default=None, - metavar="FREQ", help=( - "Only generate time series for these frequencies, e.g. '--frequency " - "6hr' or '--frequency mon day'. Frequencies are those defined for the " - "realm in _include_patterns.yaml. " - "(Default: every frequency the realm defines.)" + "Ice sheet for the landIce realm: 'gris' (Greenland) or 'ais' " + "(Antarctica). Required when --realm landIce; ignored otherwise." ), ) - parser.add_argument( + selection.add_argument( + "--years-spec", + help="colon separated specification of years to process \n" + " in format of year-first,year-last,year-increments \n " + " where year-increments specifies how many years to user for each time series file \n" + " (default: all files in inputdir are placed in one time series file)", + ) + + output = parser.add_argument_group("output") + output.add_argument( "--outputdir", type=str, - help="Full path to directory where output time series data will be placed (optional) " + help="Full path to directory where output time series data will be placed " "(default: inputdir/../time_series)", ) - parser.add_argument( + output.add_argument( "--overwrite_timeseries", action="store_true", help="Overwrite existing timeseries outputs (default: False)", ) - parser.add_argument( - "--years-spec", - help="colon separated specification of years to process \n" - " in format of year-first,year-last,year-increments \n " - " where year-increments specifies how many years to user for each time series file \n" - " (default: all files in inputdir are placed in one time series file)", - ) - parser.add_argument( + + behaviour = parser.add_argument_group("how to run") + behaviour.add_argument( "--workers", type=int, default=32, help="Number of workers (default: 32)", ) - parser.add_argument( - "--model", - choices=["cesm", "noresm"], - default="cesm", - help="Model to use, default: cesm", + behaviour.add_argument( + "--debug", action="store_true", help="Turn on debug output (False by default)." ) - # Parse Argument inputs - args = parser.parse_args() - - # Error checks - return args + return parser.parse_args() # ++++++++++++++++++++++++++++++ diff --git a/scripts/run_reference_case.py b/scripts/run_reference_case.py index f063ab1..37af09e 100644 --- a/scripts/run_reference_case.py +++ b/scripts/run_reference_case.py @@ -30,13 +30,8 @@ # pylint: disable=wrong-import-position from cmip7_prep.cv_lookup import validate as validate_cv -from cmip7_prep.reference_run import ( - ATM_RESOLUTIONS, - STAGES, - Plan, - Step, - build_plan, -) +from cmip7_prep.grids import ATM_RESOLUTIONS +from cmip7_prep.reference_run import STAGES, Plan, Step, build_plan logger = logging.getLogger("run_reference_case") @@ -73,16 +68,6 @@ def parse_arguments(): "The wrong one produces plausible but wrong output." ), ) - required.add_argument( - "--atmos-resolution", - choices=list(ATM_RESOLUTIONS), - required=True, - help=( - "Grid the atmosphere and land were run on. Ocean and sea ice are " - "always on the model's own tripolar grid, and land ice is written " - "on its native projected grid, so those are derived from this." - ), - ) required.add_argument( "--experiment", required=True, @@ -93,6 +78,16 @@ def parse_arguments(): ) selection = parser.add_argument_group("selecting what to run") + selection.add_argument( + "--atmos-res", + choices=list(ATM_RESOLUTIONS), + default=None, + help=( + "Grid the atmosphere and land were run on. Required unless the " + "realms being processed are all ocean, sea ice or land ice, which " + "derive their own grid." + ), + ) selection.add_argument( "--realms", nargs="+", @@ -264,20 +259,24 @@ def main(): logger.error("%s", exc) sys.exit(2) - plan = build_plan( - case_dir=case_dir, - outdir=args.outdir, - model=args.model, - realms=args.realms, - frequencies=args.frequencies, - years=args.years, - stages=args.stages, - atmos_resolution=args.atmos_resolution, - experiment=args.experiment, - workers=args.workers, - ice_sheet=args.ice_sheet, - scripts_dir=_LOCAL_PATH, - ) + try: + plan = build_plan( + case_dir=case_dir, + outdir=args.outdir, + model=args.model, + realms=args.realms, + frequencies=args.frequencies, + years=args.years, + stages=args.stages, + atmos_res=args.atmos_res, + experiment=args.experiment, + workers=args.workers, + ice_sheet=args.ice_sheet, + scripts_dir=_LOCAL_PATH, + ) + except ValueError as exc: + logger.error("%s", exc) + sys.exit(2) if not plan.steps: logger.error("Nothing to do. Skipped: %s", plan.skipped or "(nothing)") diff --git a/src/cmip7_prep/grids.py b/src/cmip7_prep/grids.py new file mode 100644 index 0000000..c07476c --- /dev/null +++ b/src/cmip7_prep/grids.py @@ -0,0 +1,84 @@ +"""The grid each realm's data is on. + +Model components run on different grids, so the grid to regrid from depends on +which realm is being processed: + + atmos, atmosChem, aerosol, land the atmosphere/land grid the case was + run at: ne30, ne16, or 'regular' when + the output already carries lat/lon + ocean, ocnBgchem, seaIce the model's tripolar grid: tnx1v4 for + NorESM, tx2_3v2 for CESM + landIce its own projected grid; CISM output is + georeferenced from the x/y coordinates + in the files, not regridded + +:func:`resolution_for` returns it, so a caller gives only the atmosphere/land +resolution -- and only when processing a realm that uses it. +""" + +from __future__ import annotations + +# Grids the atmosphere and land may be run on. 'regular' means the output +# already carries lat/lon and is not regridded. +ATM_RESOLUTIONS = ("ne30", "ne16", "regular") + +# Realms whose input grid is the atmosphere/land grid the case was run on. +ATM_GRID_REALMS = frozenset({"atmos", "atmosChem", "aerosol", "land"}) + +# The ocean/sea-ice grid each model writes on. These realms never share the +# atmosphere's grid, so passing one an atmosphere resolution would regrid +# through the wrong weights. +OCEAN_GRID = {"noresm": "tnx1v4", "cesm": "tx2_3v2"} + +OCEAN_GRID_REALMS = frozenset({"ocean", "ocnBgchem", "seaIce"}) + +# Realms written on their native grid, with no ESMF regridding. CISM land-ice +# output carries projected x/y coordinates, which cmor_writer georeferences with +# the ice sheet's own map projection (see cism_grid.project_xy_to_latlon), so no +# weight files are involved. The value below is a pass-through: it names the +# entry that builds a regridder and discards it, which is what the unregridded +# path expects. +NATIVE_GRID_REALMS = frozenset({"landIce"}) +NATIVE_GRID = "regular" + + +def needs_atmos_res(realm: str) -> bool: + """Return whether a realm's grid depends on the atmosphere resolution. + + Only the realms sharing the atmosphere and land grid do. Sea ice, ocean + and land ice have their own, so asking for an atmosphere resolution to + process them would be asking for something that cannot matter. + """ + return realm in ATM_GRID_REALMS + + +def resolution_for(model: str, realm: str, atmos_res: str | None = None) -> str: + """Return the input grid name for one realm. + + ``atmos_res`` is the grid the atmosphere and land were run on. It is + required only for the realms on that grid, and ignored for the rest, so a + sea-ice run need not supply one. + """ + if realm in ATM_GRID_REALMS: + if atmos_res is None: + raise ValueError( + f"Realm {realm!r} is on the atmosphere/land grid, so its " + "resolution must be given" + ) + if atmos_res not in ATM_RESOLUTIONS: + raise ValueError( + f"Unknown atmosphere resolution {atmos_res!r}; " + f"choose from {list(ATM_RESOLUTIONS)}" + ) + return atmos_res + if realm in OCEAN_GRID_REALMS: + try: + return OCEAN_GRID[model] + except KeyError: + raise ValueError( + f"No ocean grid known for model={model!r}; " + f"known models: {sorted(OCEAN_GRID)}" + ) from None + if realm in NATIVE_GRID_REALMS: + return NATIVE_GRID + raise ValueError(f"No input grid known for realm {realm!r}") diff --git a/src/cmip7_prep/reference_run.py b/src/cmip7_prep/reference_run.py index cf325ce..56f7466 100644 --- a/src/cmip7_prep/reference_run.py +++ b/src/cmip7_prep/reference_run.py @@ -21,6 +21,7 @@ from pathlib import Path from typing import Sequence +from .grids import ATM_RESOLUTIONS, needs_atmos_res from .include_patterns import load_include_patterns # Which component directory holds each realm's history files. Mirrors @@ -39,29 +40,6 @@ # The ice sheet to process when a realm is per-ice-sheet and none was named. DEFAULT_ICE_SHEET = "gris" -# Grids the atmosphere and land may be run on. 'regular' means the output -# already carries lat/lon and is not regridded. -ATM_RESOLUTIONS = ("ne30", "ne16", "regular") - -# Realms whose input grid is the atmosphere/land grid the case was run on. -ATM_GRID_REALMS = frozenset({"atmos", "atmosChem", "aerosol", "land"}) - -# Realms written on their native grid, with no ESMF regridding. CISM land-ice -# output carries projected x/y coordinates which cmor_writer georeferences with -# the ice sheet's own map projection (see cism_grid.project_xy_to_latlon), so no -# weight files are involved and the resolution would be ignored. The value -# below is a pass-through: it names the entry that builds a regridder and -# discards it, which is what the unregridded path expects. -NATIVE_GRID_REALMS = frozenset({"landIce"}) -NATIVE_GRID = "regular" - -# The ocean/sea-ice grid each model writes on. These realms never share the -# atmosphere's grid, so passing one an atmosphere resolution would regrid -# through the wrong weights and produce plausible but wrong output. -OCEAN_GRID = {"noresm": "tnx1v4", "cesm": "tx2_3v2"} - -OCEAN_GRID_REALMS = frozenset({"ocean", "ocnBgchem", "seaIce"}) - STAGES = ("timeseries", "cmor", "validate") @@ -109,35 +87,6 @@ def frequencies_for(model: str, realm: str) -> list[str]: return sorted(load_include_patterns(model)[realm]) -def resolution_for(model: str, realm: str, atm_resolution: str) -> str: - """Return the input grid name for one realm. - - Every realm's grid is derived, so no caller has to know the model's grid - layout. ``atm_resolution`` is the grid the atmosphere and land were run on - and applies only to those realms; ocean and sea ice are always on the - model's tripolar grid; land ice is written on its native projected grid and - is not regridded at all. - """ - if atm_resolution not in ATM_RESOLUTIONS: - raise ValueError( - f"Unknown atmosphere resolution {atm_resolution!r}; " - f"choose from {list(ATM_RESOLUTIONS)}" - ) - if realm in ATM_GRID_REALMS: - return atm_resolution - if realm in OCEAN_GRID_REALMS: - try: - return OCEAN_GRID[model] - except KeyError: - raise ValueError( - f"No ocean grid known for model={model!r}; " - f"known models: {sorted(OCEAN_GRID)}" - ) from None - if realm in NATIVE_GRID_REALMS: - return NATIVE_GRID - raise ValueError(f"No input grid known for realm {realm!r}") - - def history_dir(case_dir: os.PathLike | str, realm: str) -> Path: """Return the history directory holding one realm's output.""" component = REALM_COMPONENT[realm] @@ -153,7 +102,7 @@ def build_plan( frequencies: Sequence[str] | None = None, years: str | None = None, stages: Sequence[str] = STAGES, - atmos_resolution: str, + atmos_res: str | None = None, experiment: str, workers: int = 4, ice_sheet: str | None = None, @@ -165,7 +114,7 @@ def build_plan( history directory is absent is recorded in ``Plan.skipped`` rather than failing the run, since an archived case need not hold every component. - ``atmos_resolution`` is the grid the atmosphere and land were run on; every + ``atmos_res`` is the grid the atmosphere and land were run on; every other realm's grid is derived from it or from the model. ``years`` is passed through to gen_timeseries.py as ``--years-spec`` and so @@ -174,6 +123,18 @@ def build_plan( unknown = [stage for stage in stages if stage not in STAGES] if unknown: raise ValueError(f"Unknown stage(s) {unknown}; choose from {list(STAGES)}") + wanted_realms = list(realms) if realms else realms_for(model) + on_atmos_grid = [realm for realm in wanted_realms if needs_atmos_res(realm)] + if on_atmos_grid and atmos_res is None: + raise ValueError( + "An atmosphere resolution is needed for " + f"{on_atmos_grid}; the other realms derive their own grid" + ) + if atmos_res is not None and atmos_res not in ATM_RESOLUTIONS: + raise ValueError( + f"Unknown atmosphere resolution {atmos_res!r}; " + f"choose from {list(ATM_RESOLUTIONS)}" + ) case_dir = Path(case_dir) outdir = Path(outdir) @@ -182,7 +143,6 @@ def build_plan( cmor_root = outdir / "cmor" plan = Plan(timeseries_root=timeseries_root, cmor_root=cmor_root) - wanted_realms = list(realms) if realms else realms_for(model) for realm in wanted_realms: if realm not in REALM_COMPONENT: @@ -205,8 +165,6 @@ def build_plan( sheet = ice_sheet or DEFAULT_ICE_SHEET if realm == "landIce" else None ts_dir = timeseries_root / realm - realm_resolution = resolution_for(model, realm, atmos_resolution) - if "timeseries" in stages: plan.steps.append( _timeseries_step( @@ -232,7 +190,7 @@ def build_plan( ts_dir=ts_dir, cmor_root=cmor_root, model=model, - resolution=realm_resolution, + atmos_res=(atmos_res if needs_atmos_res(realm) else None), experiment=experiment, workers=workers, sheet=sheet, @@ -300,7 +258,7 @@ def _cmor_step( ts_dir, cmor_root, model, - resolution, + atmos_res, experiment, workers, sheet, @@ -319,13 +277,13 @@ def _cmor_step( str(cmor_root), "--model", model, - "--resolution", - resolution, "--experiment", experiment, "--workers", str(workers), ] + if atmos_res: + command += ["--atmos-res", atmos_res] if sheet: command += ["--ice-sheet", sheet] return Step( diff --git a/tests/test_grids.py b/tests/test_grids.py new file mode 100644 index 0000000..c1134c3 --- /dev/null +++ b/tests/test_grids.py @@ -0,0 +1,61 @@ +"""Tests for deriving each realm's input grid.""" + +import pytest + +from cmip7_prep.grids import ( + ATM_RESOLUTIONS, + needs_atmos_res, + NATIVE_GRID, + OCEAN_GRID, + resolution_for, +) + + +class TestResolutionPerRealm: + """Tests for choosing each realm's input grid. + + The grid is not a free choice: sea ice is on the model's tripolar grid + whatever the atmosphere was run on, so passing one realm another's grid + would regrid through the wrong weights and yield plausible wrong output. + """ + + @pytest.mark.parametrize("realm", ["atmos", "atmosChem", "aerosol", "land"]) + @pytest.mark.parametrize("atm", ATM_RESOLUTIONS) + def test_atmosphere_realms_take_the_given_grid(self, realm, atm): + """Atmosphere and land use the resolution the case was run at.""" + assert resolution_for("noresm", realm, atm) == atm + + @pytest.mark.parametrize("realm", ["seaIce", "ocean", "ocnBgchem"]) + @pytest.mark.parametrize("model", ["noresm", "cesm"]) + def test_ocean_realms_take_the_model_grid(self, realm, model): + """Ocean and sea ice ignore the atmosphere resolution entirely.""" + assert resolution_for(model, realm, "ne16") == OCEAN_GRID[model] + + @pytest.mark.parametrize("realm", ["seaIce", "ocean", "ocnBgchem", "landIce"]) + def test_realms_with_their_own_grid_need_no_resolution(self, realm): + """A sea-ice or land-ice run need not supply an atmosphere resolution.""" + assert resolution_for("noresm", realm) + assert not needs_atmos_res(realm) + + @pytest.mark.parametrize("realm", ["atmos", "atmosChem", "aerosol", "land"]) + def test_atmosphere_realms_say_so_when_it_is_missing(self, realm): + """A realm on the atmosphere grid refuses to guess its resolution.""" + assert needs_atmos_res(realm) + with pytest.raises(ValueError, match="resolution must be given"): + resolution_for("noresm", realm) + + def test_unknown_atmosphere_resolution_is_rejected(self): + """A resolution with no weight files fails before any run starts.""" + with pytest.raises(ValueError, match="Unknown atmosphere resolution"): + resolution_for("noresm", "atmos", "ne120") + + @pytest.mark.parametrize("model", ["noresm", "cesm"]) + @pytest.mark.parametrize("atm", ATM_RESOLUTIONS) + def test_landice_is_never_regridded(self, model, atm): + """land ice takes the pass-through grid whatever else was asked for. + + CISM output is written on its native projected grid, georeferenced from + its x/y coordinates, so no weight files apply and the atmosphere's + resolution is irrelevant to it. + """ + assert resolution_for(model, "landIce", atm) == NATIVE_GRID diff --git a/tests/test_reference_run.py b/tests/test_reference_run.py index 424467e..6de6a5f 100644 --- a/tests/test_reference_run.py +++ b/tests/test_reference_run.py @@ -3,16 +3,12 @@ import pytest from cmip7_prep.reference_run import ( - ATM_RESOLUTIONS, - NATIVE_GRID, - OCEAN_GRID, REALM_COMPONENT, STAGES, build_plan, frequencies_for, history_dir, realms_for, - resolution_for, ) @@ -31,7 +27,7 @@ def _plan(case_dir, outdir, **kwargs): wrong output, so a caller must say which it means. """ kwargs.setdefault("model", "noresm") - kwargs.setdefault("atmos_resolution", "ne16") + kwargs.setdefault("atmos_res", "ne16") kwargs.setdefault("experiment", "piControl") return build_plan(case_dir, outdir, **kwargs) @@ -223,20 +219,21 @@ def test_validate_is_strict(self, case_dir, tmp_path): assert "--strict" in _command_of(plan, "validate-seaIce-mon") def test_resolution_and_experiment_reach_cmor(self, case_dir, tmp_path): - """Resolution and experiment are passed to the CMOR step. + """The atmosphere resolution and experiment are passed to the CMOR step. - Resolution is checked on an atmosphere realm, since that is the only - kind that takes the value given; see TestResolutionPerRealm. + The CMOR step is given the atmosphere resolution rather than the grid + it will regrid from: cmor_driver.py derives that itself from --model + and --realm, so the derivation has one home. See test_grids.py. """ plan = _plan( case_dir, tmp_path / "out", realms=["atmos"], - atmos_resolution="ne30", + atmos_res="ne30", experiment="historical", ) command = _command_of(plan, "cmor-atmos-mon") - assert command[command.index("--resolution") + 1] == "ne30" + assert command[command.index("--atmos-res") + 1] == "ne30" assert command[command.index("--experiment") + 1] == "historical" def test_every_step_has_a_distinct_log_name(self, case_dir, tmp_path): @@ -246,16 +243,18 @@ def test_every_step_has_a_distinct_log_name(self, case_dir, tmp_path): assert len(names) == len(set(names)) -@pytest.mark.parametrize("omitted", ["model", "atmos_resolution", "experiment"]) +@pytest.mark.parametrize("omitted", ["model", "experiment"]) def test_case_properties_are_required(omitted): """build_plan refuses to guess a property of the case. - Model, resolution and experiment each change what the output means while - leaving it looking valid, so none of them may be defaulted. + Model and experiment change what the output means while leaving it looking + valid, so neither may be defaulted. The atmosphere resolution is not here + because it is only needed by some realms; see + TestAtmosResolutionIsConditional. """ supplied = { "model": "noresm", - "atmos_resolution": "ne16", + "atmos_res": "ne16", "experiment": "piControl", } del supplied[omitted] @@ -263,58 +262,73 @@ def test_case_properties_are_required(omitted): build_plan("case", "out", **supplied) -# ----------------------------------------------------------- input grid names +class TestAtmosResolutionIsConditional: + """The atmosphere resolution is asked for only when a realm uses it.""" + def test_not_needed_for_realms_with_their_own_grid(self, case_dir, tmp_path): + """Sea ice and land ice plan fully without one.""" + plan = build_plan( + case_dir, + tmp_path / "out", + model="noresm", + experiment="piControl", + realms=["seaIce", "landIce"], + ) + assert plan.steps + for step in plan.for_stage("cmor"): + assert "--atmos-res" not in step.command + + @pytest.mark.parametrize("realm", ["atmos", "land"]) + def test_needed_for_realms_on_the_atmosphere_grid(self, case_dir, tmp_path, realm): + """A realm on that grid refuses to plan without one, and names itself.""" + with pytest.raises(ValueError, match=f"needed for \\['{realm}'\\]"): + build_plan( + case_dir, + tmp_path / "out", + model="noresm", + experiment="piControl", + realms=[realm], + ) + + def test_the_whole_matrix_names_every_realm_that_needs_it(self, case_dir, tmp_path): + """Planning every realm without one lists the four that require it.""" + with pytest.raises(ValueError, match="atmos.*atmosChem.*aerosol.*land"): + build_plan( + case_dir, tmp_path / "out", model="noresm", experiment="piControl" + ) + + def test_only_atmosphere_realms_carry_the_flag(self, case_dir, tmp_path): + """With one supplied, only the realms that use it are given it.""" + plan = _plan(case_dir, tmp_path / "out", realms=["atmos", "seaIce"]) + atm = _command_of(plan, "cmor-atmos-mon") + sea = _command_of(plan, "cmor-seaIce-mon") + assert atm[atm.index("--atmos-res") + 1] == "ne16" + assert "--atmos-res" not in sea -class TestResolutionPerRealm: - """Tests for choosing each realm's input grid. - The grid is not a free choice: sea ice is on the model's tripolar grid - whatever the atmosphere was run on, so passing one realm another's grid - would regrid through the wrong weights and yield plausible wrong output. - """ +# ------------------------------------------- forwarding the grid decision - @pytest.mark.parametrize("realm", ["atmos", "atmosChem", "aerosol", "land"]) - @pytest.mark.parametrize("atm", ATM_RESOLUTIONS) - def test_atmosphere_realms_take_the_given_grid(self, realm, atm): - """Atmosphere and land use the resolution the case was run at.""" - assert resolution_for("noresm", realm, atm) == atm - - @pytest.mark.parametrize("realm", ["seaIce", "ocean", "ocnBgchem"]) - @pytest.mark.parametrize("model", ["noresm", "cesm"]) - def test_ocean_realms_take_the_model_grid(self, realm, model): - """Ocean and sea ice ignore the atmosphere resolution entirely.""" - assert resolution_for(model, realm, "ne16") == OCEAN_GRID[model] - - def test_unknown_atmosphere_resolution_is_rejected(self): - """A resolution with no weight files fails before any run starts.""" - with pytest.raises(ValueError, match="Unknown atmosphere resolution"): - resolution_for("noresm", "atmos", "ne120") - - @pytest.mark.parametrize("model", ["noresm", "cesm"]) - @pytest.mark.parametrize("atm", ATM_RESOLUTIONS) - def test_landice_is_never_regridded(self, model, atm): - """land ice takes the pass-through grid whatever else was asked for. - - CISM output is written on its native projected grid, georeferenced from - its x/y coordinates, so no weight files apply and the atmosphere's - resolution is irrelevant to it. - """ - assert resolution_for(model, "landIce", atm) == NATIVE_GRID - def test_cmor_step_carries_the_realm_grid(self, case_dir, tmp_path): - """Each CMOR step is told its own realm's grid, not the atmosphere's.""" - plan = _plan(case_dir, tmp_path / "out", realms=["seaIce", "atmos"]) - sea = _command_of(plan, "cmor-seaIce-mon") - atm = _command_of(plan, "cmor-atmos-mon") - assert sea[sea.index("--resolution") + 1] == "tnx1v4" - assert atm[atm.index("--resolution") + 1] == "ne16" +class TestGridIsNotDecidedHere: + """The plan never decides which grid a run regrids from.""" + + def test_no_step_carries_a_derived_grid(self, case_dir, tmp_path): + """No command names an input grid such as tnx1v4. + + cmor_driver.py derives that from --model and --realm, so a plan that + passed one would mean two places deciding the same thing. + """ + plan = _plan(case_dir, tmp_path / "out") + for step in plan.steps: + assert "--resolution" not in step.command + assert "tnx1v4" not in step.command + assert "tx2_3v2" not in step.command def test_landice_runs_every_stage_without_extra_input(self, case_dir, tmp_path): - """land ice needs nothing beyond the atmosphere resolution.""" + """land ice needs no grid input of its own.""" plan = _plan(case_dir, tmp_path / "out", realms=["landIce"]) assert {step.stage for step in plan.steps} == set(STAGES) assert not plan.skipped command = _command_of(plan, "cmor-landIce-yr") - assert command[command.index("--resolution") + 1] == NATIVE_GRID + assert "--atmos-res" not in command assert command[command.index("--ice-sheet") + 1] == "gris" From 2f0aabf6808e0994c42c67d22f6834e1c43ba3d8 Mon Sep 17 00:00:00 2001 From: mvertens Date: Mon, 5 Oct 2026 11:23:23 +0200 Subject: [PATCH 06/16] better error handling in gen_timeseries.py --- scripts/gen_timeseries.py | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/scripts/gen_timeseries.py b/scripts/gen_timeseries.py index a155458..2d723ee 100755 --- a/scripts/gen_timeseries.py +++ b/scripts/gen_timeseries.py @@ -282,21 +282,30 @@ def main(): # Determine number of files used in time series creation cnt = 0 filtered = [] + all_patterns = list(include_patterns) for include_pattern in include_patterns: num = len(glob.glob(os.path.join(inputdir, include_pattern))) - logger.info(f"include pattern {include_pattern} has num {num}") if num == 0: - logger.info(f"removing {include_pattern}") + logger.info( + "No files match %s; this stream was not written, skipping it", + include_pattern, + ) else: cnt += num - logger.info(f"Processing {num} files with {include_pattern}") + logger.info("Processing %d file(s) matching %s", num, include_pattern) filtered.append(include_pattern) include_patterns = filtered if cnt == 0: - logger.warning( - f"No input files to process in {inputdir} with {include_patterns}" + # Being asked to process a realm and finding nothing is a failure, not + # a no-op: the caller named this realm, so an empty result means the + # input directory, the realm or the frequencies were wrong. + logger.error( + "No input files in %s for realm %s with patterns %s", + inputdir, + args.realm, + all_patterns, ) - sys.exit(0) + sys.exit(1) logger.info(f"include patterns are {include_patterns}") varlist = ( From 103cf69dcdf122e4a61c214a7dd4ac62a1ccc4ed Mon Sep 17 00:00:00 2001 From: mvertens Date: Mon, 5 Oct 2026 12:33:51 +0200 Subject: [PATCH 07/16] make an interrupted end-to-end run resumable --- scripts/cmor_driver.py | 13 ++- scripts/gen_timeseries.py | 63 +++++++++++--- scripts/run_reference_case.py | 85 ++++++++++++++++-- src/cmip7_prep/reference_run.py | 115 +++++++++++++++++++++--- tests/test_reference_run.py | 150 +++++++++++++++++++++++++++++++- 5 files changed, 392 insertions(+), 34 deletions(-) diff --git a/scripts/cmor_driver.py b/scripts/cmor_driver.py index ea12f71..52cd249 100755 --- a/scripts/cmor_driver.py +++ b/scripts/cmor_driver.py @@ -173,10 +173,19 @@ def parse_args(): ), ) case.add_argument( + "--variant-label", "--realization-initialization-physics-forcing", + dest="variant_label", type=str, default="r1i1p1f1", - help="Variant indices in the format rXiYpZfW (default: r1i1p1f1)", + metavar="rXiYpZfW", + help=( + "Ensemble member this output belongs to: realization, " + "initialization, physics and forcing indices, e.g. r1i1p1f1. This " + "is the CMIP variant_label, and appears in the output paths and " + "file names. The longer spelling is accepted as an alias. " + "(default: r1i1p1f1)" + ), ) selection = parser.add_argument_group("which variables to produce") @@ -722,7 +731,7 @@ def main(): frequency = args.frequency realm = args.realm logger.debug("Realm is %s", realm) - ripf_index = args.realization_initialization_physics_forcing + ripf_index = args.variant_label # Ocean only: ocean fx fields (deptho, areacello, sftof) are read # from the MOM6 static file and merged into the native data and diff --git a/scripts/gen_timeseries.py b/scripts/gen_timeseries.py index 2d723ee..8c9243c 100755 --- a/scripts/gen_timeseries.py +++ b/scripts/gen_timeseries.py @@ -73,6 +73,37 @@ def _time_record_counts(paths, workers): return [time_record_count(path) for path in paths] +def _overwrite_or_resume(ts_collection, overwrite, logger): + """Return the collection set either to overwrite output or to resume. + + Asked to overwrite, GenTS rewrites every time series. Otherwise the + collection is narrowed to the timesteps no output covers yet, so a rerun + after an interrupted run continues instead of repeating work. + """ + if overwrite: + logger.info("Overwriting any existing time series") + return ts_collection.apply_overwrite("*") + resumed = ts_collection.skip_existing() + logger.info( + "Skipping time series that already exist; pass --overwrite_timeseries " + "to rewrite them" + ) + return resumed + + +def _execute_or_report_complete(ts_collection, logger): + """Generate the time series, or report that there is nothing left to do. + + Narrowing a finished run to its missing timesteps leaves nothing, and GenTS + treats an empty collection as an error. Having already produced everything + asked for is a success, so it is reported as one. + """ + if len(ts_collection) == 0: + logger.info("Every requested time series already exists; nothing to generate") + return + ts_collection.execute() + + def find_unusable_files(inputdir, include_patterns, workers, logger): """ Returns the history files GenTS cannot derive time bounds from: those with @@ -207,7 +238,11 @@ def parse_arguments(): output.add_argument( "--overwrite_timeseries", action="store_true", - help="Overwrite existing timeseries outputs (default: False)", + help=( + "Rewrite time series that already exist. Without this, existing " + "output is kept and only the missing timesteps are produced, so " + "an interrupted run can simply be repeated." + ), ) behaviour = parser.add_argument_group("how to run") @@ -331,14 +366,19 @@ def main(): # Create base TSCollection logger.info("Starting ts_collection") ts_collection = TSCollection(hf_collection, outputdir, num_processes=workers) - ts_collection = ts_collection.apply_overwrite("*") if varlist: ts_collection = ts_collection.include("*", var_glob=varlist) + # Empty here means nothing matched at all, which is a mistake worth + # stopping for. Empty only after narrowing below means the work was + # already done, which is not. if len(ts_collection) == 0: raise RuntimeError( "No matching variables/files found for time series generation" ) - ts_collection.execute() + ts_collection = _overwrite_or_resume( + ts_collection, args.overwrite_timeseries, logger + ) + _execute_or_report_complete(ts_collection, logger) logger.info("Finished ts_collection") else: @@ -377,11 +417,6 @@ def main(): ) logger.info("Finished ts_collection") - # Apply overwrite if requested: - # If --overwrite flag was passed, tells GenTS to overwrite existing time series files - if args.overwrite_timeseries: - ts_collection = ts_collection.apply_overwrite("*") - # Perform the time series generation for this pattern if varlist: ts_collection = ts_collection.include("*", var_glob=varlist) @@ -390,12 +425,16 @@ def main(): "No matching variables/files found for time series generation" ) - logger.info( - "Variables scheduled: %s", - [order["primary_var"] for order in ts_collection], + ts_collection = _overwrite_or_resume( + ts_collection, args.overwrite_timeseries, logger ) + if len(ts_collection): + logger.info( + "Variables scheduled: %s", + [order["primary_var"] for order in ts_collection], + ) - ts_collection.execute() + _execute_or_report_complete(ts_collection, logger) logger.info("Timeseries processing complete") diff --git a/scripts/run_reference_case.py b/scripts/run_reference_case.py index 37af09e..c39138b 100644 --- a/scripts/run_reference_case.py +++ b/scripts/run_reference_case.py @@ -9,6 +9,10 @@ frequencies, and prints one table at the end saying what succeeded, what failed, and how much was produced. +Steps that succeed are recorded in the output directory, so repeating the +command after an interruption continues rather than starting over; --force +runs everything again. + It is meant to be run on the machine holding the data, interactively -- there is no batch submission here. A full run takes a while, so start it under ``tmux`` or ``screen`` if the connection might drop. Use ``--dry-run`` first to @@ -19,6 +23,7 @@ """ import argparse +import json import logging import subprocess import sys @@ -117,6 +122,17 @@ def parse_arguments(): default=list(STAGES), help="Stages to run (default: all three)", ) + selection.add_argument( + "--variant-label", + default=None, + metavar="rXiYpZfW", + help=( + "Ensemble member to label the output with, e.g. r1i1p1f1 " + "(cmor_driver.py's --realization-initialization-physics-forcing). " + "Needed when producing output for submission rather than just " + "exercising the chain, since it distinguishes members." + ), + ) selection.add_argument( "--ice-sheet", choices=["gris", "ais"], @@ -130,8 +146,18 @@ def parse_arguments(): type=int, default=4, help=( - "Worker processes per stage (default: 4, chosen to be polite on a " - "shared interactive node rather than to be fast)" + "Worker processes for time series generation (default: 4, chosen " + "to be polite on a shared interactive node). CMORization always " + "runs with one worker." + ), + ) + behaviour.add_argument( + "--force", + action="store_true", + help=( + "Run every step again, including those a previous run in this " + "output directory already completed, and rewrite time series " + "that already exist" ), ) behaviour.add_argument( @@ -148,14 +174,41 @@ def parse_arguments(): "--tables-root", default=None, help=( - "cmip7-cmor-tables checkout, used to check --experiment against the " - "controlled vocabulary (default: the one beside this repo)" + "cmip7-cmor-tables checkout. Used to check --experiment against " + "the controlled vocabulary before starting, and passed on to " + "CMORization (default: the one beside this repo)" ), ) behaviour.add_argument("--log-level", default="INFO", help="Default: INFO") return parser.parse_args() +# Steps that finished are recorded here so a repeated run continues instead of +# starting over. One file per output directory, next to the logs. +COMPLETED_NAME = "completed_steps.json" + + +def read_completed(outdir: Path) -> set[str]: + """Return the keys of steps a previous run completed in this directory.""" + path = outdir / COMPLETED_NAME + if not path.is_file(): + return set() + try: + return set(json.loads(path.read_text(encoding="utf-8"))) + except (ValueError, OSError) as exc: + logger.warning("Ignoring %s: %s", path, exc) + return set() + + +def record_completed(outdir: Path, done: set[str]) -> None: + """Write the set of completed step keys, sorted so the file diffs cleanly.""" + path = outdir / COMPLETED_NAME + try: + path.write_text(json.dumps(sorted(done), indent=2) + "\n", encoding="utf-8") + except OSError as exc: + logger.warning("Could not record progress in %s: %s", path, exc) + + def run_step(step: Step, log_dir: Path) -> dict: """Run one step, tee its output to a log file, and return the outcome.""" log_path = log_dir / step.log_name @@ -272,6 +325,9 @@ def main(): experiment=args.experiment, workers=args.workers, ice_sheet=args.ice_sheet, + overwrite_timeseries=args.force, + variant_label=args.variant_label, + tables_root=args.tables_root, scripts_dir=_LOCAL_PATH, ) except ValueError as exc: @@ -296,14 +352,33 @@ def main(): log_dir = outdir / "logs" log_dir.mkdir(parents=True, exist_ok=True) + completed = set() if args.force else read_completed(outdir) + if completed: + logger.info( + "%d step(s) already completed in %s; pass --force to run them again", + len(completed), + outdir, + ) + outcomes = [] for step in plan.steps: + if step.key in completed: + logger.info("[%s] already done, skipping", step.key) + continue outcome = run_step(step, log_dir) outcomes.append(outcome) - if outcome["returncode"] != 0 and not args.keep_going: + if outcome["returncode"] == 0: + completed.add(step.key) + # Written after every step, so an interrupted run still knows what + # it finished. + record_completed(outdir, completed) + elif not args.keep_going: logger.error("Stopping at the first failure; pass --keep-going to continue") break + if not outcomes: + logger.info("Nothing left to do: every planned step was already completed") + return print_report(outcomes, plan) failed = any(outcome["returncode"] != 0 for outcome in outcomes) diff --git a/src/cmip7_prep/reference_run.py b/src/cmip7_prep/reference_run.py index 56f7466..b839c03 100644 --- a/src/cmip7_prep/reference_run.py +++ b/src/cmip7_prep/reference_run.py @@ -16,13 +16,15 @@ from __future__ import annotations +import glob import os +import sys from dataclasses import dataclass, field from pathlib import Path from typing import Sequence from .grids import ATM_RESOLUTIONS, needs_atmos_res -from .include_patterns import load_include_patterns +from .include_patterns import all_include_patterns, load_include_patterns # Which component directory holds each realm's history files. Mirrors # REALM_COMPONENT_MAP in scripts/cmor_driver.py. @@ -42,6 +44,17 @@ STAGES = ("timeseries", "cmor", "validate") +# The steps run with the interpreter running this, not whatever "python" means +# on PATH: that keeps every stage in the environment the run was started from, +# and works when the runner was invoked as python3 or by absolute path. +PYTHON = sys.executable or "python" + +# CMORization is run with one worker: running it in parallel has not been +# tested. Its --workers are Dask workers in any case, which is a different +# thing from the processes gen_timeseries.py spreads file reads over, so the +# two are not set from one number. +CMOR_WORKERS = 1 + @dataclass class Step: @@ -87,6 +100,32 @@ def frequencies_for(model: str, realm: str) -> list[str]: return sorted(load_include_patterns(model)[realm]) +def written_frequencies( + case_dir: os.PathLike | str, + model: str, + realm: str, + frequencies: Sequence[str], + ice_sheet: str | None = None, +) -> list[str]: + """Return those frequencies the case actually wrote history files for. + + A table declares every frequency a realm *can* produce, but a given run + writes only the streams it was configured for. Planning a CMORization step + for a frequency with no history files behind it guarantees a failure that + says nothing about the code, so those frequencies are dropped from the plan + instead. + """ + inputdir = history_dir(case_dir, realm) + written = [] + for frequency in frequencies: + patterns = all_include_patterns(model, realm, ice_sheet, None, [frequency]) + if any( + glob.glob(os.path.join(inputdir, f"*{pattern}*")) for pattern in patterns + ): + written.append(frequency) + return written + + def history_dir(case_dir: os.PathLike | str, realm: str) -> Path: """Return the history directory holding one realm's output.""" component = REALM_COMPONENT[realm] @@ -106,6 +145,9 @@ def build_plan( experiment: str, workers: int = 4, ice_sheet: str | None = None, + overwrite_timeseries: bool = False, + variant_label: str | None = None, + tables_root: os.PathLike | str | None = None, scripts_dir: os.PathLike | str = "scripts", ) -> Plan: """Return the full plan for one reference-case run. @@ -117,6 +159,13 @@ def build_plan( ``atmos_res`` is the grid the atmosphere and land were run on; every other realm's grid is derived from it or from the model. + ``workers`` applies to time series generation only; CMORization always runs + with one worker (see CMOR_WORKERS). + + ``variant_label`` and ``tables_root`` are passed to CMORization unchanged, + so the same command can produce output for submission rather than only + exercising the chain. + ``years`` is passed through to gen_timeseries.py as ``--years-spec`` and so uses its format, ``first:last:increment``. """ @@ -155,13 +204,28 @@ def build_plan( plan.skipped.append(f"{realm}: no history directory at {inputdir}") continue - realm_frequencies = _realm_frequencies(model, realm, frequencies) - if not realm_frequencies: + declared = _realm_frequencies(model, realm, frequencies) + if not declared: plan.skipped.append( f"{realm}: none of the requested frequencies are declared" ) continue + sheet_for_patterns = ( + ice_sheet or DEFAULT_ICE_SHEET if realm == "landIce" else None + ) + realm_frequencies = written_frequencies( + case_dir, model, realm, declared, sheet_for_patterns + ) + absent = [f for f in declared if f not in realm_frequencies] + if absent: + plan.skipped.append( + f"{realm}: no history files for {absent}; this case did not " + "write those streams" + ) + if not realm_frequencies: + continue + sheet = ice_sheet or DEFAULT_ICE_SHEET if realm == "landIce" else None ts_dir = timeseries_root / realm @@ -177,6 +241,7 @@ def build_plan( years=years, workers=workers, sheet=sheet, + overwrite=overwrite_timeseries, ) ) @@ -192,8 +257,9 @@ def build_plan( model=model, atmos_res=(atmos_res if needs_atmos_res(realm) else None), experiment=experiment, - workers=workers, sheet=sheet, + variant_label=variant_label, + tables_root=tables_root, ) ) if "validate" in stages: @@ -205,6 +271,7 @@ def build_plan( cmor_root=cmor_root, model=model, experiment=experiment, + variant_label=variant_label, ) ) @@ -222,11 +289,21 @@ def _realm_frequencies( def _timeseries_step( - *, scripts_dir, realm, inputdir, ts_dir, model, frequencies, years, workers, sheet + *, + scripts_dir, + realm, + inputdir, + ts_dir, + model, + frequencies, + years, + workers, + sheet, + overwrite=False, ) -> Step: """Return the gen_timeseries.py step for one realm.""" command = [ - "python", + PYTHON, str(Path(scripts_dir) / "gen_timeseries.py"), "--inputdir", str(inputdir), @@ -245,6 +322,8 @@ def _timeseries_step( command += ["--years-spec", years] if sheet: command += ["--ice-sheet", sheet] + if overwrite: + command.append("--overwrite_timeseries") return Step( key=f"timeseries-{realm}", stage="timeseries", realm=realm, command=command ) @@ -260,12 +339,13 @@ def _cmor_step( model, atmos_res, experiment, - workers, sheet, + variant_label=None, + tables_root=None, ) -> Step: """Return the cmor_driver.py step for one realm and frequency.""" command = [ - "python", + PYTHON, str(Path(scripts_dir) / "cmor_driver.py"), "--realm", realm, @@ -280,12 +360,16 @@ def _cmor_step( "--experiment", experiment, "--workers", - str(workers), + str(CMOR_WORKERS), ] if atmos_res: command += ["--atmos-res", atmos_res] if sheet: command += ["--ice-sheet", sheet] + if variant_label: + command += ["--variant-label", variant_label] + if tables_root: + command += ["--tables-root", str(tables_root)] return Step( key=f"cmor-{realm}-{frequency}", stage="cmor", @@ -296,11 +380,18 @@ def _cmor_step( def _validate_step( - *, scripts_dir, realm, frequency, cmor_root, model, experiment + *, + scripts_dir, + realm, + frequency, + cmor_root, + model, + experiment, + variant_label=None, ) -> Step: """Return the validate_cmor_output.py step for one realm and frequency.""" command = [ - "python", + PYTHON, str(Path(scripts_dir) / "validate_cmor_output.py"), "--model", model, @@ -314,6 +405,8 @@ def _validate_step( str(cmor_root), "--strict", ] + if variant_label: + command += ["--ensemble-member", variant_label] return Step( key=f"validate-{realm}-{frequency}", stage="validate", diff --git a/tests/test_reference_run.py b/tests/test_reference_run.py index 6de6a5f..c8720b1 100644 --- a/tests/test_reference_run.py +++ b/tests/test_reference_run.py @@ -1,7 +1,11 @@ """Tests for planning a full-chain reference-case run.""" +import sys +from pathlib import Path + import pytest +from cmip7_prep.include_patterns import all_include_patterns from cmip7_prep.reference_run import ( REALM_COMPONENT, STAGES, @@ -12,11 +16,42 @@ ) +def _write_history_files(case_dir, model, realms=None, frequencies=None): + """Create empty history files for the streams a case would have written. + + Only the names matter here: the plan is built by matching include patterns + against file names, never by opening them. + """ + for realm in realms or realms_for(model): + component = REALM_COMPONENT[realm] + inputdir = Path(case_dir) / component / "hist" + inputdir.mkdir(parents=True, exist_ok=True) + sheet = "gris" if realm == "landIce" else None + for frequency in frequencies_for(model, realm): + if frequencies and frequency not in frequencies: + continue + for pattern in all_include_patterns(model, realm, sheet, None, [frequency]): + name = f"case.{pattern.strip('.')}.1526-01.nc" + (inputdir / name).touch() + + +def _write_history_files_for_ice_sheet(case_dir, ice_sheet): + """Create land-ice history files for one ice sheet.""" + inputdir = Path(case_dir) / REALM_COMPONENT["landIce"] / "hist" + inputdir.mkdir(parents=True, exist_ok=True) + for frequency in frequencies_for("noresm", "landIce"): + for pattern in all_include_patterns( + "noresm", "landIce", ice_sheet, None, [frequency] + ): + (inputdir / f"case.{pattern.strip('.')}.1526-01.nc").touch() + + @pytest.fixture(name="case_dir") def case_dir_fixture(tmp_path): - """Return a case directory with every component's history directory.""" + """Return a case that wrote every stream its tables declare.""" for component in sorted(set(REALM_COMPONENT.values())): - (tmp_path / component / "hist").mkdir(parents=True) + (tmp_path / component / "hist").mkdir(parents=True, exist_ok=True) + _write_history_files(tmp_path, "noresm") return tmp_path @@ -134,12 +169,39 @@ def test_missing_component_is_skipped_not_fatal(self, tmp_path): is reported rather than ending the run. """ case = tmp_path / "case" - (case / "ice" / "hist").mkdir(parents=True) + _write_history_files(case, "noresm", realms=["seaIce"]) plan = _plan(case, tmp_path / "out") assert {step.realm for step in plan.steps} == {"seaIce"} assert len(plan.skipped) == 5 assert all("no history directory" in reason for reason in plan.skipped) + def test_only_frequencies_the_case_wrote_are_planned(self, tmp_path): + """A case that wrote only monthly output gets only monthly steps. + + A table lists every frequency a realm can produce, but a run writes + only the streams it was configured for. Planning the rest guarantees + a CMORization failure saying nothing about the code -- which is what + happened on a real ne30 case with no hourly output. + """ + case = tmp_path / "case" + _write_history_files(case, "noresm", realms=["atmos"], frequencies=["mon"]) + plan = _plan(case, tmp_path / "out", realms=["atmos"]) + assert [step.frequency for step in plan.for_stage("cmor")] == ["mon"] + assert any("no history files" in reason for reason in plan.skipped) + + def test_the_absent_frequencies_are_named(self, tmp_path): + """The report says which streams the case did not write.""" + case = tmp_path / "case" + _write_history_files( + case, "noresm", realms=["atmos"], frequencies=["mon", "day"] + ) + plan = _plan(case, tmp_path / "out", realms=["atmos"]) + reason = next(r for r in plan.skipped if "no history files" in r) + for absent in ("1hr", "3hr", "6hr"): + assert absent in reason + for present in ("'mon'", "'day'"): + assert present not in reason + def test_all_stages_present_by_default(self, case_dir, tmp_path): """The default plan runs all three stages.""" plan = _plan(case_dir, tmp_path / "out", realms=["seaIce"]) @@ -188,7 +250,13 @@ def test_ice_sheet_only_for_landice(self, case_dir, tmp_path): assert "--ice-sheet" not in _command_of(plan, "timeseries-seaIce") def test_ice_sheet_can_be_chosen(self, case_dir, tmp_path): - """The ice sheet given is the one passed on.""" + """The ice sheet given is the one passed on. + + The Antarctic files must exist for the plan to include them: the + patterns carry the ice sheet name, so a Greenland-only case offers no + Antarctic streams to process. + """ + _write_history_files_for_ice_sheet(case_dir, "ais") plan = _plan(case_dir, tmp_path / "out", realms=["landIce"], ice_sheet="ais") command = _command_of(plan, "timeseries-landIce") assert command[command.index("--ice-sheet") + 1] == "ais" @@ -236,6 +304,80 @@ def test_resolution_and_experiment_reach_cmor(self, case_dir, tmp_path): assert command[command.index("--atmos-res") + 1] == "ne30" assert command[command.index("--experiment") + 1] == "historical" + def test_variant_label_reaches_both_later_stages(self, case_dir, tmp_path): + """The ensemble member labels the output and the validation of it. + + Without it every member would be written to the same path, so a real + cmorisation run needs it passed through rather than defaulted. + """ + plan = _plan( + case_dir, tmp_path / "out", realms=["seaIce"], variant_label="r2i1p1f1" + ) + cmor = _command_of(plan, "cmor-seaIce-mon") + validate = _command_of(plan, "validate-seaIce-mon") + assert cmor[cmor.index("--variant-label") + 1] == "r2i1p1f1" + assert validate[validate.index("--ensemble-member") + 1] == "r2i1p1f1" + + def test_variant_label_is_absent_when_not_given(self, case_dir, tmp_path): + """Without one, the flag is omitted so each script's default applies.""" + plan = _plan(case_dir, tmp_path / "out", realms=["seaIce"]) + assert "--variant-label" not in _command_of(plan, "cmor-seaIce-mon") + assert "--ensemble-member" not in _command_of(plan, "validate-seaIce-mon") + + def test_tables_root_reaches_cmorization(self, case_dir, tmp_path): + """A chosen tables checkout is used for the writing, not just the check.""" + plan = _plan( + case_dir, + tmp_path / "out", + realms=["seaIce"], + tables_root="/some/cmip7-cmor-tables", + ) + command = _command_of(plan, "cmor-seaIce-mon") + assert command[command.index("--tables-root") + 1] == "/some/cmip7-cmor-tables" + + def test_cmorization_always_gets_one_worker(self, case_dir, tmp_path): + """Time series generation is parallel; CMORization is not. + + The two --workers mean different things -- processes reading files + versus Dask workers -- so a single number cannot serve both. + """ + plan = _plan(case_dir, tmp_path / "out", realms=["seaIce"], workers=16) + timeseries = _command_of(plan, "timeseries-seaIce") + cmor = _command_of(plan, "cmor-seaIce-mon") + assert timeseries[timeseries.index("--workers") + 1] == "16" + assert cmor[cmor.index("--workers") + 1] == "1" + + def test_timeseries_resumes_unless_overwriting_was_asked_for( + self, case_dir, tmp_path + ): + """Overwriting is requested explicitly, so a repeat run resumes. + + gen_timeseries.py keeps existing output and produces only the missing + timesteps unless told otherwise, which is what makes repeating an + interrupted run cheap. + """ + resuming = _plan(case_dir, tmp_path / "out", realms=["seaIce"]) + rewriting = _plan( + case_dir, + tmp_path / "out2", + realms=["seaIce"], + overwrite_timeseries=True, + ) + assert "--overwrite_timeseries" not in _command_of( + resuming, "timeseries-seaIce" + ) + assert "--overwrite_timeseries" in _command_of(rewriting, "timeseries-seaIce") + + def test_steps_run_with_the_interpreter_running_the_plan(self, case_dir, tmp_path): + """Every step uses this interpreter, not whatever 'python' means. + + The stages have to run in the environment the run was started from, + and 'python' may not be on PATH at all. + """ + plan = _plan(case_dir, tmp_path / "out", realms=["seaIce"]) + for step in plan.steps: + assert step.command[0] == sys.executable + def test_every_step_has_a_distinct_log_name(self, case_dir, tmp_path): """Logs cannot overwrite one another.""" plan = _plan(case_dir, tmp_path / "out") From 103d04ea20d6035ad7f8e7a087c83e49eb366bed Mon Sep 17 00:00:00 2001 From: mvertens Date: Mon, 5 Oct 2026 14:13:02 +0200 Subject: [PATCH 08/16] fixes to show the real run order --- scripts/run_reference_case.py | 23 ++++++++++++++--------- src/cmip7_prep/reference_run.py | 16 ++++++++++++++-- tests/test_reference_run.py | 22 +++++++++++++++++++--- 3 files changed, 47 insertions(+), 14 deletions(-) diff --git a/scripts/run_reference_case.py b/scripts/run_reference_case.py index c39138b..e78be1f 100644 --- a/scripts/run_reference_case.py +++ b/scripts/run_reference_case.py @@ -247,15 +247,20 @@ def run_step(step: Step, log_dir: Path) -> dict: def print_plan(plan: Plan) -> None: - """Print the commands a run would issue, without running them.""" - for stage in STAGES: - steps = plan.for_stage(stage) - if not steps: - continue - print(f"\n=== {stage} ({len(steps)} step(s)) ===") - for step in steps: - print(f" {step.key}") - print(f" {' '.join(step.command)}") + """Print the commands a run would issue, in the order it would issue them. + + Grouped by realm, because that is how a run proceeds: one realm's time + series, CMORization and validation complete before the next realm starts. + """ + realm = None + for step in plan.steps: + if step.realm != realm: + realm = step.realm + count = sum(1 for s in plan.steps if s.realm == realm) + print(f"\n=== {realm} ({count} step(s)) ===") + print(f" {step.key}") + print(f" {' '.join(step.command)}") + print(f"\n{len(plan.steps)} step(s) in total") if plan.skipped: print("\n=== skipped ===") for reason in plan.skipped: diff --git a/src/cmip7_prep/reference_run.py b/src/cmip7_prep/reference_run.py index b839c03..a25b0f1 100644 --- a/src/cmip7_prep/reference_run.py +++ b/src/cmip7_prep/reference_run.py @@ -44,6 +44,11 @@ STAGES = ("timeseries", "cmor", "validate") +# Frequencies are processed coarsest first, so the cheapest and most used +# output exists soonest. Alphabetical order would put 'day' before 'mon' and +# '1hr' first of all, leaving the monthly output until last. +FREQUENCY_ORDER = ("yr", "mon", "day", "6hr", "3hr", "1hr") + # The steps run with the interpreter running this, not whatever "python" means # on PATH: that keeps every stage in the environment the run was started from, # and works when the runner was invoked as python3 or by absolute path. @@ -96,8 +101,15 @@ def realms_for(model: str) -> list[str]: def frequencies_for(model: str, realm: str) -> list[str]: - """Return the frequencies declared for one model and realm.""" - return sorted(load_include_patterns(model)[realm]) + """Return the frequencies declared for one model and realm, coarsest first. + + The order is the order a run will process them in, so it is chosen rather + than alphabetical: see FREQUENCY_ORDER. + """ + declared = load_include_patterns(model)[realm] + known = [f for f in FREQUENCY_ORDER if f in declared] + rest = sorted(f for f in declared if f not in FREQUENCY_ORDER) + return known + rest def written_frequencies( diff --git a/tests/test_reference_run.py b/tests/test_reference_run.py index c8720b1..285fe55 100644 --- a/tests/test_reference_run.py +++ b/tests/test_reference_run.py @@ -124,12 +124,28 @@ def test_cmor_and_validate_are_per_frequency(self, case_dir, tmp_path): """CMOR and validation run once per realm and frequency.""" plan = _plan(case_dir, tmp_path / "out", realms=["seaIce"]) assert [step.key for step in plan.for_stage("cmor")] == [ - "cmor-seaIce-day", "cmor-seaIce-mon", + "cmor-seaIce-day", ] assert [step.key for step in plan.for_stage("validate")] == [ - "validate-seaIce-day", "validate-seaIce-mon", + "validate-seaIce-day", + ] + + def test_coarse_frequencies_are_processed_first(self, case_dir, tmp_path): + """Monthly output is produced before daily, and daily before hourly. + + A run that stops partway should have produced the most useful output, + and CMORizing daily data took 24 minutes on a real two-year case while + monthly took a fraction of that. + """ + plan = _plan(case_dir, tmp_path / "out", realms=["atmos"]) + assert [step.frequency for step in plan.for_stage("cmor")] == [ + "mon", + "day", + "6hr", + "3hr", + "1hr", ] def test_realms_can_be_narrowed(self, case_dir, tmp_path): @@ -142,7 +158,7 @@ def test_frequencies_can_be_narrowed(self, case_dir, tmp_path): plan = _plan( case_dir, tmp_path / "out", realms=["atmos"], frequencies=["mon", "day"] ) - assert [step.frequency for step in plan.for_stage("cmor")] == ["day", "mon"] + assert [step.frequency for step in plan.for_stage("cmor")] == ["mon", "day"] def test_frequency_a_realm_does_not_declare_is_skipped(self, case_dir, tmp_path): """Asking seaIce for 1hr skips it rather than planning a doomed step.""" From 91aee3a46d19e10a1def6d2eca6bc0c2b071d498 Mon Sep 17 00:00:00 2001 From: mvertens Date: Tue, 6 Oct 2026 09:57:05 +0200 Subject: [PATCH 09/16] frequencies now run monthly first, then daily, then hourly --- src/cmip7_prep/reference_run.py | 21 ++++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/src/cmip7_prep/reference_run.py b/src/cmip7_prep/reference_run.py index a25b0f1..781fd26 100644 --- a/src/cmip7_prep/reference_run.py +++ b/src/cmip7_prep/reference_run.py @@ -160,6 +160,9 @@ def build_plan( overwrite_timeseries: bool = False, variant_label: str | None = None, tables_root: os.PathLike | str | None = None, + plots: bool = False, + html: bool = False, + max_plots: int | None = None, scripts_dir: os.PathLike | str = "scripts", ) -> Plan: """Return the full plan for one reference-case run. @@ -176,7 +179,9 @@ def build_plan( ``variant_label`` and ``tables_root`` are passed to CMORization unchanged, so the same command can produce output for submission rather than only - exercising the chain. + exercising the chain. ``plots`` and ``html`` ask validation for figures + and a browsable index of them, which is how a person judges whether the + output is sensible rather than merely complete. ``years`` is passed through to gen_timeseries.py as ``--years-spec`` and so uses its format, ``first:last:increment``. @@ -284,6 +289,9 @@ def build_plan( model=model, experiment=experiment, variant_label=variant_label, + plots=plots, + html=html, + max_plots=max_plots, ) ) @@ -400,6 +408,9 @@ def _validate_step( model, experiment, variant_label=None, + plots=False, + html=False, + max_plots=None, ) -> Step: """Return the validate_cmor_output.py step for one realm and frequency.""" command = [ @@ -419,6 +430,14 @@ def _validate_step( ] if variant_label: command += ["--ensemble-member", variant_label] + if plots: + # Both kinds: a time series per variable says whether the values are + # plausible over time, a map says whether they are in the right place. + command += ["--plot-timeseries", "--plot-maps"] + if max_plots is not None: + command += ["--max-plots", str(max_plots)] + if html: + command.append("--html") return Step( key=f"validate-{realm}-{frequency}", stage="validate", From 12436be72bdcd17ffa8c7d963d55f5c932ff699f Mon Sep 17 00:00:00 2001 From: mvertens Date: Tue, 6 Oct 2026 09:57:31 +0200 Subject: [PATCH 10/16] --dry-run now lists steps in the order they will actually run --- scripts/cmor_driver.py | 1 + scripts/run_reference_case.py | 28 +++++++++ src/cmip7_prep/cmor_writer.py | 114 ++++++++++++++++++++++++++++------ tests/test_cmor_writer.py | 77 +++++++++++++++++++++++ tests/test_reference_run.py | 28 +++++++++ 5 files changed, 230 insertions(+), 18 deletions(-) diff --git a/scripts/cmor_driver.py b/scripts/cmor_driver.py index 52cd249..5506b5e 100755 --- a/scripts/cmor_driver.py +++ b/scripts/cmor_driver.py @@ -638,6 +638,7 @@ def process_one_var( dataset_attrs={"institution_id": "NCC", "GLOBAL_IS_CMIP7": True}, outdir=outdir, ice_sheet=ice_sheet, + realm=realm, ) as cm: set_cur_dataset_attribute("frequency", frequency) set_cur_dataset_attribute("realization_index", realization_index) diff --git a/scripts/run_reference_case.py b/scripts/run_reference_case.py index e78be1f..eb27e90 100644 --- a/scripts/run_reference_case.py +++ b/scripts/run_reference_case.py @@ -140,6 +140,31 @@ def parse_arguments(): help="Ice sheet for the landIce realm (default: gris)", ) + reporting = parser.add_argument_group("what to report") + reporting.add_argument( + "--plots", + action="store_true", + help=( + "Have validation plot each variable: a mean time series, and a " + "time-mean map where the data are on a plottable grid" + ), + ) + reporting.add_argument( + "--html", + action="store_true", + help=( + "Have validation rebuild its static HTML index, so the reports " + "and plots for every realm can be browsed in one place" + ), + ) + reporting.add_argument( + "--max-plots", + type=int, + default=None, + metavar="N", + help="Plot at most N variables per plot type (default: no limit)", + ) + behaviour = parser.add_argument_group("how to run") behaviour.add_argument( "--workers", @@ -331,6 +356,9 @@ def main(): workers=args.workers, ice_sheet=args.ice_sheet, overwrite_timeseries=args.force, + plots=args.plots, + html=args.html, + max_plots=args.max_plots, variant_label=args.variant_label, tables_root=args.tables_root, scripts_dir=_LOCAL_PATH, diff --git a/src/cmip7_prep/cmor_writer.py b/src/cmip7_prep/cmor_writer.py index 549f1cc..002dd0a 100644 --- a/src/cmip7_prep/cmor_writer.py +++ b/src/cmip7_prep/cmor_writer.py @@ -68,6 +68,60 @@ def _horizontal_only(da_coord: xr.DataArray, keep) -> xr.DataArray: # --------------------------------------------------------------------- # CMOR session # --------------------------------------------------------------------- + +# One unit of time is written differently in different places: a year is +# "common_year" in a CISM units string but "year" in a time_period_freq +# attribute. Each variant maps to a single name below, so that the unit a +# period is given in can be compared with the unit an axis counts in. +_TIME_UNIT_NAMES = { + "common_year": "year", + "common_years": "year", + "year": "year", + "years": "year", + "month": "month", + "months": "month", + "day": "day", + "days": "day", + "hour": "hour", + "hours": "hour", +} + + +def _time_unit_name(word: str) -> str | None: + """Return the one name for a unit of time, or None if it is unrecognised.""" + return _TIME_UNIT_NAMES.get(str(word).strip().lower()) + + +def period_in_axis_units(dataset, units: str) -> float | None: + """Return the averaging period, in the unit the time axis counts in. + + The model records its averaging period in a ``time_period_freq`` + attribute, written as a unit and a count: ``"year_1"`` for an annual mean, + ``"month_1"`` for a monthly one. The time axis counts in a unit of its + own, named at the start of its units string -- ``"common_year since + 0000-01-01"`` counts in years. + + The period is returned only when those two units are the same, since + converting between them is not always possible: a year is 360, 365 or + 365.25 days depending on the calendar, and a month has no fixed length at + all. Anything else returns None, leaving the caller to say so rather than + record a period the data does not have. + """ + period = getattr(dataset, "attrs", {}).get("time_period_freq") + if not period or not units: + return None + + name, separator, count = str(period).strip().rpartition("_") + if not separator or not count.isdigit(): + return None + + period_unit = _time_unit_name(name) + axis_unit = _time_unit_name(str(units).split()[0]) + if period_unit is None or period_unit != axis_unit: + return None + return float(count) + + class CmorSession( AbstractContextManager ): # pylint: disable=too-many-instance-attributes @@ -96,6 +150,7 @@ def __init__( log_name: str | None = None, outdir: Path | str | None = None, ice_sheet: str | None = None, + realm: str | None = None, ) -> None: self.tables_root = tables_root self.dataset_attrs = dict(dataset_attrs or {}) @@ -105,6 +160,7 @@ def __init__( # Ice sheet ('gris'/'ais') selecting the CISM map projection; only used # for the native land-ice grid (see _define_cism_grid). self._ice_sheet = ice_sheet + self._realm = realm # logging config self._log_dir = Path(log_dir) if log_dir is not None else None self._log_name = log_name @@ -625,25 +681,47 @@ def _get_time_and_bounds(dsi: xr.Dataset, want_bounds: bool = True): return tvals, None, str(units) tbnum = encode_time_to_num(tb, units, cal) if tb is not None else None - # Some model output (e.g. CISM land-ice) carries no time bounds, but - # CMOR requires them for time-averaged variables. When absent, - # synthesize bounds from the numeric time centers using the midpoints - # between consecutive steps (edges extrapolated). For evenly-spaced - # data this reproduces the averaging period (e.g. the calendar year - # for annual means). - if tbnum is None and tvals is not None and np.size(tvals) >= 2: - t = np.asarray(tvals, dtype="f8").reshape(-1) - mids = 0.5 * (t[:-1] + t[1:]) - tbnum = np.empty((t.size, 2), dtype="f8") - tbnum[1:, 0] = mids - tbnum[:-1, 1] = mids - tbnum[0, 0] = t[0] - (mids[0] - t[0]) - tbnum[-1, 1] = t[-1] + (t[-1] - mids[-1]) - logger.info( - "time coordinate had no bounds; synthesized %d bounds from " - "centers", - tbnum.shape[0], + # If CISM land-ice output carries no time bounds, it is built + # from the time values: each period is taken to END at its own time + # value, which is how the output is labelled -- a file stamped 1527 + # holds the mean for 1526. + # + # Only land ice is treated this way. Every other component writes + # bounds, so their absence means something went wrong earlier and is + # reported rather than papered over. + if tbnum is None and tvals is not None and self._realm != "landIce": + raise ValueError( + f"Time coordinate has no bounds and realm is " + f"{self._realm!r}, which is expected to write them. " + "Check the history files and the time series built from " + "them; bounds are only synthesized for landIce." ) + if tbnum is None and tvals is not None and np.size(tvals) >= 1: + t = np.asarray(tvals, dtype="f8").reshape(-1) + left = None + if t.size >= 2: + # Each period starts where the previous one ended; the + # first is extrapolated back by its own length. + left = np.empty_like(t) + left[1:] = t[:-1] + left[0] = t[0] - (t[1] - t[0]) + else: + period = period_in_axis_units(dsi, units) + if period is None: + logger.warning( + "time coordinate has no bounds, a single step, and " + "no usable time_period_freq attribute; CMOR will " + "reject the axis if the table requires bounds" + ) + else: + left = np.array([t[0] - period], dtype="f8") + if left is not None: + tbnum = np.column_stack([left, t]) + logger.info( + "time coordinate had no bounds; synthesized %d, each " + "ending at its time value", + tbnum.shape[0], + ) return tvals, tbnum, str(units) diff --git a/tests/test_cmor_writer.py b/tests/test_cmor_writer.py index 3371100..f83c9e2 100644 --- a/tests/test_cmor_writer.py +++ b/tests/test_cmor_writer.py @@ -3,6 +3,7 @@ import json from pathlib import Path import numpy as np +import pytest import xarray as xr from cmip7_prep.cmor_writer import CmorSession @@ -232,3 +233,79 @@ class BrandedName: log_name="cmor_test_zonal_mean.log", ) as session: session.write_variable(ds, CMIPVar(), VDef()) + + +def _dataset_without_time_bounds(): + """Return a dataset whose time coordinate carries no bounds.""" + data = np.random.rand(2, 4, 8) + ds = xr.Dataset( + { + "tas_tmin-h2m-hxy-u": (("time", "lat", "lon"), data), + "lat": ("lat", np.linspace(-90, 90, 4)), + "lon": ("lon", np.linspace(0, 360, 8, endpoint=False)), + "time": ("time", np.array([15, 45])), + } + ) + ds["lat"].attrs["units"] = "degrees_north" + ds["lon"].attrs["units"] = "degrees_east" + ds["time"].attrs["units"] = "days since 2000-01-01" + ds["time"].attrs["calendar"] = "noleap" + ds.attrs["branded_variable"] = "tas_tmin-h2m-hxy-u" + return ds + + +class _VDef: # pylint: disable=too-few-public-methods + """Minimal variable definition.""" + + name = "tas" + branded_variable_name = "tas_tmin-h2m-hxy-u" + units = "K" + table = "atmos" + levels: dict = {} + + +class _CMIPVar: # pylint: disable=too-few-public-methods + """Minimal CMIP variable wrapper.""" + + class BrandedName: # pylint: disable=too-few-public-methods + """Branded variable name wrapper.""" + + name = "tas_tmin-h2m-hxy-u" + + branded_variable_name = BrandedName() + + +def _session(tmp_path, realm): + """Return a CmorSession for one realm, writing into tmp_path.""" + return CmorSession( + tables_root=Path(__file__).parent.parent / "cmip7-cmor-tables", + dataset_json=Path(__file__).parent.parent / "data" / "cmor_dataset.json", + dataset_attrs={"institution_id": "NCAR", "GLOBAL_IS_CMIP7": True}, + log_dir=tmp_path, + log_name=f"cmor_{realm}.log", + realm=realm, + ) + + +def test_missing_time_bounds_are_an_error_outside_land_ice(tmp_path): + """A realm expected to write time bounds fails loudly without them. + + Only CISM land-ice output lacks bounds, so their absence anywhere else + means something went wrong before CMORization and should be reported, not + replaced with a plausible guess. + """ + with _session(tmp_path, "atmos") as session: + with pytest.raises(ValueError, match="expected to write them"): + session.write_variable(_dataset_without_time_bounds(), _CMIPVar(), _VDef()) + + +def test_missing_time_bounds_are_synthesized_for_land_ice(tmp_path): + """Land ice is the exception: its bounds are built from the time values.""" + with _session(tmp_path, "landIce") as session: + try: + session.write_variable(_dataset_without_time_bounds(), _CMIPVar(), _VDef()) + except Exception as exc: # pylint: disable=broad-except + # Getting as far as CMOR's own checks means the bounds were built. + # What CMOR then makes of an atmosphere variable in the land-ice + # table is not what this test is about. + assert "expected to write them" not in str(exc) diff --git a/tests/test_reference_run.py b/tests/test_reference_run.py index 285fe55..0a0365b 100644 --- a/tests/test_reference_run.py +++ b/tests/test_reference_run.py @@ -394,6 +394,34 @@ def test_steps_run_with_the_interpreter_running_the_plan(self, case_dir, tmp_pat for step in plan.steps: assert step.command[0] == sys.executable + def test_plots_are_not_produced_unless_asked_for(self, case_dir, tmp_path): + """Plotting every variable is expensive, so it is opt-in.""" + plan = _plan(case_dir, tmp_path / "out", realms=["seaIce"]) + command = _command_of(plan, "validate-seaIce-mon") + for flag in ("--plot-timeseries", "--plot-maps", "--html", "--max-plots"): + assert flag not in command + + def test_plots_are_requested_when_asked_for(self, case_dir, tmp_path): + """Both kinds of plot are asked for together.""" + plan = _plan(case_dir, tmp_path / "out", realms=["seaIce"], plots=True) + command = _command_of(plan, "validate-seaIce-mon") + assert "--plot-timeseries" in command + assert "--plot-maps" in command + + def test_html_and_plot_limit_are_forwarded(self, case_dir, tmp_path): + """The browsable index and the plot limit reach validation.""" + plan = _plan( + case_dir, + tmp_path / "out", + realms=["seaIce"], + plots=True, + html=True, + max_plots=12, + ) + command = _command_of(plan, "validate-seaIce-mon") + assert "--html" in command + assert command[command.index("--max-plots") + 1] == "12" + def test_every_step_has_a_distinct_log_name(self, case_dir, tmp_path): """Logs cannot overwrite one another.""" plan = _plan(case_dir, tmp_path / "out") From 0292f3385b8c179140c3496b72e336453686823e Mon Sep 17 00:00:00 2001 From: mvertens Date: Tue, 6 Oct 2026 10:01:09 +0200 Subject: [PATCH 11/16] new test for time bounds --- tests/test_time_period_bounds.py | 126 +++++++++++++++++++++++++++++++ 1 file changed, 126 insertions(+) create mode 100644 tests/test_time_period_bounds.py diff --git a/tests/test_time_period_bounds.py b/tests/test_time_period_bounds.py new file mode 100644 index 0000000..d821c4f --- /dev/null +++ b/tests/test_time_period_bounds.py @@ -0,0 +1,126 @@ +"""Tests for deriving an averaging period from a file's own attributes. + +Model output without time bounds has to have them synthesized, since CMOR +refuses a time axis without bounds when the table expects a mean. With two or +more time steps the spacing gives the period; with one step it can only come +from the file, via the ``time_period_freq`` attribute the output +carries. +""" + +import numpy as np +import xarray as xr + +from cmip7_prep.cmor_writer import period_in_axis_units + +YEAR_AXIS = "common_year since 0000-01-01 0:0:0" +DAY_AXIS = "days since 1850-01-01" + + +def _dataset(period=None): + """Return a dataset carrying one time value and an optional period.""" + attrs = {"time_period_freq": period} if period else {} + return xr.Dataset( + {"acab": ("time", np.array([1.0]))}, + coords={"time": ("time", np.array([1527.0]))}, + attrs=attrs, + ) + + +def _bounds_from(times, period=None): + """Return the bounds the writer builds for a set of time values. + + Mirrors the arithmetic in cmor_writer: each period ends at its own time + value, because that is how the output is labelled. + """ + times = np.asarray(times, dtype="f8") + if times.size >= 2: + left = np.empty_like(times) + left[1:] = times[:-1] + left[0] = times[0] - (times[1] - times[0]) + else: + left = np.array([times[0] - period], dtype="f8") + return np.column_stack([left, times]).tolist() + + +class TestEndStampedBounds: + """Each period ends at its time value, not centred on it. + + A CISM file stamped 1527 holds the annual mean for 1526, so its bounds run + 1526 to 1527. Centring would give 1526.5 to 1527.5: the right length, the + wrong year. + """ + + def test_one_annual_step(self): + """One step uses the period the file records.""" + assert _bounds_from([1527.0], period=1.0) == [[1526.0, 1527.0]] + + def test_two_annual_steps(self): + """Consecutive periods meet at the time values, on calendar years.""" + assert _bounds_from([1527.0, 1528.0]) == [ + [1526.0, 1527.0], + [1527.0, 1528.0], + ] + + def test_periods_are_contiguous(self): + """No gaps or overlaps between consecutive periods.""" + bounds = _bounds_from([1527.0, 1528.0, 1529.0]) + for earlier, later in zip(bounds, bounds[1:]): + assert earlier[1] == later[0] + + def test_one_step_and_many_steps_agree(self): + """The first period is the same whichever path produced it.""" + alone = _bounds_from([1527.0], period=1.0)[0] + together = _bounds_from([1527.0, 1528.0])[0] + assert alone == together + + +class TestMatchingUnits: + """Periods expressed in the unit the axis already counts in.""" + + def test_annual_period_on_a_year_axis(self): + """The land-ice case: one annual step on an axis counted in years.""" + assert period_in_axis_units(_dataset("year_1"), YEAR_AXIS) == 1.0 + + def test_multi_year_period(self): + """A five-year mean is five years long.""" + assert period_in_axis_units(_dataset("year_5"), YEAR_AXIS) == 5.0 + + def test_daily_period_on_a_day_axis(self): + """A daily mean is one day long on an axis counted in days.""" + assert period_in_axis_units(_dataset("day_1"), DAY_AXIS) == 1.0 + + +class TestRefusals: + """Cases with no exact answer, where guessing would misstate the period.""" + + def test_years_against_a_day_axis(self): + """A year is 360, 365 or 365.25 days depending on the calendar.""" + assert period_in_axis_units(_dataset("year_1"), DAY_AXIS) is None + + def test_months_against_a_day_axis(self): + """A month has no fixed length at all.""" + assert period_in_axis_units(_dataset("month_1"), DAY_AXIS) is None + + def test_hours_against_a_day_axis(self): + """Exact conversions are not attempted either, only matching units. + + Six hours is a quarter of a day in every calendar, but no output this + is used on needs that, and converting is where mistakes would hide. + """ + assert period_in_axis_units(_dataset("hour_6"), DAY_AXIS) is None + + def test_no_period_attribute(self): + """Output that does not record its period yields nothing.""" + assert period_in_axis_units(_dataset(), YEAR_AXIS) is None + + def test_unparseable_period(self): + """An attribute that is not _ is not guessed at.""" + assert period_in_axis_units(_dataset("annual"), YEAR_AXIS) is None + + def test_unknown_unit(self): + """A unit not in the table is refused rather than assumed.""" + assert period_in_axis_units(_dataset("fortnight_1"), YEAR_AXIS) is None + + def test_empty_units_string(self): + """An axis with no units cannot be converted to.""" + assert period_in_axis_units(_dataset("year_1"), "") is None From 4c67d6a746bba6abf8ffac521425c58c3a3b4b13 Mon Sep 17 00:00:00 2001 From: mvertens Date: Tue, 6 Oct 2026 12:13:48 +0200 Subject: [PATCH 12/16] fix testing issues --- tests/test_cv_lookup.py | 41 +++++++++++++++++++++++++++++------------ 1 file changed, 29 insertions(+), 12 deletions(-) diff --git a/tests/test_cv_lookup.py b/tests/test_cv_lookup.py index f4fab3b..316e6bc 100644 --- a/tests/test_cv_lookup.py +++ b/tests/test_cv_lookup.py @@ -1,6 +1,7 @@ """Tests for reading the controlled vocabulary to reject bad values early.""" import json +from pathlib import Path import pytest @@ -50,15 +51,27 @@ def test_a_valid_experiment_passes(self): """A registered experiment validates silently.""" cv_lookup.validate(TABLES_ROOT, "experiment_id", "piControl") - def test_the_license_we_fixed_is_registered(self): - """CC-BY-4.0 is the registered spelling, not CC-BY-4-0. - - A typo here is what stopped CMOR writing anything at all, so it is - worth asserting the vocabulary still agrees with data/cmor_dataset.json. + @pytest.mark.parametrize( + "dataset_file", ["cmor_dataset.json", "cmor_dataset_noresm.json"] + ) + def test_the_license_we_request_is_registered(self, dataset_file): + """Whatever a dataset file asks for must be in the vocabulary. + + CMOR refuses to write anything at all when it is not, and the message + it gives names only the value, not the file it came from. The spelling + itself is not asserted: it has differed between tables versions + ("CC-BY-4.0" against "CC-BY-4-0"), so pinning one here would fail on a + tables update rather than on a real mismatch. """ - licenses = cv_lookup.allowed_values(TABLES_ROOT, "license.license_id") - assert "CC-BY-4.0" in licenses - assert "CC-BY-4-0" not in licenses + path = Path(__file__).parent.parent / "data" / dataset_file + if not path.is_file(): + pytest.skip(f"{dataset_file} is not part of this branch") + requested = json.loads(path.read_text(encoding="utf-8"))["license_id"] + registered = cv_lookup.allowed_values(TABLES_ROOT, "license.license_id") + assert requested in registered, ( + f"{dataset_file} asks for license_id {requested!r}, which the " + f"controlled vocabulary does not register; it allows {registered}" + ) # --------------------------------------------------------- the stand-in tables @@ -80,10 +93,14 @@ def test_a_list_valued_key_works(self, fake_tables): assert cv_lookup.allowed_values(fake_tables, "license") == ["CC-BY-4.0"] def test_a_dotted_path_walks_into_nested_keys(self): - """The licences nest under license.license_id in the real vocabulary.""" - assert "CC-BY-4.0" in cv_lookup.allowed_values( - TABLES_ROOT, "license.license_id" - ) + """The licences nest under license.license_id in the real vocabulary. + + What they are called is not this test's business -- only that the path + reaches them and returns something. + """ + licenses = cv_lookup.allowed_values(TABLES_ROOT, "license.license_id") + assert licenses + assert all(isinstance(name, str) and name for name in licenses) def test_a_dotted_path_that_goes_nowhere_raises(self): """A wrong path says what is controlled at the level it reached.""" From 1bab9a474d0c7709f14583c07160ab888de53068 Mon Sep 17 00:00:00 2001 From: mvertens Date: Tue, 6 Oct 2026 13:10:21 +0200 Subject: [PATCH 13/16] fixed more problem in CI testing --- .github/workflows/pytest.yaml | 7 ++----- data/cmor_dataset.json | 2 +- 2 files changed, 3 insertions(+), 6 deletions(-) diff --git a/.github/workflows/pytest.yaml b/.github/workflows/pytest.yaml index a91dc7e..7444e40 100644 --- a/.github/workflows/pytest.yaml +++ b/.github/workflows/pytest.yaml @@ -12,12 +12,9 @@ jobs: steps: - name: Checkout uses: actions/checkout@v4 - - name: Checkout cmip7-cmor-tables from CESM-Development fork and cesm-dev branch - uses: actions/checkout@v4 with: - repository: CESM-Development/cmip7-cmor-tables - ref: cesm-dev - path: cmip7-cmor-tables + # The CMOR tables come from the submodule. + submodules: recursive - name: Set up Python uses: actions/setup-python@v5 with: diff --git a/data/cmor_dataset.json b/data/cmor_dataset.json index fcd2666..449b1cc 100644 --- a/data/cmor_dataset.json +++ b/data/cmor_dataset.json @@ -31,7 +31,7 @@ "sub_experiment": "none", "archive_id": "WCRP", "area_label" : "u", - "license_id": "CC-BY-4.0", + "license_id": "CC-BY-4-0", "realization_index": "r1", "initialization_index": "i1", "physics_index": "p1", From e69d1d131dce6ec5aa661722f2f96cd5c9372658 Mon Sep 17 00:00:00 2001 From: mvertens Date: Tue, 6 Oct 2026 13:18:29 +0200 Subject: [PATCH 14/16] Trigger CI with the submodule-based tables From ea0b0d2bf4be5f0d37d4ac16cc0c905d61e3a287 Mon Sep 17 00:00:00 2001 From: mvertens Date: Tue, 6 Oct 2026 13:19:06 +0200 Subject: [PATCH 15/16] Show which CMOR tables CI is using --- .github/workflows/pytest.yaml | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/.github/workflows/pytest.yaml b/.github/workflows/pytest.yaml index 7444e40..d45e4b8 100644 --- a/.github/workflows/pytest.yaml +++ b/.github/workflows/pytest.yaml @@ -15,6 +15,12 @@ jobs: with: # The CMOR tables come from the submodule. submodules: recursive + - name: Show which CMOR tables are in use + run: | + git -C cmip7-cmor-tables log -1 --format='%H %s' \ + || echo "cmip7-cmor-tables is not a git checkout" + grep -o '"CC-BY-4[^"]*"' cmip7-cmor-tables/tables-cvs/cmor-cvs.json \ + | sort -u - name: Set up Python uses: actions/setup-python@v5 with: From d80c99e6b20273af607939027273856ab7f5a6a6 Mon Sep 17 00:00:00 2001 From: mvertens Date: Tue, 6 Oct 2026 13:36:42 +0200 Subject: [PATCH 16/16] more fixes to get testing to work --- .github/workflows/pytest.yaml | 6 +++--- data/cmor_dataset.json | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/.github/workflows/pytest.yaml b/.github/workflows/pytest.yaml index d45e4b8..d4aedaa 100644 --- a/.github/workflows/pytest.yaml +++ b/.github/workflows/pytest.yaml @@ -51,14 +51,14 @@ jobs: /home/runner/cmor /home/runner/cmor3-source /home/runner/.cache/pip - key: cmor-${{ runner.os }}-${{ hashFiles('requirements.txt') }}-${{ hashFiles('**/cmor_version.txt') }}-cmor3.14.0 + key: cmor-${{ runner.os }}-${{ hashFiles('requirements.txt') }}-${{ hashFiles('**/cmor_version.txt') }}-cmor3.15.3 restore-keys: | - cmor-${{ runner.os }}-3.14.0 + cmor-${{ runner.os }}-3.15.3 - name: build and install cmor run: | if [ ! -d "$HOME/cmor" ]; then - git clone https://github.com/PCMDI/cmor.git $HOME/cmor3-source -b 3.14.0 + git clone https://github.com/PCMDI/cmor.git $HOME/cmor3-source -b 3.15.3 cd $HOME/cmor3-source ./configure --prefix=$HOME/cmor --with-netcdf=/usr/local --with-udunits2=/usr/local make diff --git a/data/cmor_dataset.json b/data/cmor_dataset.json index 449b1cc..fcd2666 100644 --- a/data/cmor_dataset.json +++ b/data/cmor_dataset.json @@ -31,7 +31,7 @@ "sub_experiment": "none", "archive_id": "WCRP", "area_label" : "u", - "license_id": "CC-BY-4-0", + "license_id": "CC-BY-4.0", "realization_index": "r1", "initialization_index": "i1", "physics_index": "p1",