Modernize derived coordinates

esmvalcore/cmor/_fixes/cmip6/cesm2.py

@@ -2,13 +2,13 @@
 
 from __future__ import annotations
 
-from shutil import copyfile
 from typing import TYPE_CHECKING
 
 import iris
 import iris.coords
+import ncdata
+import ncdata.netcdf4
 import numpy as np
-from netCDF4 import Dataset
 
 from esmvalcore.cmor._fixes.common import SiconcFixScalarCoord
 from esmvalcore.cmor._fixes.fix import Fix
@@ -19,40 +19,36 @@
     add_scalar_typesea_coord,
     fix_ocean_depth_coord,
 )
+from esmvalcore.iris_helpers import dataset_to_iris
 
 if TYPE_CHECKING:
+    from collections.abc import Sequence
     from pathlib import Path
 
+    from iris.cube import Cube
+
 
 class Cl(Fix):
     """Fixes for ``cl``."""
 
-    def _fix_formula_terms(
-        self,
-        file: str | Path,
-        output_dir: str | Path,
-        add_unique_suffix: bool = False,
-    ) -> Path:
+    @staticmethod
+    def _fix_formula_terms(dataset: ncdata.NcData) -> None:
         """Fix ``formula_terms`` attribute."""
-        new_path = self.get_fixed_filepath(
-            output_dir,
-            file,
-            add_unique_suffix=add_unique_suffix,
+        lev = dataset.variables["lev"]
+        lev.set_attrval("formula_terms", "p0: p0 a: a b: b ps: ps")
+        lev.set_attrval(
+            "standard_name",
+            "atmosphere_hybrid_sigma_pressure_coordinate",
         )
-        copyfile(file, new_path)
-        with Dataset(new_path, mode="a") as dataset:
-            dataset.variables["lev"].formula_terms = "p0: p0 a: a b: b ps: ps"
-            dataset.variables[
-                "lev"
-            ].standard_name = "atmosphere_hybrid_sigma_pressure_coordinate"
-        return new_path
+        lev.set_attrval("units", "1")
+        dataset.variables["lev_bnds"].attributes.pop("units")
 
     def fix_file(
         self,
-        file: str | Path,
-        output_dir: str | Path,
-        add_unique_suffix: bool = False,
-    ) -> Path:
+        file: Path,
+        output_dir: Path,  # noqa: ARG002
+        add_unique_suffix: bool = False,  # noqa: ARG002
+    ) -> Path | Sequence[Cube]:
         """Fix hybrid pressure coordinate.
 
         Adds missing ``formula_terms`` attribute to file.
@@ -79,45 +75,38 @@ def fix_file(
             Path to the fixed file.
 
         """
-        new_path = self._fix_formula_terms(
+        dataset = ncdata.netcdf4.from_nc4(
             file,
-            output_dir,
-            add_unique_suffix=add_unique_suffix,
+            # Use iris-style chunks to avoid mismatching chunks between data
+            # and derived coordinates, as the latter are automatically rechunked
+            # by iris.
+            dim_chunks={
+                "time": "auto",
+                "lev": None,
+                "lat": None,
+                "lon": None,
+                "nbnd": None,
+            },
         )
-        with Dataset(new_path, mode="a") as dataset:
-            dataset.variables["a_bnds"][:] = dataset.variables["a_bnds"][
-                ::-1,
-                :,
-            ]
-            dataset.variables["b_bnds"][:] = dataset.variables["b_bnds"][
-                ::-1,
-                :,
-            ]
-        return new_path
-
-    def fix_metadata(self, cubes):
-        """Fix ``atmosphere_hybrid_sigma_pressure_coordinate``.
-
-        See discussion in #882 for more details on that.
-
-        Parameters
-        ----------
-        cubes : iris.cube.CubeList
-            Input cubes.
-
-        Returns
-        -------
-        iris.cube.CubeList
-
-        """
-        cube = self.get_cube_from_list(cubes)
-        lev_coord = cube.coord(var_name="lev")
-        a_coord = cube.coord(var_name="a")
-        b_coord = cube.coord(var_name="b")
-        lev_coord.points = a_coord.core_points() + b_coord.core_points()
-        lev_coord.bounds = a_coord.core_bounds() + b_coord.core_bounds()
-        lev_coord.units = "1"
-        return cubes
+        self._fix_formula_terms(dataset)
+
+        # Correct order of bounds data
+        a_bnds = dataset.variables["a_bnds"]
+        a_bnds.data = a_bnds.data[::-1, :]
+        b_bnds = dataset.variables["b_bnds"]
+        b_bnds.data = b_bnds.data[::-1, :]
+
+        # Correct lev and lev_bnds data
+        lev = dataset.variables["lev"]
+        lev.data = dataset.variables["a"].data + dataset.variables["b"].data
+        lev_bnds = dataset.variables["lev_bnds"]
+        lev_bnds.data = (
+            dataset.variables["a_bnds"].data + dataset.variables["b_bnds"].data
+        )
+        # Remove 'title' attribute that duplicates long name
+        for var_name in dataset.variables:
+            dataset.variables[var_name].attributes.pop("title", None)
+        return [self.get_cube_from_list(dataset_to_iris(dataset, file))]
 
 
 Cli = Cl

esmvalcore/cmor/_fixes/cmip6/cesm2_waccm.py

@@ -1,8 +1,13 @@
 """Fixes for CESM2-WACCM model."""
 
-from netCDF4 import Dataset
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+import ncdata.netcdf4
 
 from esmvalcore.cmor._fixes.common import SiconcFixScalarCoord
+from esmvalcore.iris_helpers import dataset_to_iris
 
 from .cesm2 import Cl as BaseCl
 from .cesm2 import Fgco2 as BaseFgco2
@@ -12,11 +17,22 @@
 from .cesm2 import Tasmax as BaseTasmax
 from .cesm2 import Tasmin as BaseTasmin
 
+if TYPE_CHECKING:
+    from collections.abc import Sequence
+    from pathlib import Path
+
+    from iris.cube import Cube
+
 
 class Cl(BaseCl):
     """Fixes for cl."""
 
-    def fix_file(self, file, output_dir, add_unique_suffix=False):
+    def fix_file(
+        self,
+        file: Path,
+        output_dir: Path,  # noqa: ARG002
+        add_unique_suffix: bool = False,  # noqa: ARG002
+    ) -> Path | Sequence[Cube]:
         """Fix hybrid pressure coordinate.
 
         Adds missing ``formula_terms`` attribute to file.
@@ -43,21 +59,27 @@ def fix_file(self, file, output_dir, add_unique_suffix=False):
             Path to the fixed file.
 
         """
-        new_path = self._fix_formula_terms(
+        dataset = ncdata.netcdf4.from_nc4(
             file,
-            output_dir,
-            add_unique_suffix=add_unique_suffix,
+            # Use iris-style chunks to avoid mismatching chunks between data
+            # and derived coordinates, as the latter are automatically rechunked
+            # by iris.
+            dim_chunks={
+                "time": "auto",
+                "lev": None,
+                "lat": None,
+                "lon": None,
+                "nbnd": None,
+            },
         )
-        with Dataset(new_path, mode="a") as dataset:
-            dataset.variables["a_bnds"][:] = dataset.variables["a_bnds"][
-                :,
-                ::-1,
-            ]
-            dataset.variables["b_bnds"][:] = dataset.variables["b_bnds"][
-                :,
-                ::-1,
-            ]
-        return new_path
+        self._fix_formula_terms(dataset)
+
+        # Correct order of bounds data
+        a_bnds = dataset.variables["a_bnds"]
+        a_bnds.data = a_bnds.data[:, ::-1]
+        b_bnds = dataset.variables["b_bnds"]
+        b_bnds.data = b_bnds.data[:, ::-1]
+        return [self.get_cube_from_list(dataset_to_iris(dataset, file))]
 
 
 Cli = Cl

esmvalcore/config/__init__.py

@@ -31,6 +31,7 @@
 for attr, value in {
     "save_split_attrs": True,
     "date_microseconds": True,
+    "derived_bounds": True,
 }.items():
     with contextlib.suppress(AttributeError):
         setattr(iris.FUTURE, attr, value)

esmvalcore/preprocessor/_concatenate.py

@@ -13,7 +13,6 @@
 from esmvalcore.cmor.check import CheckLevels
 from esmvalcore.io.esgf.facets import FACETS
 from esmvalcore.iris_helpers import merge_cube_attributes
-from esmvalcore.preprocessor._shared import _rechunk_aux_factory_dependencies
 
 if TYPE_CHECKING:
     from collections.abc import Iterable, Sequence
@@ -282,7 +281,6 @@ def concatenate(
     cubes = _sort_cubes_by_time(cubes)
     _fix_calendars(cubes)
     cubes = _remove_time_overlaps(cubes)
-    cubes = [_rechunk_aux_factory_dependencies(cube) for cube in cubes]
     result = _concatenate_cubes(cubes, check_level=check_level)
 
     if len(result) == 1:

esmvalcore/preprocessor/_regrid.py

@@ -35,7 +35,6 @@
 )
 from esmvalcore.iris_helpers import has_irregular_grid, has_unstructured_grid
 from esmvalcore.preprocessor._shared import (
-    _rechunk_aux_factory_dependencies,
     get_array_module,
     get_dims_along_axes,
     preserve_float_dtype,
@@ -1324,18 +1323,14 @@ def extract_levels(
 
     # Add extra coordinates
     coord_names = [coord.name() for coord in cube.coords()]
-    if coordinate in coord_names:
-        cube = _rechunk_aux_factory_dependencies(cube, coordinate)
-    else:
+    if coordinate not in coord_names:
         # Try to calculate air_pressure from altitude coordinate or
         # vice versa using US standard atmosphere for conversion.
         if coordinate == "air_pressure" and "altitude" in coord_names:
             # Calculate pressure level coordinate from altitude.
-            cube = _rechunk_aux_factory_dependencies(cube, "altitude")
             add_plev_from_altitude(cube)
         if coordinate == "altitude" and "air_pressure" in coord_names:
             # Calculate altitude coordinate from pressure levels.
-            cube = _rechunk_aux_factory_dependencies(cube, "air_pressure")
             add_altitude_from_plev(cube)
 
     src_levels = cube.coord(coordinate)

esmvalcore/preprocessor/_shared.py

@@ -648,39 +648,3 @@ def apply_mask(
 
     array_module = get_array_module(mask, array)
     return array_module.ma.masked_where(mask, array)
-
-
-def _rechunk_aux_factory_dependencies(
-    cube: iris.cube.Cube,
-    coord_name: str | None = None,
-) -> iris.cube.Cube:
-    """Rechunk coordinate aux factory dependencies.
-
-    This ensures that the resulting coordinate has reasonably sized
-    chunks that are aligned with the cube data for optimal computational
-    performance.
-    """
-    # Workaround for https://github.com/SciTools/iris/issues/5457
-    if coord_name is None:
-        factories = cube.aux_factories
-    else:
-        try:
-            factories = [cube.aux_factory(coord_name)]
-        except iris.exceptions.CoordinateNotFoundError:
-            return cube
-
-    cube = cube.copy()
-    cube_chunks = cube.lazy_data().chunks
-    for factory in factories:
-        for orig_coord in factory.dependencies.values():
-            coord_dims = cube.coord_dims(orig_coord)
-            if coord_dims:
-                coord = orig_coord.copy()
-                chunks = tuple(cube_chunks[i] for i in coord_dims)
-                coord.points = coord.lazy_points().rechunk(chunks)
-                if coord.has_bounds():
-                    coord.bounds = coord.lazy_bounds().rechunk(
-                        (*chunks, None),
-                    )
-                cube.replace_coord(coord)
-    return cube

pyproject.toml

@@ -63,7 +63,7 @@ dependencies = [
     "requests",
     "rich",
     "scipy>=1.6",
-    "scitools-iris",
+    "scitools-iris>=3.13",
     "shapely>=2.0.0",
     "stratify>=0.3",
     "xarray",
@@ -223,7 +223,7 @@ geopy = "*"
 humanfriendly = "*"
 intake-esgf = ">=2025.10.22"
 intake-esm = "*"
-iris = "*"
+iris = ">=3.13"
 iris-esmf-regrid = ">=0.11.0"
 iris-grib = ">=0.20.0"
 isodate = ">=0.7.0"