From a3a8752e370f21da77ac255e309328ca15039f12 Mon Sep 17 00:00:00 2001
From: Bouwe Andela <b.andela@esciencecenter.nl>
Date: Wed, 22 Apr 2026 17:50:01 +0200
Subject: [PATCH 1/8] Copy additional coordinates when regridding

---
 esmvalcore/preprocessor/_regrid.py | 79 ++++++++++++++++++++++++------
 1 file changed, 65 insertions(+), 14 deletions(-)

diff --git a/esmvalcore/preprocessor/_regrid.py b/esmvalcore/preprocessor/_regrid.py
index 461e75002a..aa0103f4dc 100644
--- a/esmvalcore/preprocessor/_regrid.py
+++ b/esmvalcore/preprocessor/_regrid.py
@@ -776,6 +776,64 @@ def _get_name_and_shape_key(
     return (name, *shapes)
 
 
+def _copy_horizontal_coords(
+    src: Cube,
+    tgt: Cube,
+    overwrite: bool = False,
+) -> None:
+    """Copy horizontal coordinates from source to target cube.
+
+    Parameters
+    ----------
+    src:
+        Source cube to copy the coordinates from.
+    tgt:
+        Target cube to copy the coordinates to.
+    overwrite:
+        Whether to overwrite existing coordinates on the target cube.
+
+    """
+    src_horizontal_dims = get_dims_along_axes(src, ("X", "Y"))
+    tgt_horizontal_dims = get_dims_along_axes(tgt, ("X", "Y"))
+    for coord_name in [
+        "latitude",
+        "longitude",
+        "grid_latitude",
+        "grid_longitude",
+        "projection_y_coordinate",
+        "projection_x_coordinate",
+    ]:
+        if src.coords(coord_name):
+            if tgt.coords(coord_name):
+                if overwrite:
+                    coord_dims = tgt.coord_dims(coord_name)
+                    tgt.remove_coord(coord_name)
+                else:
+                    continue
+            elif src_horizontal_dims == tgt_horizontal_dims:
+                # Unable to determine coordinate dimensions: assume that the
+                # coordinate dimensions are the same as in the source cube if
+                # the horizontal dimensions of source and target cube match.
+                #
+                # This assumption may be too restrictive when the target cube
+                # has additional non-horizontal dimensions and will fail when the
+                # horizontal dimensions are transposed between source and target
+                # cube, but should work fine in the typical use case where the
+                # source cube (i.e. the target grid) is another, similar dataset.
+                coord_dims = src.coord_dims(coord_name)
+            else:
+                continue
+
+            target_coord = src.coord(coord_name).copy()
+            if (
+                isinstance(target_coord, iris.coords.DimCoord)
+                and len(coord_dims) == 1
+            ):
+                tgt.add_dim_coord(target_coord, coord_dims)
+            else:
+                tgt.add_aux_coord(target_coord, coord_dims)
+
+
 @preserve_float_dtype
 def regrid(
     cube: Cube,
@@ -791,17 +849,16 @@ def regrid(
     Note that the target grid can be a :class:`~iris.cube.Cube`, a
     :class:`~esmvalcore.dataset.Dataset`, a path to a cube
     (:class:`~pathlib.Path` or :obj:`str`), a grid spec (:obj:`str`) in the
-    form of `MxN`, or a :obj:`dict` specifying the target grid.
+    form of `MxN` degrees, or a :obj:`dict` specifying the target grid.
 
     For the latter, the `target_grid` should be a :obj:`dict` with the
     following keys:
 
     - ``start_longitude``: longitude at the center of the first grid cell.
     - ``end_longitude``: longitude at the center of the last grid cell.
-    - ``step_longitude``: constant longitude distance between grid cell
-        centers.
+    - ``step_longitude``: constant longitude distance between grid cell centers.
     - ``start_latitude``: latitude at the center of the first grid cell.
-    - ``end_latitude``: longitude at the center of the last grid cell.
+    - ``end_latitude``: latitude at the center of the last grid cell.
     - ``step_latitude``: constant latitude distance between grid cell centers.
 
     Parameters
@@ -912,15 +969,7 @@ def regrid(
     # -> Return source cube with target coordinates
     if cube.coords("latitude") and cube.coords("longitude"):
         if _horizontal_grid_is_close(cube, target_grid_cube):
-            for coord in ["latitude", "longitude"]:
-                is_dim_coord = cube.coords(coord, dim_coords=True)
-                coord_dims = cube.coord_dims(coord)
-                cube.remove_coord(coord)
-                target_coord = target_grid_cube.coord(coord).copy()
-                if is_dim_coord:
-                    cube.add_dim_coord(target_coord, coord_dims)
-                else:
-                    cube.add_aux_coord(target_coord, coord_dims)
+            _copy_horizontal_coords(target_grid_cube, cube, overwrite=True)
             return cube
 
     # Load scheme and reuse existing regridder if possible
@@ -930,7 +979,9 @@ def regrid(
 
     # Rechunk and actually perform the regridding
     cube = _rechunk(cube, target_grid_cube)
-    return regridder(cube)
+    result = regridder(cube)
+    _copy_horizontal_coords(target_grid_cube, result, overwrite=False)
+    return result
 
 
 def _cache_clear():

From 8d24d14ede7a2f1ed2eefe97338d04e9ec23ed11 Mon Sep 17 00:00:00 2001
From: Bouwe Andela <b.andela@esciencecenter.nl>
Date: Thu, 23 Apr 2026 21:21:16 +0200
Subject: [PATCH 2/8] Improve coord comparison

---
 esmvalcore/preprocessor/_regrid.py | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/esmvalcore/preprocessor/_regrid.py b/esmvalcore/preprocessor/_regrid.py
index aa0103f4dc..0b65a9b7d2 100644
--- a/esmvalcore/preprocessor/_regrid.py
+++ b/esmvalcore/preprocessor/_regrid.py
@@ -1062,10 +1062,16 @@ def _horizontal_grid_is_close(cube1: Cube, cube2: Cube) -> bool:
         coord1 = cube1.coord(coord)
         coord2 = cube2.coord(coord)
 
-        if coord1.shape != coord2.shape:
-            return False
+        if coord1.has_bounds() and coord2.has_bounds():
+            array1 = coord1.core_bounds()
+            array2 = coord2.core_bounds()
+        else:
+            array1 = coord1.core_points()
+            array2 = coord2.core_points()
 
-        if not np.allclose(coord1.bounds, coord2.bounds):
+        if array1.shape != array2.shape:
+            return False
+        if not np.allclose(array1, array2):
             return False
 
     return True

From 2e8ef06844f0a9e0d3fd95c03381d7003a9e1df8 Mon Sep 17 00:00:00 2001
From: Bouwe Andela <b.andela@esciencecenter.nl>
Date: Thu, 11 Jun 2026 12:49:15 +0200
Subject: [PATCH 3/8] Improve copying coordinates and comparing grids

---
 esmvalcore/preprocessor/_regrid.py | 161 +++++++++++++++++------------
 1 file changed, 94 insertions(+), 67 deletions(-)

diff --git a/esmvalcore/preprocessor/_regrid.py b/esmvalcore/preprocessor/_regrid.py
index 0b65a9b7d2..ac9ba9e631 100644
--- a/esmvalcore/preprocessor/_regrid.py
+++ b/esmvalcore/preprocessor/_regrid.py
@@ -17,6 +17,7 @@
 import dask.array as da
 import iris
 import iris.coords
+import iris.exceptions
 import numpy as np
 import stratify
 from geopy.geocoders import Nominatim
@@ -776,64 +777,6 @@ def _get_name_and_shape_key(
     return (name, *shapes)
 
 
-def _copy_horizontal_coords(
-    src: Cube,
-    tgt: Cube,
-    overwrite: bool = False,
-) -> None:
-    """Copy horizontal coordinates from source to target cube.
-
-    Parameters
-    ----------
-    src:
-        Source cube to copy the coordinates from.
-    tgt:
-        Target cube to copy the coordinates to.
-    overwrite:
-        Whether to overwrite existing coordinates on the target cube.
-
-    """
-    src_horizontal_dims = get_dims_along_axes(src, ("X", "Y"))
-    tgt_horizontal_dims = get_dims_along_axes(tgt, ("X", "Y"))
-    for coord_name in [
-        "latitude",
-        "longitude",
-        "grid_latitude",
-        "grid_longitude",
-        "projection_y_coordinate",
-        "projection_x_coordinate",
-    ]:
-        if src.coords(coord_name):
-            if tgt.coords(coord_name):
-                if overwrite:
-                    coord_dims = tgt.coord_dims(coord_name)
-                    tgt.remove_coord(coord_name)
-                else:
-                    continue
-            elif src_horizontal_dims == tgt_horizontal_dims:
-                # Unable to determine coordinate dimensions: assume that the
-                # coordinate dimensions are the same as in the source cube if
-                # the horizontal dimensions of source and target cube match.
-                #
-                # This assumption may be too restrictive when the target cube
-                # has additional non-horizontal dimensions and will fail when the
-                # horizontal dimensions are transposed between source and target
-                # cube, but should work fine in the typical use case where the
-                # source cube (i.e. the target grid) is another, similar dataset.
-                coord_dims = src.coord_dims(coord_name)
-            else:
-                continue
-
-            target_coord = src.coord(coord_name).copy()
-            if (
-                isinstance(target_coord, iris.coords.DimCoord)
-                and len(coord_dims) == 1
-            ):
-                tgt.add_dim_coord(target_coord, coord_dims)
-            else:
-                tgt.add_aux_coord(target_coord, coord_dims)
-
-
 @preserve_float_dtype
 def regrid(
     cube: Cube,
@@ -967,10 +910,9 @@ def regrid(
 
     # Horizontal grids from source and target (almost) match
     # -> Return source cube with target coordinates
-    if cube.coords("latitude") and cube.coords("longitude"):
-        if _horizontal_grid_is_close(cube, target_grid_cube):
-            _copy_horizontal_coords(target_grid_cube, cube, overwrite=True)
-            return cube
+    if _horizontal_grid_is_close(cube, target_grid_cube):
+        _update_horizontal_coords(target_grid_cube, cube)
+        return cube
 
     # Load scheme and reuse existing regridder if possible
     if isinstance(scheme, str):
@@ -980,7 +922,12 @@ def regrid(
     # Rechunk and actually perform the regridding
     cube = _rechunk(cube, target_grid_cube)
     result = regridder(cube)
-    _copy_horizontal_coords(target_grid_cube, result, overwrite=False)
+    # 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
+    # longitude AuxCoords are present, so we need to copy the 2D lat/lon
+    # AuxCoords from the target grid if they are not present on the result.
+    _update_2d_geographical_coords(target_grid_cube, result, overwrite=False)
     return result
 
 
@@ -1057,10 +1004,23 @@ def _horizontal_grid_is_close(cube1: Cube, cube2: Cube) -> bool:
     bool
         ``True`` if grids are close; ``False`` if not.
     """
-    # Go through the 2 expected horizontal coordinates longitude and latitude.
-    for coord in ["latitude", "longitude"]:
-        coord1 = cube1.coord(coord)
-        coord2 = cube2.coord(coord)
+    # Check that both cubes have an X and Y axis so we can compare.
+    if not all(
+        cube1.coords(axis=axis) and cube2.coords(axis=axis)
+        for axis in ("X", "Y")
+    ):
+        return False
+
+    for axis in ("X", "Y"):
+        # DimCoords are kept in-memory, so prefer those for performance.
+        try:
+            coord1 = cube1.coord(axis=axis, dim_coords=True)
+        except iris.exceptions.CoordinateNotFoundError:
+            coord1 = cube1.coord(axis=axis, dim_coords=False)
+        try:
+            coord2 = cube2.coord(axis=axis, dim_coords=True)
+        except iris.exceptions.CoordinateNotFoundError:
+            coord2 = cube2.coord(axis=axis, dim_coords=False)
 
         if coord1.has_bounds() and coord2.has_bounds():
             array1 = coord1.core_bounds()
@@ -1077,6 +1037,73 @@ def _horizontal_grid_is_close(cube1: Cube, cube2: Cube) -> bool:
     return True
 
 
+def _update_2d_geographical_coords(
+    src: Cube,
+    tgt: Cube,
+    overwrite: bool,
+) -> None:
+    """Update 2D latitude and longitude coordinates.
+
+    Parameters
+    ----------
+    src:
+        Source cube to copy the coordinates from.
+    tgt:
+        Target cube to copy the coordinates to.
+    overwrite:
+        If ``True``, overwrite existing coordinates on the target cube. If
+        ``False``, only add the coordinates if they are not already present on the
+        target cube.
+    """
+    tgt_horizontal_dims = get_dims_along_axes(tgt, ("X", "Y"))
+    if len(tgt_horizontal_dims) == 2:
+        for coord_name in [
+            "latitude",
+            "longitude",
+        ]:
+            if (
+                src.coords(coord_name, dim_coords=False)
+                and src.coord(coord_name).ndim == 2
+            ):
+                if overwrite and tgt.coords(coord_name):
+                    tgt.remove_coord(coord_name)
+                if not tgt.coords(coord_name):
+                    tgt.add_aux_coord(
+                        src.coord(coord_name).copy(),
+                        tgt_horizontal_dims,
+                    )
+
+
+def _update_horizontal_coords(src: Cube, tgt: Cube) -> None:
+    """Update horizontal coordinates.
+
+    Parameters
+    ----------
+    src:
+        Source cube to copy the coordinates from.
+    tgt:
+        Target cube to copy the coordinates to.
+
+    """
+    # Copy the horizontal DimCoords.
+    for axis in ("X", "Y"):
+        for new_coord in src.coords(axis=axis, dim_coords=True):
+            if tgt.coords(
+                standard_name=new_coord.standard_name,
+                dim_coords=True,
+            ):
+                old_coord = tgt.coord(
+                    standard_name=new_coord.standard_name,
+                    dim_coords=True,
+                )
+                tgt_dims = tgt.coord_dims(old_coord)
+                tgt.remove_coord(old_coord)
+                tgt.add_dim_coord(new_coord.copy(), tgt_dims)
+
+    # Copy the 2D latitude and longitude coordinates.
+    _update_2d_geographical_coords(src, tgt, overwrite=True)
+
+
 def _create_cube(
     src_cube: Cube,
     data: np.ndarray | da.Array,

From 1616e20e428e8f3e03837a9ec6d99c31565ad429 Mon Sep 17 00:00:00 2001
From: Bouwe Andela <b.andela@esciencecenter.nl>
Date: Thu, 11 Jun 2026 16:55:33 +0200
Subject: [PATCH 4/8] Add some tests for _horizontal_grid_is_close

---
 tests/unit/preprocessor/_regrid/__init__.py   |  7 +-
 .../unit/preprocessor/_regrid/test_regrid.py  | 71 +++++++++++++++++++
 2 files changed, 75 insertions(+), 3 deletions(-)

diff --git a/tests/unit/preprocessor/_regrid/__init__.py b/tests/unit/preprocessor/_regrid/__init__.py
index c9938e051d..211f90ddee 100644
--- a/tests/unit/preprocessor/_regrid/__init__.py
+++ b/tests/unit/preprocessor/_regrid/__init__.py
@@ -3,6 +3,7 @@
 from typing import Literal
 
 import iris
+import iris.coord_systems
 import iris.fileformats
 import numpy as np
 from iris.coords import AuxCoord, CellMethod, DimCoord
@@ -77,9 +78,9 @@ def _make_cube(  # noqa: PLR0915,C901
         cube.add_dim_coord(_make_vcoord(z, dtype=dtype), 0)
 
     if grid == "rotated":
-        # Create a synthetic test latitude coordinate.
+        # Create a synthetic test grid_latitude coordinate.
         data = np.arange(y, dtype=dtype) + 1
-        cs = iris.coord_systems.GeogCS(iris.fileformats.pp.EARTH_RADIUS)
+        cs = iris.coord_systems.RotatedGeogCS(6.08, -166.92)
         kwargs = {
             "standard_name": "grid_latitude",
             "long_name": "latitude in rotated pole grid",
@@ -93,7 +94,7 @@ def _make_cube(  # noqa: PLR0915,C901
             ycoord.guess_bounds()
         cube.add_dim_coord(ycoord, 1)
 
-        # Create a synthetic test longitude coordinate.
+        # Create a synthetic test grid_longitude coordinate.
         data = np.arange(x, dtype=dtype) + 1
         kwargs = {
             "standard_name": "grid_longitude",
diff --git a/tests/unit/preprocessor/_regrid/test_regrid.py b/tests/unit/preprocessor/_regrid/test_regrid.py
index 7b41305c7e..bbfd0f2382 100644
--- a/tests/unit/preprocessor/_regrid/test_regrid.py
+++ b/tests/unit/preprocessor/_regrid/test_regrid.py
@@ -3,6 +3,9 @@
 import dask
 import dask.array as da
 import iris
+import iris.coord_systems
+import iris.coords
+import iris.cube
 import numpy as np
 import pytest
 
@@ -238,6 +241,74 @@ def test_horizontal_grid_is_close(cube2_spec: dict, expected: bool) -> None:
     assert _horizontal_grid_is_close(cube1, cube2) == expected
 
 
+def test_horizontal_grid_is_close_aux_coords():
+    """Test for `_horizontal_grid_is_close`."""
+    cube1 = _make_cube(lat=LAT_SPEC1, lon=LON_SPEC1)
+    cube2 = _make_cube(lat=LAT_SPEC1, lon=LON_SPEC1)
+    lat = cube1.coord("latitude")
+    lon = cube1.coord("longitude")
+    cube1.remove_coord(lat)
+    cube1.remove_coord(lon)
+    cube1.add_aux_coord(lat, 0)
+    cube1.add_aux_coord(lon, 1)
+
+    lat = cube2.coord("latitude")
+    lon = cube2.coord("longitude")
+    cube2.remove_coord(lat)
+    cube2.remove_coord(lon)
+    cube2.add_aux_coord(lat, 0)
+    cube2.add_aux_coord(lon, 1)
+
+    assert _horizontal_grid_is_close(cube1, cube2) is True
+
+
+@pytest.mark.parametrize("close", [True, False])
+def test_horizontal_grid_is_close_rotated(close):
+    """Test for `_horizontal_grid_is_close`."""
+    cs = iris.coord_systems.RotatedGeogCS(6.08, -166.92)
+    cube1 = iris.cube.Cube(
+        np.zeros((2, 3)),
+        dim_coords_and_dims=[
+            (
+                iris.coords.DimCoord(
+                    [1, 2],
+                    standard_name="grid_latitude",
+                    units="degrees",
+                    coord_system=cs,
+                ),
+                0,
+            ),
+            (
+                iris.coords.DimCoord(
+                    [1, 2, 3],
+                    standard_name="grid_longitude",
+                    units="degrees",
+                    coord_system=cs,
+                ),
+                1,
+            ),
+        ],
+    )
+    cube2 = cube1.copy()
+    if not close:
+        cube2.coord("grid_latitude").points = (
+            cube2.coord("grid_latitude").points + 1.0
+        )
+
+    assert _horizontal_grid_is_close(cube1, cube2) is close
+
+
+def test_horizontal_grid_is_close_no_coords():
+    """Test for `_horizontal_grid_is_close`."""
+    cube1 = _make_cube(lat=LAT_SPEC1, lon=LON_SPEC1)
+    cube2 = _make_cube(lat=LAT_SPEC1, lon=LON_SPEC1)
+    cube1.remove_coord("latitude")
+    cube2.remove_coord("latitude")
+    cube2.remove_coord("longitude")
+
+    assert _horizontal_grid_is_close(cube1, cube2) is False
+
+
 def test_regrid_is_skipped_if_grids_are_the_same_dim_coord(mocker):
     """Test that regridding is skipped if the grids are the same."""
     mock_get_regridder = mocker.patch(

From 8ee06358ff77927d18bffa4bf207ce894ceac805 Mon Sep 17 00:00:00 2001
From: Bouwe Andela <b.andela@esciencecenter.nl>
Date: Fri, 12 Jun 2026 13:35:55 +0200
Subject: [PATCH 5/8] Identify coordinates by dimension

---
 esmvalcore/preprocessor/_regrid.py            | 94 ++++++++-----------
 .../preprocessor/_regrid/test_regrid.py       | 42 +++++++++
 2 files changed, 80 insertions(+), 56 deletions(-)

diff --git a/esmvalcore/preprocessor/_regrid.py b/esmvalcore/preprocessor/_regrid.py
index ac9ba9e631..ff7de49873 100644
--- a/esmvalcore/preprocessor/_regrid.py
+++ b/esmvalcore/preprocessor/_regrid.py
@@ -911,7 +911,7 @@ def regrid(
     # Horizontal grids from source and target (almost) match
     # -> Return source cube with target coordinates
     if _horizontal_grid_is_close(cube, target_grid_cube):
-        _update_horizontal_coords(target_grid_cube, cube)
+        _update_horizontal_coords(target_grid_cube, cube, overwrite=True)
         return cube
 
     # Load scheme and reuse existing regridder if possible
@@ -927,7 +927,15 @@ def regrid(
     # projection_x_coordinate/projection_y_coordinate DimCoords and latitude
     # longitude AuxCoords are present, so we need to copy the 2D lat/lon
     # AuxCoords from the target grid if they are not present on the result.
-    _update_2d_geographical_coords(target_grid_cube, result, overwrite=False)
+    #
+    # Similarly, when the target grid is a 2D latitude/longitude grid,
+    # regridding does not use the DimCoord because it has no meaning. This
+    # is typical in CMIP6 ocean data (e.g. coordinates named i/j, cell index
+    # along first/second dimension), so we need to copy the dummy DimCoords
+    # from the target grid in this case to ensure the resulting cube has the
+    # 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)
     return result
 
 
@@ -1037,12 +1045,8 @@ def _horizontal_grid_is_close(cube1: Cube, cube2: Cube) -> bool:
     return True
 
 
-def _update_2d_geographical_coords(
-    src: Cube,
-    tgt: Cube,
-    overwrite: bool,
-) -> None:
-    """Update 2D latitude and longitude coordinates.
+def _update_horizontal_coords(src: Cube, tgt: Cube, overwrite: bool) -> None:
+    """Update horizontal coordinates.
 
     Parameters
     ----------
@@ -1052,56 +1056,34 @@ def _update_2d_geographical_coords(
         Target cube to copy the coordinates to.
     overwrite:
         If ``True``, overwrite existing coordinates on the target cube. If
-        ``False``, only add the coordinates if they are not already present on the
-        target cube.
+        ``False``, only add the coordinates if they are not already present on
+        the target cube.
     """
+    src_horizontal_dims = get_dims_along_axes(src, ("X", "Y"))
     tgt_horizontal_dims = get_dims_along_axes(tgt, ("X", "Y"))
-    if len(tgt_horizontal_dims) == 2:
-        for coord_name in [
-            "latitude",
-            "longitude",
-        ]:
-            if (
-                src.coords(coord_name, dim_coords=False)
-                and src.coord(coord_name).ndim == 2
-            ):
-                if overwrite and tgt.coords(coord_name):
-                    tgt.remove_coord(coord_name)
-                if not tgt.coords(coord_name):
-                    tgt.add_aux_coord(
-                        src.coord(coord_name).copy(),
-                        tgt_horizontal_dims,
-                    )
-
-
-def _update_horizontal_coords(src: Cube, tgt: Cube) -> None:
-    """Update horizontal coordinates.
-
-    Parameters
-    ----------
-    src:
-        Source cube to copy the coordinates from.
-    tgt:
-        Target cube to copy the coordinates to.
-
-    """
-    # Copy the horizontal DimCoords.
-    for axis in ("X", "Y"):
-        for new_coord in src.coords(axis=axis, dim_coords=True):
-            if tgt.coords(
-                standard_name=new_coord.standard_name,
-                dim_coords=True,
-            ):
-                old_coord = tgt.coord(
-                    standard_name=new_coord.standard_name,
-                    dim_coords=True,
-                )
-                tgt_dims = tgt.coord_dims(old_coord)
-                tgt.remove_coord(old_coord)
-                tgt.add_dim_coord(new_coord.copy(), tgt_dims)
-
-    # Copy the 2D latitude and longitude coordinates.
-    _update_2d_geographical_coords(src, tgt, overwrite=True)
+    if len(src_horizontal_dims) != len(tgt_horizontal_dims):
+        return
+
+    # Copy DimCoords. This assumes they are in the same order in the source
+    # and target cube, which seems to be required for the regridders.
+    for src_dim, tgt_dim in zip(
+        src_horizontal_dims,
+        tgt_horizontal_dims,
+        strict=True,
+    ):
+        for coord in src.coords(dimensions=src_dim, dim_coords=True):
+            if overwrite and tgt.coords(coord.name()):
+                tgt.remove_coord(coord.name())
+            if not tgt.coords(coord.name()):
+                tgt.add_dim_coord(coord.copy(), tgt_dim)
+
+    # Copy 2D latitude and longitude coordinates.
+    if len(src_horizontal_dims) == 2 and len(tgt_horizontal_dims) == 2:
+        for coord in src.coords(dimensions=src_horizontal_dims):
+            if overwrite and tgt.coords(coord.name()):
+                tgt.remove_coord(coord.name())
+            if not tgt.coords(coord.name()):
+                tgt.add_aux_coord(coord.copy(), tgt_horizontal_dims)
 
 
 def _create_cube(
diff --git a/tests/integration/preprocessor/_regrid/test_regrid.py b/tests/integration/preprocessor/_regrid/test_regrid.py
index e374262619..7877b0bdcf 100644
--- a/tests/integration/preprocessor/_regrid/test_regrid.py
+++ b/tests/integration/preprocessor/_regrid/test_regrid.py
@@ -166,6 +166,48 @@ def test_regrid__linear_multidim(self):
         result.data = np.round(result.data, 1)
         assert_array_equal(result.data, expected)
 
+    @pytest.mark.parametrize("close", [True, False])
+    def test_regrid__nearest_cordex_result_has_all_coords(self, close):
+        """Test that the result of regridding has both 1D and 2D lat and lon."""
+        target_grid = self.multidim_cube.copy()
+        if not close:
+            target_grid.coord("grid_longitude").points = (
+                target_grid.coord("grid_longitude").points + 0.1
+            )
+            target_grid.coord("longitude").points = (
+                target_grid.coord("longitude").points + 0.1
+            )
+        result = regrid(
+            self.multidim_cube,
+            target_grid,
+            "nearest",
+        )
+        expected = np.ma.masked_array(
+            [
+                [
+                    [3.0, 3.0],
+                    [3.0, 3.0],
+                ],
+                [
+                    [7.0, 7.0],
+                    [7.0, 7.0],
+                ],
+                [
+                    [11.0, 11.0],
+                    [11.0, 11.0],
+                ],
+            ],
+            mask=False,
+        )
+        for coord in (
+            "latitude",
+            "longitude",
+            "grid_latitude",
+            "grid_longitude",
+        ):
+            assert result.coord(coord) == target_grid.coord(coord)
+        assert_array_equal(result.data, expected)
+
     @pytest.mark.parametrize("cache_weights", [True, False])
     def test_regrid__linear_file(self, tmp_path, cache_weights):
         file = tmp_path / "file.nc"

From 2f205d9026d12d3650e4fd672dac293f543c9f91 Mon Sep 17 00:00:00 2001
From: Bouwe Andela <b.andela@esciencecenter.nl>
Date: Fri, 12 Jun 2026 14:39:47 +0200
Subject: [PATCH 6/8] Better coordinate matching and tests

---
 esmvalcore/preprocessor/_regrid.py            | 24 ++++----
 .../preprocessor/_regrid/test_regrid.py       | 60 ++++++++++---------
 .../unit/preprocessor/_regrid/test_regrid.py  |  2 +-
 3 files changed, 44 insertions(+), 42 deletions(-)

diff --git a/esmvalcore/preprocessor/_regrid.py b/esmvalcore/preprocessor/_regrid.py
index ff7de49873..6628e6e48d 100644
--- a/esmvalcore/preprocessor/_regrid.py
+++ b/esmvalcore/preprocessor/_regrid.py
@@ -1012,23 +1012,23 @@ def _horizontal_grid_is_close(cube1: Cube, cube2: Cube) -> bool:
     bool
         ``True`` if grids are close; ``False`` if not.
     """
-    # Check that both cubes have an X and Y axis so we can compare.
-    if not all(
-        cube1.coords(axis=axis) and cube2.coords(axis=axis)
-        for axis in ("X", "Y")
-    ):
-        return False
-
     for axis in ("X", "Y"):
         # DimCoords are kept in-memory, so prefer those for performance.
         try:
             coord1 = cube1.coord(axis=axis, dim_coords=True)
         except iris.exceptions.CoordinateNotFoundError:
-            coord1 = cube1.coord(axis=axis, dim_coords=False)
-        try:
-            coord2 = cube2.coord(axis=axis, dim_coords=True)
-        except iris.exceptions.CoordinateNotFoundError:
-            coord2 = cube2.coord(axis=axis, dim_coords=False)
+            try:
+                coord1 = cube1.coord(axis=axis, dim_coords=False)
+            except iris.exceptions.CoordinateNotFoundError:
+                # No horizontal coordinate found.
+                return False
+
+        for coord2 in cube2.coords(axis=axis):
+            if coord2.shape == coord1.shape:
+                break
+        else:
+            # No matching coordinate found.
+            return False
 
         if coord1.has_bounds() and coord2.has_bounds():
             array1 = coord1.core_bounds()
diff --git a/tests/integration/preprocessor/_regrid/test_regrid.py b/tests/integration/preprocessor/_regrid/test_regrid.py
index 7877b0bdcf..6f0e772d55 100644
--- a/tests/integration/preprocessor/_regrid/test_regrid.py
+++ b/tests/integration/preprocessor/_regrid/test_regrid.py
@@ -42,10 +42,15 @@ def setUp(self):
         )
 
         # Setup cube with multiple horizontal dimensions
-        self.multidim_cube = _make_cube(data, grid="rotated", aux_coord=False)
+        self.multidim_cube = _make_cube(
+            data,
+            grid="rotated",
+            dtype=np.float64,
+            aux_coord=False,
+        )
         lats, lons = np.meshgrid(
-            np.arange(1, data.shape[-2] + 1),
-            np.arange(1, data.shape[-1] + 1),
+            np.arange(1, data.shape[-2] + 1).astype(np.float64),
+            np.arange(1, data.shape[-1] + 1).astype(np.float64),
         )
         self.multidim_cube.add_aux_coord(
             iris.coords.AuxCoord(
@@ -166,38 +171,36 @@ def test_regrid__linear_multidim(self):
         result.data = np.round(result.data, 1)
         assert_array_equal(result.data, expected)
 
+    @pytest.mark.parametrize(
+        "use_src_coords",
+        [
+            ["latitude", "longitude"],
+            ["grid_latitude", "grid_longitude"],
+        ],
+    )
     @pytest.mark.parametrize("close", [True, False])
-    def test_regrid__nearest_cordex_result_has_all_coords(self, close):
+    def test_regrid__nearest_cordex_result_has_all_coords(
+        self,
+        use_src_coords,
+        close,
+    ):
         """Test that the result of regridding has both 1D and 2D lat and lon."""
         target_grid = self.multidim_cube.copy()
-        if not close:
-            target_grid.coord("grid_longitude").points = (
-                target_grid.coord("grid_longitude").points + 0.1
-            )
-            target_grid.coord("longitude").points = (
-                target_grid.coord("longitude").points + 0.1
-            )
+        offset = 1e-9 if close else 0.1
+        target_grid.coord("grid_longitude").points = (
+            target_grid.coord("grid_longitude").points + offset
+        )
+        target_grid.coord("grid_longitude").bounds = (
+            target_grid.coord("grid_longitude").bounds + offset
+        )
+        target_grid.coord("longitude").points = (
+            target_grid.coord("longitude").points + offset
+        )
         result = regrid(
             self.multidim_cube,
             target_grid,
             "nearest",
-        )
-        expected = np.ma.masked_array(
-            [
-                [
-                    [3.0, 3.0],
-                    [3.0, 3.0],
-                ],
-                [
-                    [7.0, 7.0],
-                    [7.0, 7.0],
-                ],
-                [
-                    [11.0, 11.0],
-                    [11.0, 11.0],
-                ],
-            ],
-            mask=False,
+            use_src_coords=use_src_coords,
         )
         for coord in (
             "latitude",
@@ -206,7 +209,6 @@ def test_regrid__nearest_cordex_result_has_all_coords(self, close):
             "grid_longitude",
         ):
             assert result.coord(coord) == target_grid.coord(coord)
-        assert_array_equal(result.data, expected)
 
     @pytest.mark.parametrize("cache_weights", [True, False])
     def test_regrid__linear_file(self, tmp_path, cache_weights):
diff --git a/tests/unit/preprocessor/_regrid/test_regrid.py b/tests/unit/preprocessor/_regrid/test_regrid.py
index bbfd0f2382..b2f1826dbb 100644
--- a/tests/unit/preprocessor/_regrid/test_regrid.py
+++ b/tests/unit/preprocessor/_regrid/test_regrid.py
@@ -39,7 +39,7 @@ def _make_coord(
     coord = iris.coords.DimCoord(
         np.linspace(start, stop, step),
         standard_name=name,
-        units="degrees",
+        units="degrees_north" if name == "latitude" else "degrees_east",
     )
     coord.guess_bounds()
     return coord

From 93725aa39881cf5200e6beb5254f27a68dbda59d Mon Sep 17 00:00:00 2001
From: Bouwe Andela <b.andela@esciencecenter.nl>
Date: Fri, 12 Jun 2026 14:48:42 +0200
Subject: [PATCH 7/8] Update test

---
 tests/unit/preprocessor/_regrid/test_regrid.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/tests/unit/preprocessor/_regrid/test_regrid.py b/tests/unit/preprocessor/_regrid/test_regrid.py
index b2f1826dbb..9e0610b0db 100644
--- a/tests/unit/preprocessor/_regrid/test_regrid.py
+++ b/tests/unit/preprocessor/_regrid/test_regrid.py
@@ -303,8 +303,6 @@ def test_horizontal_grid_is_close_no_coords():
     cube1 = _make_cube(lat=LAT_SPEC1, lon=LON_SPEC1)
     cube2 = _make_cube(lat=LAT_SPEC1, lon=LON_SPEC1)
     cube1.remove_coord("latitude")
-    cube2.remove_coord("latitude")
-    cube2.remove_coord("longitude")
 
     assert _horizontal_grid_is_close(cube1, cube2) is False
 

From b6a1e2358d49d1a0f443695d116a83b97fbcede7 Mon Sep 17 00:00:00 2001
From: Bouwe Andela <b.andela@esciencecenter.nl>
Date: Fri, 12 Jun 2026 15:09:26 +0200
Subject: [PATCH 8/8] Remove unnessary check

---
 esmvalcore/preprocessor/_regrid.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/esmvalcore/preprocessor/_regrid.py b/esmvalcore/preprocessor/_regrid.py
index 6628e6e48d..4a037e0c4c 100644
--- a/esmvalcore/preprocessor/_regrid.py
+++ b/esmvalcore/preprocessor/_regrid.py
@@ -1037,8 +1037,6 @@ def _horizontal_grid_is_close(cube1: Cube, cube2: Cube) -> bool:
             array1 = coord1.core_points()
             array2 = coord2.core_points()
 
-        if array1.shape != array2.shape:
-            return False
         if not np.allclose(array1, array2):
             return False