From a291dd60026922742ec9b6cd6a45b7a899aaeed1 Mon Sep 17 00:00:00 2001 From: Bouwe Andela Date: Tue, 21 Jan 2025 13:01:20 +0100 Subject: [PATCH] Add ESMValTool Transient Climate Response recipe (#62) --- changelog/62.feature.md | 1 + .../cmip_ref_metrics_esmvaltool/__init__.py | 9 +- .../metrics/__init__.py | 2 + .../metrics/tcr.py | 117 ++++++++++++++ .../cmip_ref_metrics_esmvaltool/recipes.txt | 1 + .../tests/unit/metrics/input_files_tcr.json | 146 ++++++++++++++++++ .../tests/unit/metrics/test_tcr.py | 36 +++++ .../tests/unit/test_provider.py | 10 +- 8 files changed, 317 insertions(+), 5 deletions(-) create mode 100644 changelog/62.feature.md create mode 100644 packages/ref-metrics-esmvaltool/src/cmip_ref_metrics_esmvaltool/metrics/tcr.py create mode 100644 packages/ref-metrics-esmvaltool/tests/unit/metrics/input_files_tcr.json create mode 100644 packages/ref-metrics-esmvaltool/tests/unit/metrics/test_tcr.py diff --git a/changelog/62.feature.md b/changelog/62.feature.md new file mode 100644 index 0000000..9bb1206 --- /dev/null +++ b/changelog/62.feature.md @@ -0,0 +1 @@ +Added Transient Climate Response (TCS) to the ESMValTool metrics package. diff --git a/packages/ref-metrics-esmvaltool/src/cmip_ref_metrics_esmvaltool/__init__.py b/packages/ref-metrics-esmvaltool/src/cmip_ref_metrics_esmvaltool/__init__.py index 5096bf0..740f06e 100644 --- a/packages/ref-metrics-esmvaltool/src/cmip_ref_metrics_esmvaltool/__init__.py +++ b/packages/ref-metrics-esmvaltool/src/cmip_ref_metrics_esmvaltool/__init__.py @@ -2,11 +2,12 @@ Rapid evaluating CMIP data with ESMValTool. """ +import cmip_ref_metrics_esmvaltool.metrics from cmip_ref_core.providers import MetricsProvider from cmip_ref_metrics_esmvaltool._version import __version__ -from cmip_ref_metrics_esmvaltool.metrics import EquilibriumClimateSensitivity, GlobalMeanTimeseries -# Initialise the metrics manager and register the example metric +# Initialise the metrics manager and register the metrics provider = MetricsProvider("ESMValTool", __version__) -provider.register(GlobalMeanTimeseries()) -provider.register(EquilibriumClimateSensitivity()) +for _metric_cls_name in cmip_ref_metrics_esmvaltool.metrics.__all__: + _metric_cls = getattr(cmip_ref_metrics_esmvaltool.metrics, _metric_cls_name) + provider.register(_metric_cls()) diff --git a/packages/ref-metrics-esmvaltool/src/cmip_ref_metrics_esmvaltool/metrics/__init__.py b/packages/ref-metrics-esmvaltool/src/cmip_ref_metrics_esmvaltool/metrics/__init__.py index fc35d9d..21ada81 100644 --- a/packages/ref-metrics-esmvaltool/src/cmip_ref_metrics_esmvaltool/metrics/__init__.py +++ b/packages/ref-metrics-esmvaltool/src/cmip_ref_metrics_esmvaltool/metrics/__init__.py @@ -2,8 +2,10 @@ from cmip_ref_metrics_esmvaltool.metrics.ecs import EquilibriumClimateSensitivity from cmip_ref_metrics_esmvaltool.metrics.example import GlobalMeanTimeseries +from cmip_ref_metrics_esmvaltool.metrics.tcr import TransientClimateResponse __all__ = [ "EquilibriumClimateSensitivity", "GlobalMeanTimeseries", + "TransientClimateResponse", ] diff --git a/packages/ref-metrics-esmvaltool/src/cmip_ref_metrics_esmvaltool/metrics/tcr.py b/packages/ref-metrics-esmvaltool/src/cmip_ref_metrics_esmvaltool/metrics/tcr.py new file mode 100644 index 0000000..d50b6d9 --- /dev/null +++ b/packages/ref-metrics-esmvaltool/src/cmip_ref_metrics_esmvaltool/metrics/tcr.py @@ -0,0 +1,117 @@ +from pathlib import Path + +import pandas +import xarray + +from cmip_ref_core.datasets import FacetFilter, SourceDatasetType +from cmip_ref_core.metrics import DataRequirement +from cmip_ref_metrics_esmvaltool._version import __version__ +from cmip_ref_metrics_esmvaltool.metrics.base import ESMValToolMetric +from cmip_ref_metrics_esmvaltool.recipe import dataframe_to_recipe +from cmip_ref_metrics_esmvaltool.types import OutputBundle, Recipe + + +class TransientClimateResponse(ESMValToolMetric): + """ + Calculate the global mean transient climate response for a dataset. + """ + + name = "Transient Climate Response" + slug = "esmvaltool-transient-climate-response" + base_recipe = "recipe_tcr.yml" + + data_requirements = ( + DataRequirement( + source_type=SourceDatasetType.CMIP6, + filters=( + FacetFilter( + facets={ + "variable_id": ("tas",), + "experiment_id": ( + "1pctCO2", + "piControl", + ), + }, + ), + ), + # TODO: Select only datasets that have both experiments + # TODO: Select only datasets that have a contiguous, shared timerange + # TODO: Add cell areas to the groups + # constraints=(AddCellAreas(),), + group_by=("source_id", "variant_label"), + ), + ) + + @staticmethod + def update_recipe(recipe: Recipe, input_files: pandas.DataFrame) -> None: + """Update the recipe.""" + # Only run the diagnostic that computes TCR for a single model. + recipe["diagnostics"] = { + "cmip6": { + "description": "Calculate TCR.", + "variables": { + "tas": { + "preprocessor": "spatial_mean", + }, + }, + "scripts": { + "tcr": { + "script": "climate_metrics/tcr.py", + "calculate_mmm": False, + }, + }, + }, + } + + # Prepare updated datasets section in recipe. It contains two + # datasets, one for the "1pctCO2" and one for the "piControl" + # experiment. + recipe_variables = dataframe_to_recipe(input_files) + + # Select a timerange covered by all datasets. + start_times, end_times = [], [] + for variable in recipe_variables.values(): + for dataset in variable["additional_datasets"]: + start, end = dataset["timerange"].split("/") + start_times.append(start) + end_times.append(end) + timerange = f"{max(start_times)}/{min(end_times)}" + + datasets = recipe_variables["tas"]["additional_datasets"] + for dataset in datasets: + dataset["timerange"] = timerange + + recipe["datasets"] = datasets + + @staticmethod + def format_result(result_dir: Path) -> OutputBundle: + """Format the result.""" + tcr_file = result_dir / "work/cmip6/tcr/tcr.nc" + tcr = xarray.open_dataset(tcr_file) + + source_id = tcr.dataset.values[0].decode("utf-8") + cmec_output = { + "DIMENSIONS": { + "dimensions": { + "source_id": {source_id: {}}, + "region": {"global": {}}, + "variable": {"tcr": {}}, + }, + "json_structure": [ + "model", + "region", + "statistic", + ], + }, + # Is the schema tracked? + "SCHEMA": { + "name": "CMEC-REF", + "package": "cmip_ref_metrics_esmvaltool", + "version": __version__, + }, + "RESULTS": { + source_id: {"global": {"tcr": float(tcr.tcr.values[0])}}, + }, + } + + return cmec_output diff --git a/packages/ref-metrics-esmvaltool/src/cmip_ref_metrics_esmvaltool/recipes.txt b/packages/ref-metrics-esmvaltool/src/cmip_ref_metrics_esmvaltool/recipes.txt index 9ccd5b6..fc91638 100644 --- a/packages/ref-metrics-esmvaltool/src/cmip_ref_metrics_esmvaltool/recipes.txt +++ b/packages/ref-metrics-esmvaltool/src/cmip_ref_metrics_esmvaltool/recipes.txt @@ -1,2 +1,3 @@ examples/recipe_python.yml ab3f06d269bb2c1368f4dc39da9bcb232fb2adb1fa556ba769e6c16294ffb4a3 recipe_ecs.yml 0cc57034fcb64e32015b4ff949ece5df8cdb8c6f493618b50ceded119fb37918 +recipe_tcr.yml 35f9ef035a4e71aff5cac5dd26c49da2162fc00291bf3b0bd16b661b7b2f606b diff --git a/packages/ref-metrics-esmvaltool/tests/unit/metrics/input_files_tcr.json b/packages/ref-metrics-esmvaltool/tests/unit/metrics/input_files_tcr.json new file mode 100644 index 0000000..5a44afb --- /dev/null +++ b/packages/ref-metrics-esmvaltool/tests/unit/metrics/input_files_tcr.json @@ -0,0 +1,146 @@ +{ + "start_time": { + "94": "0101-01-16T12:00:00.000", + "92": "0101-01-16T12:00:00.000" + }, + "end_time": { + "94": "0250-12-16T12:00:00.000", + "92": "0600-12-16T12:00:00.000" + }, + "path": { + "94": "/home/bandela/climate_data/CMIP6/CMIP/CSIRO/ACCESS-ESM1-5/1pctCO2/r1i1p1f1/Amon/tas/gn/v20191115/tas_Amon_ACCESS-ESM1-5_1pctCO2_r1i1p1f1_gn_010101-025012.nc", + "92": "/home/bandela/climate_data/CMIP6/CMIP/CSIRO/ACCESS-ESM1-5/piControl/r1i1p1f1/Amon/tas/gn/v20210316/tas_Amon_ACCESS-ESM1-5_piControl_r1i1p1f1_gn_010101-060012.nc" + }, + "activity_id": { + "94": "CMIP", + "92": "CMIP" + }, + "branch_method": { + "94": "standard", + "92": "standard" + }, + "branch_time_in_child": { + "94": 0.0, + "92": 0.0 + }, + "branch_time_in_parent": { + "94": 0.0, + "92": 36524.0 + }, + "experiment": { + "94": "1 percent per year increase in CO2", + "92": "pre-industrial control" + }, + "experiment_id": { + "94": "1pctCO2", + "92": "piControl" + }, + "frequency": { + "94": "mon", + "92": "mon" + }, + "grid": { + "94": "native atmosphere N96 grid (145x192 latxlon)", + "92": "native atmosphere N96 grid (145x192 latxlon)" + }, + "grid_label": { + "94": "gn", + "92": "gn" + }, + "institution_id": { + "94": "CSIRO", + "92": "CSIRO" + }, + "nominal_resolution": { + "94": "250 km", + "92": "250 km" + }, + "parent_activity_id": { + "94": "CMIP", + "92": "CMIP" + }, + "parent_experiment_id": { + "94": "piControl", + "92": "piControl-spinup" + }, + "parent_source_id": { + "94": "ACCESS-ESM1-5", + "92": "ACCESS-ESM1-5" + }, + "parent_time_units": { + "94": "days since 0101-01-01", + "92": "days since 0001-01-01" + }, + "parent_variant_label": { + "94": "r1i1p1f1", + "92": "r1i1p1f1" + }, + "product": { + "94": "model-output", + "92": "model-output" + }, + "realm": { + "94": "atmos", + "92": "atmos" + }, + "source_id": { + "94": "ACCESS-ESM1-5", + "92": "ACCESS-ESM1-5" + }, + "source_type": { + "94": "AOGCM", + "92": "AOGCM" + }, + "sub_experiment": { + "94": "none", + "92": "none" + }, + "sub_experiment_id": { + "94": "none", + "92": "none" + }, + "table_id": { + "94": "Amon", + "92": "Amon" + }, + "variable_id": { + "94": "tas", + "92": "tas" + }, + "variant_label": { + "94": "r1i1p1f1", + "92": "r1i1p1f1" + }, + "member_id": { + "94": "r1i1p1f1", + "92": "r1i1p1f1" + }, + "standard_name": { + "94": "air_temperature", + "92": "air_temperature" + }, + "long_name": { + "94": "Near-Surface Air Temperature", + "92": "Near-Surface Air Temperature" + }, + "units": { + "94": "K", + "92": "K" + }, + "vertical_levels": { + "94": 1, + "92": 1 + }, + "init_year": { + "94": null, + "92": null + }, + "version": { + "94": "v20191115", + "92": "v20210316" + }, + "instance_id": { + "94": "CMIP6.CMIP.CSIRO.ACCESS-ESM1-5.1pctCO2.r1i1p1f1.Amon.tas.gn.v20191115", + "92": "CMIP6.CMIP.CSIRO.ACCESS-ESM1-5.piControl.r1i1p1f1.Amon.tas.gn.v20210316" + } +} diff --git a/packages/ref-metrics-esmvaltool/tests/unit/metrics/test_tcr.py b/packages/ref-metrics-esmvaltool/tests/unit/metrics/test_tcr.py new file mode 100644 index 0000000..8abf8a5 --- /dev/null +++ b/packages/ref-metrics-esmvaltool/tests/unit/metrics/test_tcr.py @@ -0,0 +1,36 @@ +from pathlib import Path + +import numpy as np +import pandas +import xarray as xr +from cmip_ref_metrics_esmvaltool.metrics import TransientClimateResponse +from cmip_ref_metrics_esmvaltool.recipe import load_recipe + + +def test_update_recipe(): + input_files = pandas.read_json(Path(__file__).parent / "input_files_tcr.json") + recipe = load_recipe("recipe_tcr.yml") + TransientClimateResponse().update_recipe(recipe, input_files) + assert len(recipe["datasets"]) == 2 + assert len(recipe["diagnostics"]) == 1 + assert set(recipe["diagnostics"]["cmip6"]["variables"]) == {"tas"} + + +def test_format_output(tmp_path): + tcr = xr.Dataset( + data_vars={ + "tcr": (["dim0"], np.array([1.0], dtype=np.float32)), + }, + coords={ + "dataset": ("dim0", np.array([b"abc"])), + }, + ) + result_dir = tmp_path + subdir = result_dir / "work" / "cmip6" / "tcr" + subdir.mkdir(parents=True) + tcr.to_netcdf(subdir / "tcr.nc") + + output_bundle = TransientClimateResponse().format_result(result_dir) + + assert isinstance(output_bundle, dict) + assert output_bundle["RESULTS"]["abc"]["global"]["tcr"] == 1.0 diff --git a/packages/ref-metrics-esmvaltool/tests/unit/test_provider.py b/packages/ref-metrics-esmvaltool/tests/unit/test_provider.py index 1159b1f..e6ad244 100644 --- a/packages/ref-metrics-esmvaltool/tests/unit/test_provider.py +++ b/packages/ref-metrics-esmvaltool/tests/unit/test_provider.py @@ -1,3 +1,5 @@ +import importlib.metadata + from cmip_ref_metrics_esmvaltool import __version__, provider @@ -6,4 +8,10 @@ def test_provider(): assert provider.slug == "esmvaltool" assert provider.version == __version__ - assert len(provider) == 2 + metric_modules = importlib.resources.files("cmip_ref_metrics_esmvaltool").glob("metrics/*.py") + ignore = { + "__init__.py", + "base.py", + } + n_metric_modules = len([f for f in metric_modules if f.name not in ignore]) + assert len(provider) == n_metric_modules