Add sea ice sensitivity diagnostic

changelog/378.feature.md

@@ -0,0 +1 @@
+Added sea ice sensitivity diagnostic.

packages/climate-ref-core/src/climate_ref_core/constraints.py

@@ -200,17 +200,27 @@ def apply(
         for facet, values in supplementary_facets.items():
             mask = supplementary_group[facet].isin(values)
             supplementary_group = supplementary_group[mask]
-
-        if not supplementary_group.empty and self.optional_matching_facets:
-            facets = list(self.matching_facets + self.optional_matching_facets)
+        if not supplementary_group.empty:
+            matching_facets = list(self.matching_facets)
+            facets = matching_facets + list(self.optional_matching_facets)
             datasets = group[facets].drop_duplicates()
             indices = set()
             for i in range(len(datasets)):
-                scores = (supplementary_group[facets] == datasets.iloc[i]).sum(axis=1)
-                matches = supplementary_group[scores == scores.max()]
-                # Select the latest version if there are multiple matches
-                matches = matches[matches["version"] == matches["version"].max()]
-                indices.add(matches.index[0])
+                dataset = datasets.iloc[i]
+                # Restrict the supplementary datasets to those that match the main dataset.
+                supplementaries = supplementary_group[
+                    (supplementary_group[matching_facets] == dataset[matching_facets]).all(1)
+                ]
+                if not supplementaries.empty:
+                    # Select the best matching supplementary dataset based on the optional matching facets.
+                    scores = (supplementaries[facets] == dataset).sum(axis=1)
+                    matches = supplementaries[scores == scores.max()]
+                    if "version" in facets:
+                        # Select the latest version if there are multiple matches
+                        matches = matches[matches["version"] == matches["version"].max()]
+                    # Select one match per dataset
+                    indices.add(matches.index[0])
+
             supplementary_group = supplementary_group.loc[list(indices)].drop_duplicates()
 
         return pd.concat([group, supplementary_group])

packages/climate-ref-core/tests/unit/test_constraints.py

@@ -60,6 +60,21 @@ class TestAddSupplementaryDataset:
                 ),
                 [0],
             ),
+            (
+                # Test that missing supplementary files are handled gracefully.
+                pd.DataFrame(
+                    {
+                        "variable_id": ["tas", "areacella", "areacella", "tas"],
+                        "source_id": ["X", "X", "X", "Y"],
+                        "grid_label": ["gn"] * 4,
+                        "table_id": ["Amon", "fx", "fx", "Amon"],
+                        "experiment_id": ["historical"] * 4,
+                        "member_id": ["r1i1p1f1", "r1i1p1f1", "r2i1p1f1", "r2i1p1f1"],
+                        "version": ["v20210316", "v20210316", "v20210317", "v20210317"],
+                    }
+                ),
+                [0, 3, 1],
+            ),
             (
                 # Test that the grid_label matches.
                 pd.DataFrame(
@@ -68,7 +83,7 @@ class TestAddSupplementaryDataset:
                         "source_id": ["ACCESS-ESM1-5"] * 3,
                         "grid_label": ["gn", "gn", "gr"],
                         "table_id": ["Amon", "fx", "fx"],
-                        "experiment_id": ["historical"] * 3,
+                        "experiment_id": ["historical", "piControl", "historical"],
                         "member_id": ["r1i1p1f1", "r2i1p1f1", "r1i1p1f1"],
                         "version": ["v20210316"] * 3,
                     }

packages/climate-ref-esmvaltool/pyproject.toml

@@ -25,7 +25,7 @@ classifiers = [
 dependencies = [
     "pooch >= 1.8",
     "climate-ref-core",
-    "ruamel.yaml >= 0.18",
+    "pyyaml",
     "xarray >= 2023.3.0",
 ]
 [project.entry-points."climate-ref.providers"]

packages/climate-ref-esmvaltool/src/climate_ref_esmvaltool/diagnostics/__init__.py

@@ -6,6 +6,7 @@
 from climate_ref_esmvaltool.diagnostics.enso import ENSOBasicClimatology, ENSOCharacteristics
 from climate_ref_esmvaltool.diagnostics.example import GlobalMeanTimeseries
 from climate_ref_esmvaltool.diagnostics.sea_ice_area_basic import SeaIceAreaBasic
+from climate_ref_esmvaltool.diagnostics.sea_ice_sensitivity import SeaIceSensitivity
 from climate_ref_esmvaltool.diagnostics.tcr import TransientClimateResponse
 from climate_ref_esmvaltool.diagnostics.tcre import TransientClimateResponseEmissions
 from climate_ref_esmvaltool.diagnostics.zec import ZeroEmissionCommitment
@@ -18,6 +19,7 @@
     "EquilibriumClimateSensitivity",
     "GlobalMeanTimeseries",
     "SeaIceAreaBasic",
+    "SeaIceSensitivity",
     "TransientClimateResponse",
     "TransientClimateResponseEmissions",
     "ZeroEmissionCommitment",

packages/climate-ref-esmvaltool/src/climate_ref_esmvaltool/diagnostics/base.py

@@ -4,8 +4,8 @@
 from typing import ClassVar
 
 import pandas
+import yaml
 from loguru import logger
-from ruamel.yaml import YAML
 
 from climate_ref_core.dataset_registry import dataset_registry_manager
 from climate_ref_core.datasets import ExecutionDatasetCollection, SourceDatasetType
@@ -19,8 +19,6 @@
 from climate_ref_esmvaltool.recipe import load_recipe, prepare_climate_data
 from climate_ref_esmvaltool.types import MetricBundleArgs, OutputBundleArgs, Recipe
 
-yaml = YAML()
-
 
 class ESMValToolDiagnostic(CommandLineDiagnostic):
     """ESMValTool Diagnostic base class."""
@@ -177,7 +175,7 @@ def build_execution_result(
         # Add the plots and data files
         plot_suffixes = {".png", ".jpg", ".pdf", ".ps"}
         for metadata_file in result_dir.glob("run/*/*/diagnostic_provenance.yml"):
-            metadata = yaml.load(metadata_file.read_text(encoding="utf-8"))
+            metadata = yaml.safe_load(metadata_file.read_text(encoding="utf-8"))
             for filename in metadata:
                 caption = metadata[filename].get("caption", "")
                 relative_path = definition.as_relative_path(filename)

packages/climate-ref-esmvaltool/src/climate_ref_esmvaltool/diagnostics/sea_ice_sensitivity.py

@@ -0,0 +1,105 @@
+from pathlib import Path
+
+import pandas
+import pandas as pd
+
+from climate_ref_core.constraints import (
+    AddSupplementaryDataset,
+    RequireContiguousTimerange,
+    RequireFacets,
+)
+from climate_ref_core.datasets import ExecutionDatasetCollection, FacetFilter, SourceDatasetType
+from climate_ref_core.diagnostics import DataRequirement
+from climate_ref_core.pycmec.metric import CMECMetric, MetricCV
+from climate_ref_core.pycmec.output import CMECOutput
+from climate_ref_esmvaltool.diagnostics.base import ESMValToolDiagnostic
+from climate_ref_esmvaltool.recipe import dataframe_to_recipe
+from climate_ref_esmvaltool.types import MetricBundleArgs, OutputBundleArgs, Recipe
+
+
+class SeaIceSensitivity(ESMValToolDiagnostic):
+    """
+    Calculate sea ice sensitivity.
+    """
+
+    name = "Sea ice sensitivity"
+    slug = "sea-ice-sensitivity"
+    base_recipe = "recipe_seaice_sensitivity.yml"
+
+    variables = (
+        "siconc",
+        "tas",
+    )
+
+    data_requirements = (
+        DataRequirement(
+            source_type=SourceDatasetType.CMIP6,
+            filters=(
+                FacetFilter(
+                    facets={
+                        "variable_id": variables,
+                        "experiment_id": "historical",
+                    },
+                ),
+            ),
+            group_by=("experiment_id",),  # this does nothing, but group_by cannot be empty
+            constraints=(
+                AddSupplementaryDataset.from_defaults("areacella", SourceDatasetType.CMIP6),
+                AddSupplementaryDataset.from_defaults("areacello", SourceDatasetType.CMIP6),
+                RequireContiguousTimerange(group_by=("instance_id",)),
+                RequireFacets("variable_id", variables),
+                # TODO: Add a constraint to ensure that tas, siconc and areacello
+                # are available for each model or alternatively filter out
+                # incomplete models below.
+            ),
+        ),
+    )
+    facets = ("experiment_id", "source_id", "region", "metric")
+
+    @staticmethod
+    def update_recipe(recipe: Recipe, input_files: pandas.DataFrame) -> None:
+        """Update the recipe."""
+        recipe_variables = dataframe_to_recipe(input_files)
+        datasets = recipe_variables["tas"]["additional_datasets"]
+        for dataset in datasets:
+            dataset.pop("mip")
+            dataset["timerange"] = "1979/2014"
+        recipe["datasets"] = datasets
+
+    @staticmethod
+    def format_result(
+        result_dir: Path,
+        execution_dataset: ExecutionDatasetCollection,
+        metric_args: MetricBundleArgs,
+        output_args: OutputBundleArgs,
+    ) -> tuple[CMECMetric, CMECOutput]:
+        """Format the result."""
+        metric_args[MetricCV.DIMENSIONS.value] = {
+            "json_structure": [
+                "source_id",
+                "region",
+                "metric",
+            ],
+            "source_id": {},
+            "region": {},
+            "metric": {},
+        }
+        for region in "antarctic", "arctic":
+            df = pd.read_csv(
+                result_dir / "work" / region / "sea_ice_sensitivity_script" / "plotted_values.csv"
+            )
+            df = df.rename(columns={"Unnamed: 0": "source_id"}).drop(columns=["label"])
+            metric_args[MetricCV.DIMENSIONS.value]["region"][region] = {}
+            for metric in df.columns[1:]:
+                metric_args[MetricCV.DIMENSIONS.value]["metric"][metric] = {}
+            for row in df.itertuples(index=False):
+                source_id = row.source_id
+                metric_args[MetricCV.DIMENSIONS.value]["source_id"][source_id] = {}
+                for metric, value in zip(df.columns[1:], row[1:]):
+                    if source_id not in metric_args[MetricCV.RESULTS.value]:
+                        metric_args[MetricCV.RESULTS.value][source_id] = {}
+                    if region not in metric_args[MetricCV.RESULTS.value][source_id]:
+                        metric_args[MetricCV.RESULTS.value][source_id][region] = {}
+                    metric_args[MetricCV.RESULTS.value][source_id][region][metric] = value
+
+        return CMECMetric.model_validate(metric_args), CMECOutput.model_validate(output_args)

packages/climate-ref-esmvaltool/src/climate_ref_esmvaltool/recipe.py

@@ -5,15 +5,13 @@
 from typing import TYPE_CHECKING, Any
 
 import pooch
-from ruamel.yaml import YAML
+import yaml
 
 from climate_ref_esmvaltool.types import Recipe
 
 if TYPE_CHECKING:
     import pandas as pd
 
-yaml = YAML()
-
 FACETS = {
     "CMIP6": {
         "activity": "activity_id",
@@ -115,8 +113,8 @@ def dataframe_to_recipe(files: pd.DataFrame) -> dict[str, Any]:
     return variables
 
 
-_ESMVALTOOL_COMMIT = "58fd0b8ece981bc97c4fbd213b11f2228d90db28"
-_ESMVALTOOL_VERSION = f"2.13.0.dev65+g{_ESMVALTOOL_COMMIT[:9]}"
+_ESMVALTOOL_COMMIT = "8f56863a70ba4df76ec501ba0372c571a0af6cf9"
+_ESMVALTOOL_VERSION = f"2.13.0.dev120+g{_ESMVALTOOL_COMMIT[:9]}"
 
 _RECIPES = pooch.create(
     path=pooch.os_cache("climate_ref_esmvaltool"),
@@ -144,7 +142,7 @@ def load_recipe(recipe: str) -> Recipe:
         The loaded recipe.
     """
     filename = _RECIPES.fetch(recipe)
-    return yaml.load(Path(filename).read_text(encoding="utf-8"))  # type: ignore[no-any-return]
+    return yaml.safe_load(Path(filename).read_text(encoding="utf-8"))  # type: ignore[no-any-return]
 
 
 def prepare_climate_data(datasets: pd.DataFrame, climate_data_dir: Path) -> None:

packages/climate-ref-esmvaltool/src/climate_ref_esmvaltool/recipes.txt

@@ -1,10 +1,11 @@
 examples/recipe_python.yml ab3f06d269bb2c1368f4dc39da9bcb232fb2adb1fa556ba769e6c16294ffb4a3
 recipe_calculate_gwl_exceedance_stats.yml 5aa266abc9a8029649b689a2b369a47623b0935d609354332ff4148994642d6b
 recipe_ecs.yml             0cc57034fcb64e32015b4ff949ece5df8cdb8c6f493618b50ceded119fb37918
+recipe_seaice_sensitivity.yml f7247c076e161c582d422947c8155f3ca98549e6f2e4c3b1c76414786d7e50c5
 recipe_tcr.yml             35f9ef035a4e71aff5cac5dd26c49da2162fc00291bf3b0bd16b661b7b2f606b
 recipe_tcre.yml            48fc9e3baf541bbcef7491853ea3a774053771dca33352b41466425faeaa38af
 recipe_zec.yml             b0af7f789b7610ab3f29a6617124aa40c40866ead958204fc199eaf82863de51
 ref/recipe_enso_basicclimatology.yml 9ea7deb7ee668e39ac44618b96496d898bd82285c22dcee4fce4695e0c9fa82b
 ref/recipe_enso_characteristics.yml 34c2518b138068ac96d212910b979d54a8fcedee2c0089b5acd56a42c41dc3e4
-ref/recipe_ref_cre.yml     4f35d9639f1008be3b5382a5bd8933a855cb5368ccf5d04a1c70227172e2e82c
+ref/recipe_ref_cre.yml     4375f262479c3b3e1b348b71080a6d758e195bda76516a591182045a3a29aa32
 ref/recipe_ref_sea_ice_area_basic.yml 7d01a8527880663ca28284772f83a8356d9972fb4f022a4000e50a56ce044b09