Copy cell measures from the target grid when regridding

esmvalcore/preprocessor/_regrid.py

@@ -986,6 +986,24 @@ def regrid(
     # Rechunk and actually perform the regridding
     cube = _rechunk(cube, target_grid_cube)
     result = regridder(cube)
+    for ancillary_var in cube.ancillary_variables():
+        ancillary_dims = cube.ancillary_variable_dims(ancillary_var)
+        ancillary_slice = tuple(
+            slice(None) if i in ancillary_dims else 0 for i in range(cube.ndim)
+        )
+        ancillary_cube = cube[ancillary_slice].copy(ancillary_var.core_data())
+        ancillary_cube = _rechunk(ancillary_cube, target_grid_cube)
+        ancillary_result = regridder(ancillary_cube)
+        if result.has_lazy_data() and ancillary_result.has_lazy_data():
+            # Keep the chunks of the ancillary variable aligned with the
+            # regridded data variable.
+            ancillary_result.data = ancillary_result.lazy_data().rechunk(
+                result[ancillary_slice].lazy_data().chunks,
+            )
+        result.add_ancillary_variable(
+            ancillary_var.copy(ancillary_result.core_data()),
+            ancillary_dims,
+        )
     # Iris only supports regridding of 1D coordinates and iris-esmf-regrid
     # uses only the DimCoords if both grid_latitude/grid_longitude or
     # projection_x_coordinate/projection_y_coordinate DimCoords and latitude
@@ -1000,6 +1018,8 @@ def regrid(
     # same coordinates when using the regridded cubes as input to the
     # multi-model statistics or similar preprocessor functions later on.
     _update_horizontal_coords(target_grid_cube, result, overwrite=False)
+    _copy_cell_measures(cube, target_grid_cube, result)
+
     return result
 
 
@@ -1156,6 +1176,35 @@ def _update_horizontal_coords(src: Cube, tgt: Cube, overwrite: bool) -> None:
                 tgt.add_aux_coord(coord.copy(), tgt_horizontal_dims)
 
 
+def _copy_cell_measures(
+    cube: Cube,
+    target_grid_cube: Cube,
+    result: Cube,
+) -> None:
+    """Copy cell measures from the target grid to the result cube.
+
+    Copy cell measures from the target grid to the result cube if they are
+    present on the source cube and target grid, to make it look like cell
+    measures are "preserved" during regridding.
+    """
+    for cell_measure in target_grid_cube.cell_measures():
+        if cube.cell_measures(
+            cell_measure.standard_name,
+        ) and target_grid_cube.cell_measures(cell_measure.standard_name):
+            cm_dims = target_grid_cube.cell_measure_dims(cell_measure)
+            cm_slice = tuple(
+                slice(None) if dim in cm_dims else 0
+                for dim in range(target_grid_cube.ndim)
+            )
+            result.add_cell_measure(
+                cell_measure.copy(),
+                tuple(
+                    result.coord_dims(dim_coord)[0]
+                    for dim_coord in target_grid_cube[cm_slice].dim_coords
+                ),
+            )
+
+
 def _create_cube(
     src_cube: Cube,
     data: np.ndarray | da.Array,

tests/integration/preprocessor/_regrid/test_regrid.py

@@ -1,6 +1,11 @@
 """Integration tests for :func:`esmvalcore.preprocessor.regrid`."""
 
+import dask.array as da
 import iris
+import iris.coord_systems
+import iris.coords
+import iris.cube
+import iris.fileformats.pp
 import numpy as np
 import pytest
 from numpy import ma
@@ -364,6 +369,62 @@ def test_regrid__linear_with_mask(self, cache_weights):
         expected[:, 1, 1] = np.array([1.5, 5.5, 9.5])
         assert_array_equal(result.data, expected)
 
+    def test_regrid__linear_with_ancillary(self) -> None:
+        """Test that ancillary coordinates are also regridded."""
+        cube = self.cube.copy()
+        cube.data = cube.lazy_data()
+        cube.add_ancillary_variable(
+            iris.coords.AncillaryVariable(
+                da.arange(2, 6).astype(np.float32).reshape(2, 2),
+                var_name="ancillary",
+            ),
+            (1, 2),
+        )
+        result = regrid(cube, self.grid_for_linear, "linear")
+        ancillary_result = result.ancillary_variable("ancillary")
+        assert isinstance(ancillary_result, iris.coords.AncillaryVariable)
+        assert ancillary_result.has_lazy_data()
+        assert_array_equal(
+            ancillary_result.data,
+            np.array([3.5], dtype=np.float32).reshape(1, 1),
+        )
+
+    @pytest.mark.parametrize("source_has_cell_measure", [True, False])
+    def test_regrid__linear_with_cell_measure(
+        self,
+        source_has_cell_measure: bool,
+    ) -> None:
+        """Test that cell measures are preserved when present on the target grid."""
+        cube = self.cube.copy()
+        if source_has_cell_measure:
+            cube.add_cell_measure(
+                iris.coords.CellMeasure(
+                    np.arange(2, 6).astype(np.float32).reshape(2, 2),
+                    standard_name="cell_area",
+                    units="m2",
+                ),
+                (1, 2),
+            )
+        target_grid = self.grid_for_linear.copy()
+        target_grid.add_cell_measure(
+            iris.coords.CellMeasure(
+                np.ones((1, 1), dtype=np.float32),
+                standard_name="cell_area",
+                units="m2",
+            ),
+            (0, 1),
+        )
+        result = regrid(cube, target_grid, "linear")
+        if source_has_cell_measure:
+            cell_measure_result = result.cell_measure("cell_area")
+            assert isinstance(cell_measure_result, iris.coords.CellMeasure)
+            assert_array_equal(
+                cell_measure_result.data,
+                np.array([1.0], dtype=np.float32).reshape(1, 1),
+            )
+        else:
+            assert not result.cell_measures("cell_area")
+
     @pytest.mark.parametrize("cache_weights", [True, False])
     def test_regrid__nearest(self, cache_weights):
         data = np.empty((1, 1))