Add CORDEX-CMIP5 evaluation recipes

[ghossh]

Continuation from #4199

Description

This pull request adds recipes to demonstrate that ESMValTool is now capable of handling CORDEX data.

Use the esmvalcore-dev environment to run the recipes.

examples/recipe_map_CORDEX-CMIP5.yml

recipe_cordex-cmip5.yml plots maps for the CORDEX data and, if an observational dataset is present that is marked with the facet reference_for_monitor_diags: true, computes the bias against it. The maps are useful to visually inspect the data and detect issues in the values or the grid. Some example plots that can be produced are the following (with reference_for_monitor_diags: true on the observational data):

and (without reference_for_monitor_diags: true on the observational data):

All figures with comparison to CCI data are available at https://esmvaltool.dkrz.de/shared/esmvaltool/cordex/recipe_map_CORDEX-CMIP5_20260713_151914/ and figures with individual plots are available at https://esmvaltool.dkrz.de/shared/esmvaltool/cordex/recipe_map_CORDEX-CMIP5_20260714_084936/

Code to generate the recipe:

Details

"""Script to generate a recipe that plots maps of CORDEX model vs reference data."""

import copy

from esmvalcore._recipe.writer import to_yaml
from esmvalcore.config import CFG
from esmvalcore.dataset import Dataset, datasets_to_recipe


def main() -> None:
    """Make a recipe with all CORDEX data we can find."""
    # Make sure to configure the data sources before running this script.
    # To configure ESGF as a data source, run:
    # `esmvaltool config copy data-esmvalcore-esgf.yml`
    CFG["search_data"] = "complete"

    # Set up the preprocessors.
    preprocessor = {
        "custom_order": True,
        "monthly_statistics": {
            "operator": "mean",
        },
        "climate_statistics": {
            "period": "full",
            "operator": "mean",
        },
        "mask_landsea": {
            "mask_out": "sea",
        },
        "regrid": {
            "target_grid": "EUR-11",
            "scheme": "linear",
        },
        "convert_units": {
            "units": "degC",
        },
        "multi_model_statistics": {
            "span": "overlap",
            "statistics": ["mean", "median"],
            "exclude": ["reference_dataset"],
        },
    }
    land_preprocessor = copy.deepcopy(preprocessor)
    land_preprocessor["mask_landsea"]["mask_out"] = "sea"
    sea_preprocessor = copy.deepcopy(preprocessor)
    sea_preprocessor["mask_landsea"]["mask_out"] = "land"
    sea_preprocessor["align_metadata"] = {
        "target_project": "CORDEX",
        "target_mip": "day",
        "target_short_name": "ts",
    }
    preprocessor.pop("mask_landsea")
    preprocessor.pop("convert_units")
    sic_preprocessor = copy.deepcopy(preprocessor)
    sic_preprocessor["align_metadata"] = {
        "target_project": "CORDEX",
        "target_mip": "day",
        "target_short_name": "sic",
    }

    documentation = {
        "title": "CORDEX-CMIP5 maps",
        "description": """This is a recipe that generates maps for CORDEX-CMIP5 models. It computes the climatology and plots the bias against an ESA CCI dataset.""",
        "authors": ["loosveldt-tomas_saskia", "andela_bouwe", "ghosh_supriyo"],
        "maintainer": ["ghosh_supriyo", "andela_bouwe"],
    }

    # Define the datasets to include in the recipe.
    reference_datasets = {
        "clivi": Dataset(
            short_name="clivi",
            mip="Amon",
            project="OBS6",
            tier=2,
            dataset="ESACCI-CLOUD",
            version="v3.0-AVHRR-AMPM",
            type="sat",
        ),
        "clt": Dataset(
            short_name="clt",
            mip="day",
            project="OBS6",
            tier=2,
            dataset="ESACCI-CLOUD",
            version="v3.0-AVHRR-AMPM",
            type="sat",
        ),
        "lwp": Dataset(
            short_name="lwp",
            mip="Amon",
            project="OBS6",
            tier=2,
            dataset="ESACCI-CLOUD",
            version="v3.0-AVHRR-AMPM",
            type="sat",
            derive=False,
        ),
        "prw": Dataset(
            short_name="prw",
            mip="Eday",
            project="OBS6",
            tier=3,
            dataset="ESACCI-WATERVAPOUR",
            version="CDR2-L3-COMBI-05deg-fv3.1",
            type="sat",
        ),
        "rlut": Dataset(
            short_name="rlut",
            mip="Amon",
            project="OBS6",
            tier=2,
            dataset="ESACCI-CLOUD",
            version="v3.0-AVHRR-AMPM",
            type="sat",
        ),
        "rsus": Dataset(
            short_name="rsus",
            mip="Amon",
            project="OBS6",
            tier=2,
            dataset="ESACCI-CLOUD",
            type="sat",
            version="v3.0-AVHRR-AMPM",
        ),
        "sic": Dataset(
            short_name="sic",
            mip="SIday",
            project="OBS6",
            tier=2,
            dataset="ESACCI-SEAICE",
            version="L4-SICONC-RE-SSMI-12.5kmEASE2-fv3.0-NH",
            type="sat",
        ),
        "snw": Dataset(
            short_name="snw",
            mip="day",
            project="OBS6",
            tier=2,
            dataset="ESACCI-SNOW",
            version="v2.0",
            type="sat",
        ),
        "lst": Dataset(
            short_name="ts",
            mip="Amon",
            project="OBS",
            tier=2,
            dataset="ESACCI-LST",
            version="1.00",
            type="sat",
        ),
        "sst": Dataset(
            short_name="tos",
            mip="Oday",
            project="OBS6",
            tier=2,
            dataset="ESACCI-SST",
            version="3.0-L4-analysis",
            type="sat",
        ),
    }
    reference_datasets["lst"].add_supplementary(
        short_name="sftlf",
        skip=True,
    )
    reference_datasets["sst"].add_supplementary(
        short_name="sftof",
        skip=True,
    )
    for variable_group, dataset in reference_datasets.items():
        dataset.set_facet("reference_for_monitor_diags", True, persist=True)
        dataset.set_facet("variable_group", variable_group, persist=False)

    datasets = list(reference_datasets.values())
    for short_name in [
        "clivi",
        "clt",
        "clwvi",  # "lwp" is derived from "clwvi" and "clivi"
        "prw",
        "rsus",
        "rlut",
        "sic",
        "snw",
        "ts",
    ]:
        template = Dataset(
            project="CORDEX",
            mip="day",
            short_name=short_name,
            domain="EUR-11",
            dataset="*",
            exp="historical",
            institute="*",
            rcm_version="*",
            driver="*",
            ensemble="*",
        )
        for dataset in template.from_files():
            if dataset.facets["dataset"] == "RegCM4-6":
                continue  # Data no longer available from Cineca ESGF node. -> Some left on Levante, but would need fixes.
            if dataset.facets["short_name"] == "clwvi":
                # "lwp" is derived from "clwvi" and "clivi", but "clwvi" is
                # not required for the recipe.
                dataset.facets["short_name"] = "lwp"
            if dataset.facets["short_name"] == "ts":
                # "ts" needs to be separated into land and sea surface
                # temperature in order to compare with LST and SST reference
                # datasets.
                dataset.add_supplementary(
                    short_name="sftlf",
                    mip="fx",
                    ensemble="*",
                    exp="*",
                    # Not entirely sure if driver and institute should be
                    # wildcards here or if grids differ between drivers and
                    # institutes.
                    driver="*",
                    institute="*",
                )
                result = next(dataset.from_files())
                if result.supplementaries:
                    if (
                        result.supplementaries[0].facets["driver"]
                        == result.facets["driver"]
                    ):
                        # do not include the supplementary in the recipe explicitly, it will be added automatically.
                        result.supplementaries.clear()
                    lst = result.copy(preprocessor="land")
                    lst.set_facet("variable_group", "lst", persist=False)
                    datasets.append(lst)
                    sst = result.copy(preprocessor="sea")
                    sst.set_facet("variable_group", "sst", persist=False)
                    datasets.append(sst)
            else:
                datasets.append(dataset)

    for dataset in datasets:
        dataset.set_facet(
            "diagnostic",
            dataset.facets.get("variable_group", dataset.facets["short_name"]),
            persist=False,
        )

    recipe = datasets_to_recipe(datasets)
    recipe["documentation"] = documentation
    recipe["preprocessors"] = {
        "default": preprocessor,
        "land": land_preprocessor,
        "sea": sea_preprocessor,
        "sic": sic_preprocessor,
    }

    for diagnostic_name, diagnostic in recipe["diagnostics"].items():
        # Configure reference datasets.
        for variable_group, variable in diagnostic["variables"].items():
            if variable_group in reference_datasets:
                variable["reference_dataset"] = reference_datasets[
                    variable_group
                ].facets["dataset"]

        # Move common facets to the variable level.
        for variable_group, variable in diagnostic["variables"].items():
            variable["project"] = "CORDEX"
            variable["mip"] = "day"
            variable["domain"] = "EUR-11"
            variable["exp"] = "historical"
            variable["ensemble"] = "r1i1p1"
            variable["rcm_version"] = "v1"
            variable["timerange"] = "2003/2005"
            if variable_group == "lwp":
                variable["derive"] = True
            if variable_group == "lst":
                variable["preprocessor"] = "land"
                variable["short_name"] = "ts"
            elif variable_group == "sst":
                variable["preprocessor"] = "sea"
                variable["short_name"] = "ts"
            elif variable_group == "sic":
                variable["preprocessor"] = "sic"
            else:
                variable["preprocessor"] = "default"
            for dataset in diagnostic["additional_datasets"]:
                for key in list(dataset.keys()):
                    if dataset[key] == variable.get(key):
                        dataset.pop(key)

        # Set up the diagnostic script.
        diagnostic["scripts"] = {
            "plot": {
                "plot_folder": "{plot_dir}",
                "plot_filename": "{plot_type}_{variable_group}_{alias}_{mip}",
                "facet_used_for_labels": "alias",
                "script": "monitor/multi_datasets.py",
                "plots": {
                    "map": {
                        "common_cbar": False,
                        "fontsize": 6,
                    },
                },
            },
        }

    # Move datasets that are common to all variables to the top level.
    first_datasets = recipe["diagnostics"]["clivi"]["additional_datasets"]
    for dataset in first_datasets:
        if all(
            dataset
            in diagnostic["additional_datasets"]
            for diagnostic in recipe["diagnostics"].values()
        ):
            recipe["datasets"].append(dataset)
            for diagnostic in recipe["diagnostics"].values():
                diagnostic["additional_datasets"].remove(dataset)

    print(to_yaml(recipe))


if __name__ == "__main__":
    main()

recipe_perfmetrics_CORDEX-CMIP5.yml

Produces the performance metrics plot for several variables in CORDEX and ESA-CCI data, which has been previously examined using the recipe to produce maps to ensure the data is on the correct domain. An attempt has been made to correct issues with the units and to leave out datasets with invalid numbers (e.g. all zeros or of the order of 1e36), but it is not straightforward in all cases what the correct interpretation of the data is. For some models, the data is orders of magnitude removed from the other models and ESA CCI data.

The run including provenance and log files is available at https://esmvaltool.dkrz.de/shared/esmvaltool/cordex/recipe_perfmetrics_CORDEX-CMIP5_20260716_083604/

Code to generate the recipe:

Details

"""Script to generate a performance metrics recipe for CORDEX-CMIP5 data."""

import copy
import textwrap

from esmvalcore._recipe.writer import to_yaml
from esmvalcore.config import CFG
from esmvalcore.dataset import Dataset, datasets_to_recipe


def main() -> None:
    """Make a recipe with all CORDEX data we can find."""
    # Make sure to configure the data sources before running this script.
    # To configure ESGF as a data source, run:
    # `esmvaltool config copy data-esmvalcore-esgf.yml`
    CFG["search_data"] = "complete"

    # Set up the documentation and preprocessors.
    preprocessor = {
        "custom_order": True,
        "monthly_statistics": {
            "operator": "mean",
        },
        "climate_statistics": {
            "period": "full",
            "operator": "mean",
        },
        "mask_landsea": {
            "mask_out": "sea",
        },
        "regrid": {
            "target_grid": "EUR-11",
            "scheme": "linear",
        },
        "convert_units": {
            "units": "degC",
        },
        "distance_metric": {
            "metric": "weighted_rmse",
        },
        "multi_model_statistics": {
            "span": "overlap",
            "statistics": [
                {"operator": "mean"},
                {"operator": "median"},
            ],
            "exclude": ["reference_dataset"],
        },
    }
    land_preprocessor = copy.deepcopy(preprocessor)
    land_preprocessor["mask_landsea"]["mask_out"] = "sea"
    sea_preprocessor = copy.deepcopy(preprocessor)
    sea_preprocessor["mask_landsea"]["mask_out"] = "land"
    preprocessor.pop("mask_landsea")
    preprocessor.pop("convert_units")

    documentation = {
        "title": "CORDEX-CMIP5 performance metrics",
        "description": """This is a recipe that computes performance metrics for CORDEX-CMIP5 models.""",
        "authors": ["loosveldt-tomas_saskia", "andela_bouwe", "ghosh_supriyo", "lauer_axel"],
        "maintainer": ["ghosh_supriyo", "andela_bouwe"],
    }

    # Define the datasets to include in the recipe.
    reference_datasets = {
        "clivi": Dataset(
            short_name="clivi",
            mip="Amon",
            project="OBS6",
            tier=2,
            dataset="ESACCI-CLOUD",
            version="v3.0-AVHRR-AMPM",
            type="sat",
        ),
        "clt": Dataset(
            short_name="clt",
            mip="day",
            project="OBS6",
            tier=2,
            dataset="ESACCI-CLOUD",
            version="v3.0-AVHRR-AMPM",
            type="sat",
        ),
        "lwp": Dataset(
            short_name="lwp",
            mip="Amon",
            project="OBS6",
            tier=2,
            dataset="ESACCI-CLOUD",
            version="v3.0-AVHRR-AMPM",
            type="sat",
            derive=False,
        ),
        "prw": Dataset(
            short_name="prw",
            mip="Eday",
            project="OBS6",
            tier=3,
            dataset="ESACCI-WATERVAPOUR",
            version="CDR2-L3-COMBI-05deg-fv3.1",
            type="sat",
        ),
        "rlut": Dataset(
            short_name="rlut",
            mip="Amon",
            project="OBS6",
            tier=2,
            dataset="ESACCI-CLOUD",
            version="v3.0-AVHRR-AMPM",
            type="sat",
        ),
        "rsus": Dataset(
            short_name="rsus",
            mip="Amon",
            project="OBS6",
            tier=2,
            dataset="ESACCI-CLOUD",
            type="sat",
            version="v3.0-AVHRR-AMPM",
        ),
        "sic": Dataset(
            short_name="sic",
            mip="SIday",
            project="OBS6",
            tier=2,
            dataset="ESACCI-SEAICE",
            version="L4-SICONC-RE-SSMI-12.5kmEASE2-fv3.0-NH",
            type="sat",
        ),
        "snw": Dataset(
            short_name="snw",
            mip="day",
            project="OBS6",
            tier=2,
            dataset="ESACCI-SNOW",
            version="v2.0",
            type="sat",
        ),
        "lst": Dataset(
            short_name="ts",
            mip="Amon",
            project="OBS",
            tier=2,
            dataset="ESACCI-LST",
            version="1.00",
            type="sat",
        ),
        "sst": Dataset(
            short_name="tos",
            mip="Oday",
            project="OBS6",
            tier=2,
            dataset="ESACCI-SST",
            version="3.0-L4-analysis",
            type="sat",
        ),
    }
    reference_datasets["lst"].add_supplementary(
        short_name="sftlf",
        skip=True,
    )
    reference_datasets["sst"].add_supplementary(
        short_name="sftof",
        skip=True,
    )
    for variable_group, dataset in reference_datasets.items():
        dataset.set_facet("reference_for_metric", True, persist=True)
        dataset.set_facet("variable_group", variable_group, persist=False)

    datasets = list(reference_datasets.values())
    for short_name in [
        "clivi",
        "clt",
        "clwvi",  # "lwp" is derived from "clwvi" and "clivi"
        "prw",
        "rsus",
        "rlut",
        "sic",
        "snw",
        "ts",
    ]:
        template = Dataset(
            project="CORDEX",
            mip="day",
            short_name=short_name,
            domain="EUR-11",
            dataset="*",
            exp="historical",
            institute="*",
            rcm_version="*",
            driver="*",
            ensemble="*",
        )
        for dataset in template.from_files():
            if dataset.facets["dataset"] == "RegCM4-6":
                continue  # Data no longer available from Cineca ESGF node.
            if (
                dataset.facets["dataset"] == "HIRHAM5"
                and dataset.facets["short_name"] in {"clivi", "prw"}
            ):
                # HIRHAM5 clivi and prw data is much higher than expected.
                continue
            if (
                dataset.facets["dataset"] == "CCLM4-8-17"
                and dataset.facets["short_name"] == "clt"
            ):
                # CCLM4-8-17 clt data is unusable.
                continue
            if (
                dataset.facets["dataset"] == "CCLM4-8-17"
                and dataset.facets["short_name"] == "rsus"
                and dataset.facets["driver"] in {"CNRM-CERFACS-CNRM-CM5", "ICHEC-EC-EARTH", "MOHC-HadGEM2", "MPI-M-MPI-ESM-LR"}
            ):
                # Most CCLM4-8-17 rsus data is unusable.
                continue
            if (
                dataset.facets["dataset"] == "ALADIN63"
                and dataset.facets["short_name"] == "sic"
            ):
                # ALADIN63 sic data is all zeros.
                continue
            if (
                dataset.facets["dataset"] == "RACMO22E"
                and dataset.facets["short_name"] == "snw"
            ):
                # RACMO22E snw data is has a 10 times higher mean than expected.
                continue
            if dataset.facets["short_name"] == "clwvi":
                # "lwp" is derived from "clwvi" and "clivi", but "clwvi" is
                # not required for the recipe.
                dataset.facets["short_name"] = "lwp"
            if dataset.facets["short_name"] == "ts":
                # "ts" needs to be separated into land and sea surface
                # temperature in order to compare with LST and SST reference
                # datasets.
                dataset.add_supplementary(
                    short_name="sftlf",
                    mip="fx",
                    ensemble="*",
                    exp="*",
                    # Not entirely sure if driver and institute should be
                    # wildcards here or if grids differ between drivers and
                    # institutes.
                    driver="*",
                    institute="*",
                )
                result = next(dataset.from_files())
                if result.supplementaries:
                    if (
                        result.supplementaries[0].facets["driver"]
                        == result.facets["driver"]
                    ):
                        # do not include the supplementary in the recipe explicitly, it will be added automatically.
                        result.supplementaries.clear()
                    lst = result.copy(preprocessor="land")
                    lst.set_facet("variable_group", "lst", persist=False)
                    datasets.append(lst)
                    sst = result.copy(preprocessor="sea")
                    sst.set_facet("variable_group", "sst", persist=False)
                    datasets.append(sst)
            else:
                datasets.append(dataset)

    for dataset in datasets:
        dataset.set_facet(
            "diagnostic",
            "diagnostic",
            persist=False,
        )

    recipe = datasets_to_recipe(datasets)
    recipe["documentation"] = documentation
    recipe["preprocessors"] = {
        "default": preprocessor,
        "land": land_preprocessor,
        "sea": sea_preprocessor,
    }

    # Configure reference datasets.
    for variable_group, variable in recipe["diagnostics"]["diagnostic"][
        "variables"
    ].items():
        if variable_group in reference_datasets:
            variable["reference_dataset"] = reference_datasets[
                variable_group
            ].facets["dataset"]

    # Move common facets to the variable level.
    for variable_group, variable in recipe["diagnostics"]["diagnostic"][
        "variables"
    ].items():
        variable["project"] = "CORDEX"
        variable["mip"] = "day"
        variable["domain"] = "EUR-11"
        variable["exp"] = "historical"
        variable["ensemble"] = "r1i1p1"
        variable["rcm_version"] = "v1"
        variable["timerange"] = "2003/2005"
        if variable_group == "lwp":
            variable["derive"] = True
        if variable_group == "lst":
            variable["preprocessor"] = "land"
            variable["short_name"] = "ts"
        elif variable_group == "sst":
            variable["preprocessor"] = "sea"
            variable["short_name"] = "ts"
        else:
            variable["preprocessor"] = "default"
        for dataset in variable["additional_datasets"]:
            for key in list(dataset.keys()):
                if dataset[key] == variable.get(key):
                    dataset.pop(key)

    # Move datasets that are common to all variables to the top level.
    first_datasets = next(
        iter(recipe["diagnostics"]["diagnostic"]["variables"].values())
    )["additional_datasets"]
    for dataset in first_datasets:
        if all(
            dataset
            in recipe["diagnostics"]["diagnostic"]["variables"][var][
                "additional_datasets"
            ]
            for var in recipe["diagnostics"]["diagnostic"]["variables"]
        ):
            recipe["datasets"].append(dataset)
            for var in recipe["diagnostics"]["diagnostic"]["variables"]:
                recipe["diagnostics"]["diagnostic"]["variables"][var][
                    "additional_datasets"
                ].remove(dataset)

    # Set up the diagnostic script.
    recipe["diagnostics"]["diagnostic"]["scripts"] = {
        "portrait": {
            "script": "portrait_plot.py",
            "x_by": "variable_group",
            "y_by": "alias",
            "group_by": "project",
            "figsize": [8, 8],
            "normalize": "centered_median",
            # "default_split": "Ref1",
            "nan_color": "#828282",
            "plot_kwargs": {
                "vmin": -0.5,
                "vmax": +0.5,
            },
            "cbar_kwargs": {
                "label": "Relative RMSE",
                "extend": "both",
                "aspect": 45,
            },
            "matplotlib_rc_params": {
                "ytick.labelsize": 4,
                "xtick.labelsize": 6,
                "axes.labelsize": 12,
                "font.size": 12,
            },
        },
    }

    print(to_yaml(recipe))


if __name__ == "__main__":
    main()

Closes #4198

Link to documentation:
https://esmvaltool--4568.org.readthedocs.build/en/4568/recipes/recipe_perfmetrics.html
https://esmvaltool--4568.org.readthedocs.build/en/4568/input.html

---

Before you get started

• ☝ Create an issue to discuss what you are going to do

Checklist

It is the responsibility of the author to make sure the pull request is ready to review. The icons indicate whether the item will be subject to the 🛠 Technical or 🧪 Scientific review.

• 🛠 This pull request has a descriptive title
• 🛠 Code is written according to the code quality guidelines
• 🛠 Documentation is available
• 🛠 The list of authors is up to date
• 🛠 All checks below this pull request were successful

New or updated recipe/diagnostic

• 🧪 Recipe runs successfully
• 🧪 Recipe is well documented
• 🧪 Figure(s) and data look as expected from literature
• 🛠 Provenance information has been added

---

To help with the number of pull requests:

• 🙏 We kindly ask you to review two other open pull requests in this repository

[valeriupredoi]

Thanks @ghossh - just to add that this contains work initially done by @sloosvel with reviews from both @bouweandela and @axel-lauer 🍺

The work was transferred from #4199 that was fully approved by Bouwe and Axel, so that we can transfer ownership from Saskia to Supriyo; Saskia is no longer a developer at ESMValTool, so we think this is the best approach to ensure further continuation of work on this.