diff --git a/doc/recipe/preprocessor.rst b/doc/recipe/preprocessor.rst index 7918f49bf8..a03d7d9502 100644 --- a/doc/recipe/preprocessor.rst +++ b/doc/recipe/preprocessor.rst @@ -803,6 +803,26 @@ or ana4MIPs datasets can be used); in this case the `scheme` is target_grid: ERA-Interim scheme: linear +Regridding on a CORDEX domain grid +---------------------------------- + +It is also possible to regrid to a standard CORDEX domain by using the +CORDEX domain name as ``target_grid``. For example, to regrid to the +``EUR-11`` domain: + +.. code-block:: yaml + + preprocessors: + regrid_preprocessor: + regrid: + target_grid: EUR-11 + scheme: linear + +Any domain name recognized by the ``cordex`` package can be used, for example +``EUR-11``. This creates the target grid from the official CORDEX +domain definition instead of interpreting the value as an ``MxN`` grid +specification. + Regridding on an ``MxN`` grid specification ------------------------------------------- diff --git a/esmvalcore/_recipe/recipe.py b/esmvalcore/_recipe/recipe.py index 8ad4080300..60a29b0672 100644 --- a/esmvalcore/_recipe/recipe.py +++ b/esmvalcore/_recipe/recipe.py @@ -45,6 +45,7 @@ _spec_to_latlonvals, get_cmor_levels, get_reference_levels, + is_cordex_domain, parse_cell_spec, ) from esmvalcore.preprocessor._shared import _group_products @@ -183,7 +184,7 @@ def _update_target_grid( else: # Check that MxN grid spec is correct target_grid = settings["regrid"]["target_grid"] - if isinstance(target_grid, str): + if isinstance(target_grid, str) and not is_cordex_domain(target_grid): parse_cell_spec(target_grid) # Check that cdo spec is correct elif isinstance(target_grid, dict): diff --git a/esmvalcore/preprocessor/_regrid.py b/esmvalcore/preprocessor/_regrid.py index 461e75002a..293eb6eb68 100644 --- a/esmvalcore/preprocessor/_regrid.py +++ b/esmvalcore/preprocessor/_regrid.py @@ -14,17 +14,21 @@ from pathlib import Path from typing import TYPE_CHECKING, Any, Literal +import cordex as cx import dask.array as da import iris import iris.coords +import ncdata.iris_xarray import numpy as np import stratify +import xarray as xr from geopy.geocoders import Nominatim from iris.analysis import ( AreaWeighted, Linear, Nearest, ) +from iris.coord_systems import RotatedGeogCS from iris.cube import Cube from iris.util import broadcast_to_shape @@ -196,6 +200,74 @@ def parse_cell_spec(spec: str) -> tuple[float, float]: return dlon, dlat +def is_cordex_domain(spec: str) -> bool: + """Return ``True`` if ``spec`` is a known CORDEX domain name. + + Parameters + ---------- + spec: + Candidate CORDEX domain identifier (e.g. ``EUR-11``). + + Returns + ------- + bool + Whether ``spec`` is recognised by :mod:`cordex`. + """ + try: + cx.domain_info(spec) + except KeyError: + return False + return True + + +@functools.lru_cache +def _cordex_stock_cube(domain_name: str) -> Cube: + """Create a stock cube for a CORDEX domain target grid. + + Parameters + ---------- + domain_name: + CORDEX domain identifier (e.g. ``EUR-11``). + + Returns + ------- + iris.cube.Cube + Dummy cube with rotated-pole dimension coordinates and geographical + auxiliary coordinates matching the official domain specification. + """ + domain = cx.cordex_domain(domain_name, bounds=True) + domain_info = cx.domain_info(domain_name) + + data = xr.DataArray( + np.zeros((domain.sizes["rlat"], domain.sizes["rlon"]), dtype=np.int32), + dims=["rlat", "rlon"], + coords={ + "rlat": domain["rlat"], + "rlon": domain["rlon"], + "lat": domain["lat"], + "lon": domain["lon"], + }, + name="grid", + ) + (cube,) = ncdata.iris_xarray.cubes_from_xarray(data.to_dataset()) + + coord_system = RotatedGeogCS( + grid_north_pole_latitude=domain_info["pollat"], + grid_north_pole_longitude=domain_info["pollon"], + ) + for dim_coord in ("rlat", "rlon"): + coord = cube.coord(var_name=dim_coord) + coord.coord_system = coord_system + if not coord.has_bounds(): + coord.guess_bounds() + + for aux_coord in ("lat", "lon"): + coord = cube.coord(var_name=aux_coord) + coord.bounds = domain[f"{aux_coord}_vertices"].data + + return cube + + def _generate_cube_from_dimcoords( latdata: np.ndarray | da.Array, londata: np.ndarray | da.Array, @@ -608,19 +680,22 @@ def _get_target_grid_cube( elif isinstance(target_grid, (str, Path)) and os.path.isfile(target_grid): target_grid_cube = iris.load_cube(target_grid) elif isinstance(target_grid, str): - # Generate a target grid from the provided cell-specification - target_grid_cube = _global_stock_cube( - target_grid, - lat_offset, - lon_offset, - ) - # Align the target grid coordinate system to the source - # coordinate system. - src_cs = cube.coord_system() - xcoord = target_grid_cube.coord(axis="x", dim_coords=True) - ycoord = target_grid_cube.coord(axis="y", dim_coords=True) - xcoord.coord_system = src_cs - ycoord.coord_system = src_cs + if is_cordex_domain(target_grid): + target_grid_cube = _cordex_stock_cube(target_grid) + else: + # Generate a target grid from the provided cell-specification + target_grid_cube = _global_stock_cube( + target_grid, + lat_offset, + lon_offset, + ) + # Align the target grid coordinate system to the source + # coordinate system. + src_cs = cube.coord_system() + xcoord = target_grid_cube.coord(axis="x", dim_coords=True) + ycoord = target_grid_cube.coord(axis="y", dim_coords=True) + xcoord.coord_system = src_cs + ycoord.coord_system = src_cs elif isinstance(target_grid, dict): # Generate a target grid from the provided specification, target_grid_cube = _regional_stock_cube(target_grid) diff --git a/tests/unit/preprocessor/_regrid/test_cordex_target_grid.py b/tests/unit/preprocessor/_regrid/test_cordex_target_grid.py new file mode 100644 index 0000000000..d3d0c8fb6e --- /dev/null +++ b/tests/unit/preprocessor/_regrid/test_cordex_target_grid.py @@ -0,0 +1,133 @@ +"""Unit tests for CORDEX domain target grids in regridding.""" + +import cordex as cx +import iris +import iris.coords +import numpy as np +import pytest + +from esmvalcore._recipe import recipe as recipe_module +from esmvalcore.dataset import Dataset +from esmvalcore.preprocessor._regrid import ( + _cordex_stock_cube, + _get_target_grid_cube, + _global_stock_cube, + is_cordex_domain, + parse_cell_spec, +) + + +@pytest.fixture(autouse=True) +def clear_lru_cache(): + """Clear LRU caches for stock cube helpers.""" + yield + _global_stock_cube.cache_clear() + _cordex_stock_cube.cache_clear() + + +@pytest.mark.parametrize( + ("spec", "expected"), + [ + ("EUR-11", True), + ("EUR-44", True), + ("not-a-grid", False), + ("1x1", False), + ("RCA4", False), + ], +) +def test_is_cordex_domain(spec, expected): + """Test CORDEX domain name detection.""" + assert is_cordex_domain(spec) is expected + + +def test_parse_cell_spec_rejects_cordex_domain(): + """CORDEX domains must not be parsed as MxN cell specifications.""" + with pytest.raises(ValueError, match="Invalid MxN cell specification"): + parse_cell_spec("EUR-11") + + +def test_cordex_stock_cube_eur11(): + """Test stock cube for EUR-11 matches the official domain grid.""" + domain = cx.cordex_domain("EUR-11", bounds=True) + cube = _cordex_stock_cube("EUR-11") + + np.testing.assert_array_equal( + cube.coord(var_name="rlat").points, + domain["rlat"].data, + ) + np.testing.assert_array_equal( + cube.coord(var_name="rlon").points, + domain["rlon"].data, + ) + np.testing.assert_array_equal( + cube.coord(var_name="lat").points, + domain["lat"].data, + ) + np.testing.assert_array_equal( + cube.coord(var_name="lon").points, + domain["lon"].data, + ) + assert cube.coord(var_name="rlat").has_bounds() + assert cube.coord(var_name="rlon").has_bounds() + assert cube.coord(var_name="lat").has_bounds() + assert cube.coord(var_name="lon").has_bounds() + + +@pytest.fixture +def global_cube(): + """Return a simple regular global cube for target-grid construction tests.""" + lat_coord = iris.coords.DimCoord( + np.linspace(-85, 85, 18), + standard_name="latitude", + units="degrees", + ) + lon_coord = iris.coords.DimCoord( + np.linspace(5, 355, 36), + standard_name="longitude", + units="degrees", + ) + lat_coord.guess_bounds() + lon_coord.guess_bounds() + return iris.cube.Cube( + np.zeros((18, 36), dtype=np.float32), + dim_coords_and_dims=[(lat_coord, 0), (lon_coord, 1)], + ) + + +def test_get_target_grid_cube_cordex_domain(global_cube): + """Test target grid cube construction for a CORDEX domain.""" + target = _get_target_grid_cube(global_cube, "EUR-11") + assert target.coord(var_name="rlat") is not None + assert target.coord(var_name="rlon") is not None + + +def test_update_target_grid_accepts_cordex_domain(): + """Test recipe preprocessing accepts CORDEX domain target grids.""" + dataset = Dataset( + dataset="RCA4", + project="CORDEX", + domain="EUR-11", + diagnostic="bias", + variable_group="ts", + preprocessor="ts_pp", + ) + settings = {"regrid": {"target_grid": "EUR-11", "scheme": "linear"}} + + recipe_module._update_target_grid(dataset, [dataset], settings) + + assert settings["regrid"]["target_grid"] == "EUR-11" + + +def test_update_target_grid_still_validates_mxn(): + """Test invalid MxN target grids are still rejected.""" + dataset = Dataset( + dataset="RCA4", + project="CORDEX", + diagnostic="bias", + variable_group="ts", + preprocessor="ts_pp", + ) + settings = {"regrid": {"target_grid": "EUR-11x", "scheme": "linear"}} + + with pytest.raises(ValueError, match="Invalid MxN cell specification"): + recipe_module._update_target_grid(dataset, [dataset], settings)