-
Notifications
You must be signed in to change notification settings - Fork 48
Add CORDEX domain target grids for regridding #3096
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
base: main
Are you sure you want to change the base?
Changes from all commits
a96e132
8e401c2
ef94822
502047a
0bd0814
1b6e678
d5db352
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change | ||||||
|---|---|---|---|---|---|---|---|---|
|
|
@@ -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. | ||||||||
|
Comment on lines
+823
to
+824
Member
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more.
Suggested change
would suggest leaving that out as it doesn't add much |
||||||||
|
|
||||||||
| Regridding on an ``MxN`` grid specification | ||||||||
| ------------------------------------------- | ||||||||
|
|
||||||||
|
|
||||||||
| Original file line number | Diff line number | Diff line change | ||||||||
|---|---|---|---|---|---|---|---|---|---|---|
|
|
@@ -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`. | ||||||||||
| """ | ||||||||||
|
Member
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more.
Suggested change
Adding this here would avoid the extra nesting in the places where you're calling this function, making that code a bit easier to read. |
||||||||||
| 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) | ||||||||||
|
Member
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Maybe you could simplify this code by using
Member
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This is a bit simpler, but not much. What I meant was something like domain = cordex.domain(domain_name, bounds=True).reset_coords()
# Add a mock variable.
domain["data"] = xr.zeros_like(domain.lon)
# Associate the grid mapping and coordinates with the data.
domain.data.attrs = {
"coordinates": "lat lon",
"grid_mapping": "rotated_latitude_longitude",
}
# Associate the bounds with the coordinates.
domain.lat.attrs["bounds"] = "lat_vertices"
domain.lon.attrs["bounds"] = "lon_vertices"
cube, = ncdata.iris_xarray.cubes_from_xarray(domain)
cube.coord("grid_latitude").guess_bounds()
cube.coord("grid_longitude").guess_bounds() |
||||||||||
| 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) | ||||||||||
|
|
||||||||||
| Original file line number | Diff line number | Diff line change | ||||
|---|---|---|---|---|---|---|
| @@ -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 | ||||||
|
Member
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more.
Suggested change
would test that there is at least one coordinate with the name "rlat". The |
||||||
| assert target.coord(var_name="rlon") is not None | ||||||
|
|
||||||
|
|
||||||
| def test_update_target_grid_accepts_cordex_domain(): | ||||||
|
Member
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Nice that you're adding tests for this! I think this would be a better place for them: ESMValCore/tests/unit/recipe/test_recipe.py Line 809 in ddd67d7
|
||||||
| """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) | ||||||
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It would be nice to add a link to this page so users can find the available domains easily.