From 0126dc5c54e09f53cdb41568e519e4ff45757019 Mon Sep 17 00:00:00 2001
From: "stephen.worsley" <stephen.worsley@metoffice.gov.uk>
Date: Mon, 29 Sep 2025 22:07:24 +0100
Subject: [PATCH 01/14] add 2D guess bounds

---
 lib/iris/analysis/_grid_angles.py             | 58 +++++++++++++++++++
 .../cartography/test_gridcell_angles.py       | 11 ++++
 2 files changed, 69 insertions(+)

diff --git a/lib/iris/analysis/_grid_angles.py b/lib/iris/analysis/_grid_angles.py
index 3ba406e02a..ca8fbd2fb7 100644
--- a/lib/iris/analysis/_grid_angles.py
+++ b/lib/iris/analysis/_grid_angles.py
@@ -453,3 +453,61 @@ def rotate_grid_vectors(u_cube, v_cube, grid_angles_cube=None, grid_angles_kwarg
     v_out.data = np.ma.masked_array(vv, mask=mask)
 
     return u_out, v_out
+
+
+def _vectorised_matmul(mats, vecs):
+    return np.einsum("ijk,ji->ki", mats, vecs)
+
+def _generate_180_mats_from_uvecs(uvecs):
+    mats = np.einsum("ji,ki->ijk", uvecs, uvecs) * 2
+    np.einsum("ijj->ij", mats)[:] -= 1
+    return mats
+
+def _2D_guess_bounds_first_pass(array):
+    # average and normalise, boundary buffer represents edges and corners
+    result_array = np.zeros((array.shape[0], array.shape[1] + 1, array.shape[2] + 1))
+    pads = ((0, 1), (1, 0))
+    for pad_i in pads:
+        for pad_j in pads:
+            result_array += np.pad(array, ((0, 0), pad_i, pad_j))
+
+    # normalise
+    result_array /= np.linalg.norm(result_array, ord=2, axis=0)[np.newaxis, ...]
+    return result_array
+
+def _2D_gb_buffer_outer(array_shape):
+    # return appropriate numpy slice for outer halo
+    _, x, y = array_shape
+    x_i = list(range(x)) + ([x - 1] * (y - 2)) + list(range(x))[::-1] + ([0] * (y - 2))
+    y_i = ([0] * (x - 1)) + list(range(y)) + ([y - 1] * (x - 2)) + list(range(1, y))[::-1]
+    return np.s_[:, x_i, y_i]
+
+def _2D_gb_buffer_inner(array_shape):
+    # return appropriate numpy slice for inner halo
+    _, x, y = array_shape
+    x_i = [1] + list(range(1, x - 1)) + ([x - 2] * y) + list(range(1, x - 1))[::-1] + ([1] * (y - 1))
+    y_i = ([1] * x) + list(range(1, y - 1)) + ([y - 1] * x) + list(range(1, y - 1))[::-1]
+    return np.s_[:, x_i, y_i]
+
+def _2D_geuss_bounds(cube):
+    lons = cube.coord(axis="X")
+    lats = cube.coord(axis="Y")
+
+    # TODO: check units and coordsystem and dims and such
+
+    lon_array = lons.points
+    lat_array = lats.points
+    xyz_array = _3d_xyz_from_latlon(lon_array, lat_array)
+
+    result_xyz = _2D_guess_bounds_first_pass(xyz_array)
+    outer_inds = _2D_gb_buffer_outer(xyz_array.shape)
+    inner_inds = _2D_gb_buffer_inner(xyz_array.shape)
+
+    mats = _generate_180_mats_from_uvecs(result_xyz[outer_inds])
+    result_xyz[outer_inds] = _vectorised_matmul(mats, result_xyz[inner_inds])
+
+    result_lon_bounds, result_lat_bounds = _latlon_from_xyz(result_xyz)
+
+    # add these bounds cf style
+    lons.bounds = np.stack([result_lon_bounds[:-1,:-1], result_lon_bounds[:-1,1:], result_lon_bounds[1:,1:], result_lon_bounds[1:,:-1]], axis=2)
+    lats.bounds = np.stack([result_lat_bounds[:-1,:-1], result_lat_bounds[:-1,1:], result_lat_bounds[1:,1:], result_lat_bounds[1:,:-1]], axis=2)
diff --git a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
index 32af84317c..41956e1410 100644
--- a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
+++ b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
@@ -12,6 +12,7 @@
 import pytest
 
 from iris.analysis.cartography import gridcell_angles
+from iris.analysis._grid_angles import _2D_geuss_bounds
 from iris.coords import AuxCoord
 from iris.cube import Cube
 from iris.tests import _shared_utils
@@ -301,3 +302,13 @@ def test_fail_noncoord_coord(self):
     def test_fail_bad_method(self):
         with pytest.raises(ValueError, match="unrecognised cell_angle_boundpoints"):
             self._check_multiple_orientations_and_latitudes(method="something_unknown")
+
+def test_2D_guess_bounds():
+    cube = _2d_multicells_testcube()
+    assert not cube.coord("latitude").is_contiguous()
+    assert not cube.coord("longitude").is_contiguous()
+
+    _2D_geuss_bounds(cube)
+    assert cube.coord("latitude").is_contiguous()
+    assert cube.coord("longitude").is_contiguous()
+

From 070006f7449612876aa2ead1f12bf7bc7860ddfe Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Thu, 13 Nov 2025 15:27:21 +0000
Subject: [PATCH 02/14] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 lib/iris/analysis/_grid_angles.py             | 41 ++++++++++++++++---
 .../cartography/test_gridcell_angles.py       |  4 +-
 2 files changed, 38 insertions(+), 7 deletions(-)

diff --git a/lib/iris/analysis/_grid_angles.py b/lib/iris/analysis/_grid_angles.py
index ca8fbd2fb7..7a579df185 100644
--- a/lib/iris/analysis/_grid_angles.py
+++ b/lib/iris/analysis/_grid_angles.py
@@ -458,11 +458,13 @@ def rotate_grid_vectors(u_cube, v_cube, grid_angles_cube=None, grid_angles_kwarg
 def _vectorised_matmul(mats, vecs):
     return np.einsum("ijk,ji->ki", mats, vecs)
 
+
 def _generate_180_mats_from_uvecs(uvecs):
     mats = np.einsum("ji,ki->ijk", uvecs, uvecs) * 2
     np.einsum("ijj->ij", mats)[:] -= 1
     return mats
 
+
 def _2D_guess_bounds_first_pass(array):
     # average and normalise, boundary buffer represents edges and corners
     result_array = np.zeros((array.shape[0], array.shape[1] + 1, array.shape[2] + 1))
@@ -475,20 +477,33 @@ def _2D_guess_bounds_first_pass(array):
     result_array /= np.linalg.norm(result_array, ord=2, axis=0)[np.newaxis, ...]
     return result_array
 
+
 def _2D_gb_buffer_outer(array_shape):
     # return appropriate numpy slice for outer halo
     _, x, y = array_shape
     x_i = list(range(x)) + ([x - 1] * (y - 2)) + list(range(x))[::-1] + ([0] * (y - 2))
-    y_i = ([0] * (x - 1)) + list(range(y)) + ([y - 1] * (x - 2)) + list(range(1, y))[::-1]
+    y_i = (
+        ([0] * (x - 1)) + list(range(y)) + ([y - 1] * (x - 2)) + list(range(1, y))[::-1]
+    )
     return np.s_[:, x_i, y_i]
 
+
 def _2D_gb_buffer_inner(array_shape):
     # return appropriate numpy slice for inner halo
     _, x, y = array_shape
-    x_i = [1] + list(range(1, x - 1)) + ([x - 2] * y) + list(range(1, x - 1))[::-1] + ([1] * (y - 1))
-    y_i = ([1] * x) + list(range(1, y - 1)) + ([y - 1] * x) + list(range(1, y - 1))[::-1]
+    x_i = (
+        [1]
+        + list(range(1, x - 1))
+        + ([x - 2] * y)
+        + list(range(1, x - 1))[::-1]
+        + ([1] * (y - 1))
+    )
+    y_i = (
+        ([1] * x) + list(range(1, y - 1)) + ([y - 1] * x) + list(range(1, y - 1))[::-1]
+    )
     return np.s_[:, x_i, y_i]
 
+
 def _2D_geuss_bounds(cube):
     lons = cube.coord(axis="X")
     lats = cube.coord(axis="Y")
@@ -509,5 +524,21 @@ def _2D_geuss_bounds(cube):
     result_lon_bounds, result_lat_bounds = _latlon_from_xyz(result_xyz)
 
     # add these bounds cf style
-    lons.bounds = np.stack([result_lon_bounds[:-1,:-1], result_lon_bounds[:-1,1:], result_lon_bounds[1:,1:], result_lon_bounds[1:,:-1]], axis=2)
-    lats.bounds = np.stack([result_lat_bounds[:-1,:-1], result_lat_bounds[:-1,1:], result_lat_bounds[1:,1:], result_lat_bounds[1:,:-1]], axis=2)
+    lons.bounds = np.stack(
+        [
+            result_lon_bounds[:-1, :-1],
+            result_lon_bounds[:-1, 1:],
+            result_lon_bounds[1:, 1:],
+            result_lon_bounds[1:, :-1],
+        ],
+        axis=2,
+    )
+    lats.bounds = np.stack(
+        [
+            result_lat_bounds[:-1, :-1],
+            result_lat_bounds[:-1, 1:],
+            result_lat_bounds[1:, 1:],
+            result_lat_bounds[1:, :-1],
+        ],
+        axis=2,
+    )
diff --git a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
index 41956e1410..bc61b39a67 100644
--- a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
+++ b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
@@ -11,8 +11,8 @@
 import numpy as np
 import pytest
 
-from iris.analysis.cartography import gridcell_angles
 from iris.analysis._grid_angles import _2D_geuss_bounds
+from iris.analysis.cartography import gridcell_angles
 from iris.coords import AuxCoord
 from iris.cube import Cube
 from iris.tests import _shared_utils
@@ -303,6 +303,7 @@ def test_fail_bad_method(self):
         with pytest.raises(ValueError, match="unrecognised cell_angle_boundpoints"):
             self._check_multiple_orientations_and_latitudes(method="something_unknown")
 
+
 def test_2D_guess_bounds():
     cube = _2d_multicells_testcube()
     assert not cube.coord("latitude").is_contiguous()
@@ -311,4 +312,3 @@ def test_2D_guess_bounds():
     _2D_geuss_bounds(cube)
     assert cube.coord("latitude").is_contiguous()
     assert cube.coord("longitude").is_contiguous()
-

From aa16d3f57bbea0055af3be951e53f671df89fd34 Mon Sep 17 00:00:00 2001
From: "stephen.worsley" <stephen.worsley@metoffice.gov.uk>
Date: Wed, 26 Nov 2025 08:38:48 +0000
Subject: [PATCH 03/14] fix _2D_guess_bounds, add reflection method

---
 lib/iris/analysis/_grid_angles.py             | 48 +++++++++++++++++--
 .../cartography/test_gridcell_angles.py       |  4 ++
 2 files changed, 47 insertions(+), 5 deletions(-)

diff --git a/lib/iris/analysis/_grid_angles.py b/lib/iris/analysis/_grid_angles.py
index 7a579df185..1177fa1850 100644
--- a/lib/iris/analysis/_grid_angles.py
+++ b/lib/iris/analysis/_grid_angles.py
@@ -464,6 +464,42 @@ def _generate_180_mats_from_uvecs(uvecs):
     np.einsum("ijj->ij", mats)[:] -= 1
     return mats
 
+def _generate_reflection_mats_from_uvecs(xyz_array):
+    vecs_0 = xyz_array[:,:,0]
+    vecs_1 = xyz_array[:,-1,:]
+    vecs_2 = xyz_array[:,:,-1]
+    vecs_3 = xyz_array[:,0,:]
+
+    def make_reflection_mat(vecs):
+        normals = np.cross(vecs[:,1:], vecs[:,:-1], axis=0)
+        normals = normals / np.linalg.norm(normals, axis=0)[np.newaxis, :]
+        # mats = mats = np.einsum("ji,ki->ijk", normals, normals) * -2
+        # np.einsum("ijj->ij", mats)[:] += 1
+        mats = -_generate_180_mats_from_uvecs(normals)
+        return mats
+
+    mats_0 = make_reflection_mat(vecs_0)
+    mats_1 = make_reflection_mat(vecs_1)
+    mats_2 = make_reflection_mat(vecs_2)
+    mats_3 = make_reflection_mat(vecs_3)
+
+    corner_mat_0 = np.matmul(mats_0[0], mats_3[0])
+    corner_mat_1 = np.matmul(mats_0[-1], mats_1[0])
+    corner_mat_2 = np.matmul(mats_2[-1], mats_1[-1])
+    corner_mat_3 = np.matmul(mats_2[0], mats_3[-1])
+
+    mats = np.concatenate((
+        corner_mat_0[np.newaxis],
+        mats_0,
+        corner_mat_1[np.newaxis],
+        mats_1,
+        corner_mat_2[np.newaxis],
+        mats_2[::-1],
+        corner_mat_3[np.newaxis],
+        mats_3[::-1],
+    ))
+    return mats
+
 
 def _2D_guess_bounds_first_pass(array):
     # average and normalise, boundary buffer represents edges and corners
@@ -504,7 +540,7 @@ def _2D_gb_buffer_inner(array_shape):
     return np.s_[:, x_i, y_i]
 
 
-def _2D_geuss_bounds(cube):
+def _2D_geuss_bounds(cube, rotate=True):
     lons = cube.coord(axis="X")
     lats = cube.coord(axis="Y")
 
@@ -515,10 +551,12 @@ def _2D_geuss_bounds(cube):
     xyz_array = _3d_xyz_from_latlon(lon_array, lat_array)
 
     result_xyz = _2D_guess_bounds_first_pass(xyz_array)
-    outer_inds = _2D_gb_buffer_outer(xyz_array.shape)
-    inner_inds = _2D_gb_buffer_inner(xyz_array.shape)
-
-    mats = _generate_180_mats_from_uvecs(result_xyz[outer_inds])
+    outer_inds = _2D_gb_buffer_outer(result_xyz.shape)
+    inner_inds = _2D_gb_buffer_inner(result_xyz.shape)
+    if rotate:
+        mats = _generate_180_mats_from_uvecs(result_xyz[outer_inds])
+    else:
+        mats = _generate_reflection_mats_from_uvecs(xyz_array)
     result_xyz[outer_inds] = _vectorised_matmul(mats, result_xyz[inner_inds])
 
     result_lon_bounds, result_lat_bounds = _latlon_from_xyz(result_xyz)
diff --git a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
index bc61b39a67..3e22d1acaa 100644
--- a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
+++ b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
@@ -312,3 +312,7 @@ def test_2D_guess_bounds():
     _2D_geuss_bounds(cube)
     assert cube.coord("latitude").is_contiguous()
     assert cube.coord("longitude").is_contiguous()
+
+    _2D_geuss_bounds(cube, rotate=False)
+    assert cube.coord("latitude").is_contiguous()
+    assert cube.coord("longitude").is_contiguous()

From 84815af47f1e3db94f64593c7565d4ccef55ee00 Mon Sep 17 00:00:00 2001
From: "stephen.worsley" <stephen.worsley@metoffice.gov.uk>
Date: Mon, 9 Mar 2026 16:48:48 +0000
Subject: [PATCH 04/14] add extrapolate and in place

---
 lib/iris/analysis/_grid_angles.py             | 75 +++++++---------
 .../cartography/test_gridcell_angles.py       | 90 +++++++++++++++++--
 2 files changed, 112 insertions(+), 53 deletions(-)

diff --git a/lib/iris/analysis/_grid_angles.py b/lib/iris/analysis/_grid_angles.py
index 1177fa1850..b5d132fc4b 100644
--- a/lib/iris/analysis/_grid_angles.py
+++ b/lib/iris/analysis/_grid_angles.py
@@ -14,6 +14,7 @@
 import numpy as np
 
 import iris
+from iris.coord_systems import GeogCS, RotatedGeogCS
 
 
 def _3d_xyz_from_latlon(lon, lat):
@@ -464,42 +465,6 @@ def _generate_180_mats_from_uvecs(uvecs):
     np.einsum("ijj->ij", mats)[:] -= 1
     return mats
 
-def _generate_reflection_mats_from_uvecs(xyz_array):
-    vecs_0 = xyz_array[:,:,0]
-    vecs_1 = xyz_array[:,-1,:]
-    vecs_2 = xyz_array[:,:,-1]
-    vecs_3 = xyz_array[:,0,:]
-
-    def make_reflection_mat(vecs):
-        normals = np.cross(vecs[:,1:], vecs[:,:-1], axis=0)
-        normals = normals / np.linalg.norm(normals, axis=0)[np.newaxis, :]
-        # mats = mats = np.einsum("ji,ki->ijk", normals, normals) * -2
-        # np.einsum("ijj->ij", mats)[:] += 1
-        mats = -_generate_180_mats_from_uvecs(normals)
-        return mats
-
-    mats_0 = make_reflection_mat(vecs_0)
-    mats_1 = make_reflection_mat(vecs_1)
-    mats_2 = make_reflection_mat(vecs_2)
-    mats_3 = make_reflection_mat(vecs_3)
-
-    corner_mat_0 = np.matmul(mats_0[0], mats_3[0])
-    corner_mat_1 = np.matmul(mats_0[-1], mats_1[0])
-    corner_mat_2 = np.matmul(mats_2[-1], mats_1[-1])
-    corner_mat_3 = np.matmul(mats_2[0], mats_3[-1])
-
-    mats = np.concatenate((
-        corner_mat_0[np.newaxis],
-        mats_0,
-        corner_mat_1[np.newaxis],
-        mats_1,
-        corner_mat_2[np.newaxis],
-        mats_2[::-1],
-        corner_mat_3[np.newaxis],
-        mats_3[::-1],
-    ))
-    return mats
-
 
 def _2D_guess_bounds_first_pass(array):
     # average and normalise, boundary buffer represents edges and corners
@@ -540,24 +505,46 @@ def _2D_gb_buffer_inner(array_shape):
     return np.s_[:, x_i, y_i]
 
 
-def _2D_geuss_bounds(cube, rotate=True):
+def _2D_guess_bounds(cube, extrapolate=True, in_place=False):
     lons = cube.coord(axis="X")
     lats = cube.coord(axis="Y")
+    h_dims = cube.coord_dims(lons)
+    assert h_dims == cube.coord_dims(lats)
+    assert len(h_dims) == 2
+
+    if lons.units != "degrees" or lats.units != "degrees":
+        msg = "Coordinate units are expected to be degrees."
+        raise ValueError(msg)
+    if not all(isinstance(coord.coord_system, GeogCS | RotatedGeogCS | None) for coord in [lats, lons]):
+        msg = "Coordinate systems are expected geodetic."
+        raise ValueError(msg)
+
+    if in_place:
+        _2D_guess_bounds_in_place(lons, lats, extrapolate=extrapolate)
 
-    # TODO: check units and coordsystem and dims and such
+    else:
+        new_lons = lons.copy()
+        new_lats = lats.copy()
+        _2D_guess_bounds_in_place(new_lons, new_lats, extrapolate=extrapolate)
+        cube.remove_coord(lons)
+        cube.remove_coord(lats)
+        cube.add_aux_coord(new_lons, h_dims)
+        cube.add_aux_coord(new_lats, h_dims)
+
+        return cube
 
+
+def _2D_guess_bounds_in_place(lons, lats, extrapolate=True):
     lon_array = lons.points
     lat_array = lats.points
     xyz_array = _3d_xyz_from_latlon(lon_array, lat_array)
 
     result_xyz = _2D_guess_bounds_first_pass(xyz_array)
-    outer_inds = _2D_gb_buffer_outer(result_xyz.shape)
-    inner_inds = _2D_gb_buffer_inner(result_xyz.shape)
-    if rotate:
+    if extrapolate:
+        outer_inds = _2D_gb_buffer_outer(result_xyz.shape)
+        inner_inds = _2D_gb_buffer_inner(result_xyz.shape)
         mats = _generate_180_mats_from_uvecs(result_xyz[outer_inds])
-    else:
-        mats = _generate_reflection_mats_from_uvecs(xyz_array)
-    result_xyz[outer_inds] = _vectorised_matmul(mats, result_xyz[inner_inds])
+        result_xyz[outer_inds] = _vectorised_matmul(mats, result_xyz[inner_inds])
 
     result_lon_bounds, result_lat_bounds = _latlon_from_xyz(result_xyz)
 
diff --git a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
index 3e22d1acaa..c8b37d4139 100644
--- a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
+++ b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
@@ -11,8 +11,9 @@
 import numpy as np
 import pytest
 
-from iris.analysis._grid_angles import _2D_geuss_bounds
-from iris.analysis.cartography import gridcell_angles
+from iris.analysis._grid_angles import _2D_guess_bounds
+from iris.analysis.cartography import gridcell_angles, _transform_xy
+from iris.coord_systems import RotatedGeogCS, Mercator
 from iris.coords import AuxCoord
 from iris.cube import Cube
 from iris.tests import _shared_utils
@@ -304,15 +305,86 @@ def test_fail_bad_method(self):
             self._check_multiple_orientations_and_latitudes(method="something_unknown")
 
 
-def test_2D_guess_bounds():
+def test_2D_guess_bounds_contiguity():
     cube = _2d_multicells_testcube()
     assert not cube.coord("latitude").is_contiguous()
     assert not cube.coord("longitude").is_contiguous()
 
-    _2D_geuss_bounds(cube)
-    assert cube.coord("latitude").is_contiguous()
-    assert cube.coord("longitude").is_contiguous()
+    result_extrap = _2D_guess_bounds(cube)
+    assert result_extrap.coord("latitude").is_contiguous()
+    assert result_extrap.coord("longitude").is_contiguous()
 
-    _2D_geuss_bounds(cube, rotate=False)
-    assert cube.coord("latitude").is_contiguous()
-    assert cube.coord("longitude").is_contiguous()
+    result_clipped = _2D_guess_bounds(cube, extrapolate=False)
+    assert result_clipped.coord("latitude").is_contiguous()
+    assert result_clipped.coord("longitude").is_contiguous()
+
+
+def test_2D_guess_bounds_rotational_equivalence():
+    # Check that _2D_guess_bounds is rotationally equivalent.
+    cube = _2d_multicells_testcube()
+
+    # Define a rotation with a pair of coordinate systems.
+    rotated_cs = RotatedGeogCS(20, 30).as_cartopy_crs()
+    unrotated_cs = RotatedGeogCS(0, 0).as_cartopy_crs()
+
+    # Guess the bounds before rotating the lat-lon points.
+    _2D_guess_bounds(cube, extrapolate=True, in_place=True)
+    lon_bounds_unrotated = cube.coord("longitude").bounds
+    lat_bounds_unrotated = cube.coord("latitude").bounds
+
+    # Rotate the guessed lat-lon bounds.
+    rotated_lon_bounds, rotated_lat_bounds = _transform_xy(
+        unrotated_cs, lon_bounds_unrotated.flatten(), lat_bounds_unrotated.flatten(), rotated_cs
+    )
+
+    # Rotate the lat-lon points.
+    lat = cube.coord("latitude")
+    lon = cube.coord("longitude")
+    lon.points, lat.points = _transform_xy(unrotated_cs, lon.points, lat.points, rotated_cs)
+
+    # guess the bounds after rotating the lat-lon points.
+    _2D_guess_bounds(cube, extrapolate=True, in_place=True)
+    lon_bounds_from_rotated_points = cube.coord("longitude").bounds
+    lat_bounds_from_rotated_points = cube.coord("latitude").bounds
+
+    # Check that the results are equivalent.
+    assert np.allclose(rotated_lon_bounds, lon_bounds_from_rotated_points.flatten())
+    assert np.allclose(rotated_lat_bounds, lat_bounds_from_rotated_points.flatten())
+
+
+def test_2D_guess_bounds_transpose_equivalence():
+    # Check that _2D_guess_bounds is transpose equivalent.
+    cube = _2d_multicells_testcube()
+    cube_transposed = _2d_multicells_testcube()
+    cube_transposed.transpose()
+
+    _2D_guess_bounds(cube, extrapolate=True, in_place=True)
+    _2D_guess_bounds(cube_transposed, extrapolate=True, in_place=True)
+
+    cube_transposed.transpose()
+
+    assert np.allclose(cube_transposed.coord("latitude").bounds, cube.coord("latitude").bounds)
+    assert np.allclose(cube_transposed.coord("longitude").bounds, cube.coord("longitude").bounds)
+
+
+def test_2D_guess_bounds_coord_systems():
+    rotated_cs = RotatedGeogCS(20, 30)
+    mercator_cs = Mercator()
+
+    rotated_cube = _2d_multicells_testcube()
+    rotated_cube.coord("latitude").coord_system = rotated_cs
+    rotated_cube.coord("latitude").standard_name = "grid_latitude"
+    rotated_cube.coord("longitude").coord_system = rotated_cs
+    rotated_cube.coord("longitude").standard_name = "grid_longitude"
+
+    _2D_guess_bounds(rotated_cube, in_place=True)
+
+    assert rotated_cube.coord("grid_latitude").is_contiguous()
+    assert rotated_cube.coord("grid_longitude").is_contiguous()
+
+    mercator_cube = _2d_multicells_testcube()
+    mercator_cube.coord("latitude").coord_system = mercator_cs
+    mercator_cube.coord("longitude").coord_system = mercator_cs
+
+    with pytest.raises(ValueError):
+        _2D_guess_bounds(mercator_cube)

From 9a6385c1b23bcd84041e781666b84e50bb8df35e Mon Sep 17 00:00:00 2001
From: "stephen.worsley" <stephen.worsley@metoffice.gov.uk>
Date: Wed, 11 Mar 2026 13:38:47 +0000
Subject: [PATCH 05/14] further transpose testing

---
 .../cartography/test_gridcell_angles.py       | 22 +++++++++++++++++--
 1 file changed, 20 insertions(+), 2 deletions(-)

diff --git a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
index c8b37d4139..72aadecc1d 100644
--- a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
+++ b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
@@ -356,15 +356,33 @@ def test_2D_guess_bounds_transpose_equivalence():
     # Check that _2D_guess_bounds is transpose equivalent.
     cube = _2d_multicells_testcube()
     cube_transposed = _2d_multicells_testcube()
+
+    def transpose_2D_coord(coord):
+        new_points = coord.points.transpose()
+        new_bounds = coord.bounds.transpose((1, 0, 2))[:,:,(0,3,2,1)]
+        new_coord = AuxCoord(new_points, bounds=new_bounds, standard_name=coord.standard_name, units=coord.units)
+        return new_coord
+
+
     cube_transposed.transpose()
 
+    new_lat = transpose_2D_coord(cube_transposed.coord("latitude"))
+    new_lon = transpose_2D_coord(cube_transposed.coord("longitude"))
+    cube_transposed.remove_coord("latitude")
+    cube_transposed.remove_coord("longitude")
+    cube_transposed.add_aux_coord(new_lat, (0, 1))
+    cube_transposed.add_aux_coord(new_lon, (0, 1))
+
     _2D_guess_bounds(cube, extrapolate=True, in_place=True)
     _2D_guess_bounds(cube_transposed, extrapolate=True, in_place=True)
 
     cube_transposed.transpose()
 
-    assert np.allclose(cube_transposed.coord("latitude").bounds, cube.coord("latitude").bounds)
-    assert np.allclose(cube_transposed.coord("longitude").bounds, cube.coord("longitude").bounds)
+    untransposed_lat = transpose_2D_coord(cube_transposed.coord("latitude"))
+    untransposed_lon = transpose_2D_coord(cube_transposed.coord("longitude"))
+
+    assert np.allclose(untransposed_lat.bounds, cube.coord("latitude").bounds)
+    assert np.allclose(untransposed_lon.bounds, cube.coord("longitude").bounds)
 
 
 def test_2D_guess_bounds_coord_systems():

From c5ee1e47e6a4220d5f742688175a7346b10b2304 Mon Sep 17 00:00:00 2001
From: "stephen.worsley" <stephen.worsley@metoffice.gov.uk>
Date: Fri, 13 Mar 2026 14:54:12 +0000
Subject: [PATCH 06/14] fix transpose test

---
 lib/iris/analysis/_grid_angles.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/iris/analysis/_grid_angles.py b/lib/iris/analysis/_grid_angles.py
index b5d132fc4b..f5dd6b805e 100644
--- a/lib/iris/analysis/_grid_angles.py
+++ b/lib/iris/analysis/_grid_angles.py
@@ -500,7 +500,7 @@ def _2D_gb_buffer_inner(array_shape):
         + ([1] * (y - 1))
     )
     y_i = (
-        ([1] * x) + list(range(1, y - 1)) + ([y - 1] * x) + list(range(1, y - 1))[::-1]
+        ([1] * x) + list(range(1, y - 1)) + ([y - 2] * x) + list(range(1, y - 1))[::-1]
     )
     return np.s_[:, x_i, y_i]
 

From e69d727ec69dc8242fb17b8cf31fdec6a0040bbc Mon Sep 17 00:00:00 2001
From: "stephen.worsley" <stephen.worsley@metoffice.gov.uk>
Date: Fri, 13 Mar 2026 15:06:12 +0000
Subject: [PATCH 07/14] ruff fixes

---
 lib/iris/analysis/_grid_angles.py             |  5 +++-
 .../cartography/test_gridcell_angles.py       | 25 +++++++++++++------
 2 files changed, 21 insertions(+), 9 deletions(-)

diff --git a/lib/iris/analysis/_grid_angles.py b/lib/iris/analysis/_grid_angles.py
index f5dd6b805e..506215be4d 100644
--- a/lib/iris/analysis/_grid_angles.py
+++ b/lib/iris/analysis/_grid_angles.py
@@ -515,7 +515,10 @@ def _2D_guess_bounds(cube, extrapolate=True, in_place=False):
     if lons.units != "degrees" or lats.units != "degrees":
         msg = "Coordinate units are expected to be degrees."
         raise ValueError(msg)
-    if not all(isinstance(coord.coord_system, GeogCS | RotatedGeogCS | None) for coord in [lats, lons]):
+    if not all(
+        isinstance(coord.coord_system, GeogCS | RotatedGeogCS | None)
+        for coord in [lats, lons]
+    ):
         msg = "Coordinate systems are expected geodetic."
         raise ValueError(msg)
 
diff --git a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
index 72aadecc1d..b5d1386854 100644
--- a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
+++ b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
@@ -12,8 +12,8 @@
 import pytest
 
 from iris.analysis._grid_angles import _2D_guess_bounds
-from iris.analysis.cartography import gridcell_angles, _transform_xy
-from iris.coord_systems import RotatedGeogCS, Mercator
+from iris.analysis.cartography import _transform_xy, gridcell_angles
+from iris.coord_systems import Mercator, RotatedGeogCS
 from iris.coords import AuxCoord
 from iris.cube import Cube
 from iris.tests import _shared_utils
@@ -334,13 +334,18 @@ def test_2D_guess_bounds_rotational_equivalence():
 
     # Rotate the guessed lat-lon bounds.
     rotated_lon_bounds, rotated_lat_bounds = _transform_xy(
-        unrotated_cs, lon_bounds_unrotated.flatten(), lat_bounds_unrotated.flatten(), rotated_cs
+        unrotated_cs,
+        lon_bounds_unrotated.flatten(),
+        lat_bounds_unrotated.flatten(),
+        rotated_cs,
     )
 
     # Rotate the lat-lon points.
     lat = cube.coord("latitude")
     lon = cube.coord("longitude")
-    lon.points, lat.points = _transform_xy(unrotated_cs, lon.points, lat.points, rotated_cs)
+    lon.points, lat.points = _transform_xy(
+        unrotated_cs, lon.points, lat.points, rotated_cs
+    )
 
     # guess the bounds after rotating the lat-lon points.
     _2D_guess_bounds(cube, extrapolate=True, in_place=True)
@@ -359,11 +364,15 @@ def test_2D_guess_bounds_transpose_equivalence():
 
     def transpose_2D_coord(coord):
         new_points = coord.points.transpose()
-        new_bounds = coord.bounds.transpose((1, 0, 2))[:,:,(0,3,2,1)]
-        new_coord = AuxCoord(new_points, bounds=new_bounds, standard_name=coord.standard_name, units=coord.units)
+        new_bounds = coord.bounds.transpose((1, 0, 2))[:, :, (0, 3, 2, 1)]
+        new_coord = AuxCoord(
+            new_points,
+            bounds=new_bounds,
+            standard_name=coord.standard_name,
+            units=coord.units,
+        )
         return new_coord
 
-
     cube_transposed.transpose()
 
     new_lat = transpose_2D_coord(cube_transposed.coord("latitude"))
@@ -404,5 +413,5 @@ def test_2D_guess_bounds_coord_systems():
     mercator_cube.coord("latitude").coord_system = mercator_cs
     mercator_cube.coord("longitude").coord_system = mercator_cs
 
-    with pytest.raises(ValueError):
+    with pytest.raises(ValueError, match="Coordinate systems are expected geodetic."):
         _2D_guess_bounds(mercator_cube)

From 1dd7bacfbf90a40542df015b46198807e4563683 Mon Sep 17 00:00:00 2001
From: "stephen.worsley" <stephen.worsley@metoffice.gov.uk>
Date: Thu, 16 Apr 2026 16:10:02 +0100
Subject: [PATCH 08/14] add public function to util

---
 lib/iris/analysis/_grid_angles.py             |  33 ----
 .../cartography/test_gridcell_angles.py       | 115 +------------
 .../tests/unit/util/test_guess_2D_bounds.py   | 160 ++++++++++++++++++
 lib/iris/util.py                              |  42 ++++-
 4 files changed, 202 insertions(+), 148 deletions(-)
 create mode 100644 lib/iris/tests/unit/util/test_guess_2D_bounds.py

diff --git a/lib/iris/analysis/_grid_angles.py b/lib/iris/analysis/_grid_angles.py
index 874fa5e751..f49d6234ba 100644
--- a/lib/iris/analysis/_grid_angles.py
+++ b/lib/iris/analysis/_grid_angles.py
@@ -14,7 +14,6 @@
 import numpy as np
 
 import iris
-from iris.coord_systems import GeogCS, RotatedGeogCS
 
 
 def _3d_xyz_from_latlon(lon, lat):
@@ -505,38 +504,6 @@ def _2D_gb_buffer_inner(array_shape):
     return np.s_[:, x_i, y_i]
 
 
-def _2D_guess_bounds(cube, extrapolate=True, in_place=False):
-    lons = cube.coord(axis="X")
-    lats = cube.coord(axis="Y")
-    h_dims = cube.coord_dims(lons)
-    assert h_dims == cube.coord_dims(lats)
-    assert len(h_dims) == 2
-
-    if lons.units != "degrees" or lats.units != "degrees":
-        msg = "Coordinate units are expected to be degrees."
-        raise ValueError(msg)
-    if not all(
-        isinstance(coord.coord_system, GeogCS | RotatedGeogCS | None)
-        for coord in [lats, lons]
-    ):
-        msg = "Coordinate systems are expected geodetic."
-        raise ValueError(msg)
-
-    if in_place:
-        _2D_guess_bounds_in_place(lons, lats, extrapolate=extrapolate)
-
-    else:
-        new_lons = lons.copy()
-        new_lats = lats.copy()
-        _2D_guess_bounds_in_place(new_lons, new_lats, extrapolate=extrapolate)
-        cube.remove_coord(lons)
-        cube.remove_coord(lats)
-        cube.add_aux_coord(new_lons, h_dims)
-        cube.add_aux_coord(new_lats, h_dims)
-
-        return cube
-
-
 def _2D_guess_bounds_in_place(lons, lats, extrapolate=True):
     lon_array = lons.points
     lat_array = lats.points
diff --git a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
index b5d1386854..4446fdb92a 100644
--- a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
+++ b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
@@ -12,8 +12,7 @@
 import pytest
 
 from iris.analysis._grid_angles import _2D_guess_bounds
-from iris.analysis.cartography import _transform_xy, gridcell_angles
-from iris.coord_systems import Mercator, RotatedGeogCS
+from iris.analysis.cartography import gridcell_angles
 from iris.coords import AuxCoord
 from iris.cube import Cube
 from iris.tests import _shared_utils
@@ -303,115 +302,3 @@ def test_fail_noncoord_coord(self):
     def test_fail_bad_method(self):
         with pytest.raises(ValueError, match="unrecognised cell_angle_boundpoints"):
             self._check_multiple_orientations_and_latitudes(method="something_unknown")
-
-
-def test_2D_guess_bounds_contiguity():
-    cube = _2d_multicells_testcube()
-    assert not cube.coord("latitude").is_contiguous()
-    assert not cube.coord("longitude").is_contiguous()
-
-    result_extrap = _2D_guess_bounds(cube)
-    assert result_extrap.coord("latitude").is_contiguous()
-    assert result_extrap.coord("longitude").is_contiguous()
-
-    result_clipped = _2D_guess_bounds(cube, extrapolate=False)
-    assert result_clipped.coord("latitude").is_contiguous()
-    assert result_clipped.coord("longitude").is_contiguous()
-
-
-def test_2D_guess_bounds_rotational_equivalence():
-    # Check that _2D_guess_bounds is rotationally equivalent.
-    cube = _2d_multicells_testcube()
-
-    # Define a rotation with a pair of coordinate systems.
-    rotated_cs = RotatedGeogCS(20, 30).as_cartopy_crs()
-    unrotated_cs = RotatedGeogCS(0, 0).as_cartopy_crs()
-
-    # Guess the bounds before rotating the lat-lon points.
-    _2D_guess_bounds(cube, extrapolate=True, in_place=True)
-    lon_bounds_unrotated = cube.coord("longitude").bounds
-    lat_bounds_unrotated = cube.coord("latitude").bounds
-
-    # Rotate the guessed lat-lon bounds.
-    rotated_lon_bounds, rotated_lat_bounds = _transform_xy(
-        unrotated_cs,
-        lon_bounds_unrotated.flatten(),
-        lat_bounds_unrotated.flatten(),
-        rotated_cs,
-    )
-
-    # Rotate the lat-lon points.
-    lat = cube.coord("latitude")
-    lon = cube.coord("longitude")
-    lon.points, lat.points = _transform_xy(
-        unrotated_cs, lon.points, lat.points, rotated_cs
-    )
-
-    # guess the bounds after rotating the lat-lon points.
-    _2D_guess_bounds(cube, extrapolate=True, in_place=True)
-    lon_bounds_from_rotated_points = cube.coord("longitude").bounds
-    lat_bounds_from_rotated_points = cube.coord("latitude").bounds
-
-    # Check that the results are equivalent.
-    assert np.allclose(rotated_lon_bounds, lon_bounds_from_rotated_points.flatten())
-    assert np.allclose(rotated_lat_bounds, lat_bounds_from_rotated_points.flatten())
-
-
-def test_2D_guess_bounds_transpose_equivalence():
-    # Check that _2D_guess_bounds is transpose equivalent.
-    cube = _2d_multicells_testcube()
-    cube_transposed = _2d_multicells_testcube()
-
-    def transpose_2D_coord(coord):
-        new_points = coord.points.transpose()
-        new_bounds = coord.bounds.transpose((1, 0, 2))[:, :, (0, 3, 2, 1)]
-        new_coord = AuxCoord(
-            new_points,
-            bounds=new_bounds,
-            standard_name=coord.standard_name,
-            units=coord.units,
-        )
-        return new_coord
-
-    cube_transposed.transpose()
-
-    new_lat = transpose_2D_coord(cube_transposed.coord("latitude"))
-    new_lon = transpose_2D_coord(cube_transposed.coord("longitude"))
-    cube_transposed.remove_coord("latitude")
-    cube_transposed.remove_coord("longitude")
-    cube_transposed.add_aux_coord(new_lat, (0, 1))
-    cube_transposed.add_aux_coord(new_lon, (0, 1))
-
-    _2D_guess_bounds(cube, extrapolate=True, in_place=True)
-    _2D_guess_bounds(cube_transposed, extrapolate=True, in_place=True)
-
-    cube_transposed.transpose()
-
-    untransposed_lat = transpose_2D_coord(cube_transposed.coord("latitude"))
-    untransposed_lon = transpose_2D_coord(cube_transposed.coord("longitude"))
-
-    assert np.allclose(untransposed_lat.bounds, cube.coord("latitude").bounds)
-    assert np.allclose(untransposed_lon.bounds, cube.coord("longitude").bounds)
-
-
-def test_2D_guess_bounds_coord_systems():
-    rotated_cs = RotatedGeogCS(20, 30)
-    mercator_cs = Mercator()
-
-    rotated_cube = _2d_multicells_testcube()
-    rotated_cube.coord("latitude").coord_system = rotated_cs
-    rotated_cube.coord("latitude").standard_name = "grid_latitude"
-    rotated_cube.coord("longitude").coord_system = rotated_cs
-    rotated_cube.coord("longitude").standard_name = "grid_longitude"
-
-    _2D_guess_bounds(rotated_cube, in_place=True)
-
-    assert rotated_cube.coord("grid_latitude").is_contiguous()
-    assert rotated_cube.coord("grid_longitude").is_contiguous()
-
-    mercator_cube = _2d_multicells_testcube()
-    mercator_cube.coord("latitude").coord_system = mercator_cs
-    mercator_cube.coord("longitude").coord_system = mercator_cs
-
-    with pytest.raises(ValueError, match="Coordinate systems are expected geodetic."):
-        _2D_guess_bounds(mercator_cube)
diff --git a/lib/iris/tests/unit/util/test_guess_2D_bounds.py b/lib/iris/tests/unit/util/test_guess_2D_bounds.py
new file mode 100644
index 0000000000..909d5694fd
--- /dev/null
+++ b/lib/iris/tests/unit/util/test_guess_2D_bounds.py
@@ -0,0 +1,160 @@
+# Copyright Iris contributors
+#
+# This file is part of Iris and is released under the BSD license.
+# See LICENSE in the root of the repository for full licensing details.
+"""Unit tests for the function :func:`iris.util.guess_2D_bounds`."""
+
+import numpy as np
+import pytest
+
+from iris.analysis.cartography import _transform_xy
+from iris.coord_systems import Mercator, RotatedGeogCS
+from iris.coords import AuxCoord
+from iris.tests.unit.analysis.cartography.test_gridcell_angles import (
+    _2d_multicells_testcube,
+)
+from iris.util import guess_2D_bounds
+
+
+def _2D_guess_bounds(cube, extrapolate=True, in_place=False):
+    lon = cube.coord(axis="X")
+    lat = cube.coord(axis="Y")
+    new_lon, new_lat = guess_2D_bounds(
+        lon, lat, extrapolate=extrapolate, in_place=in_place
+    )
+    if not in_place:
+        for old, new in [[lon, new_lon], [lat, new_lat]]:
+            dims = cube.coord_dims(old)
+            cube.remove_coord(old)
+            cube.add_aux_coord(new, dims)
+    return cube
+
+
+def test_2D_guess_bounds_contiguity():
+    cube = _2d_multicells_testcube()
+    assert not cube.coord("latitude").is_contiguous()
+    assert not cube.coord("longitude").is_contiguous()
+
+    result_extrap = _2D_guess_bounds(cube)
+    assert result_extrap.coord("latitude").is_contiguous()
+    assert result_extrap.coord("longitude").is_contiguous()
+
+    result_clipped = _2D_guess_bounds(cube, extrapolate=False)
+    assert result_clipped.coord("latitude").is_contiguous()
+    assert result_clipped.coord("longitude").is_contiguous()
+
+
+def test_2D_guess_bounds_rotational_equivalence():
+    # Check that _2D_guess_bounds is rotationally equivalent.
+    cube = _2d_multicells_testcube()
+
+    # Define a rotation with a pair of coordinate systems.
+    rotated_cs = RotatedGeogCS(20, 30).as_cartopy_crs()
+    unrotated_cs = RotatedGeogCS(0, 0).as_cartopy_crs()
+
+    # Guess the bounds before rotating the lat-lon points.
+    _2D_guess_bounds(cube, extrapolate=True, in_place=True)
+    lon_bounds_unrotated = cube.coord("longitude").bounds
+    lat_bounds_unrotated = cube.coord("latitude").bounds
+
+    # Rotate the guessed lat-lon bounds.
+    rotated_lon_bounds, rotated_lat_bounds = _transform_xy(
+        unrotated_cs,
+        lon_bounds_unrotated.flatten(),
+        lat_bounds_unrotated.flatten(),
+        rotated_cs,
+    )
+
+    # Rotate the lat-lon points.
+    lat = cube.coord("latitude")
+    lon = cube.coord("longitude")
+    lon.points, lat.points = _transform_xy(
+        unrotated_cs, lon.points, lat.points, rotated_cs
+    )
+
+    # guess the bounds after rotating the lat-lon points.
+    _2D_guess_bounds(cube, extrapolate=True, in_place=True)
+    lon_bounds_from_rotated_points = cube.coord("longitude").bounds
+    lat_bounds_from_rotated_points = cube.coord("latitude").bounds
+
+    # Check that the results are equivalent.
+    assert np.allclose(rotated_lon_bounds, lon_bounds_from_rotated_points.flatten())
+    assert np.allclose(rotated_lat_bounds, lat_bounds_from_rotated_points.flatten())
+
+
+def test_2D_guess_bounds_transpose_equivalence():
+    # Check that _2D_guess_bounds is transpose equivalent.
+    cube = _2d_multicells_testcube()
+    cube_transposed = _2d_multicells_testcube()
+
+    def transpose_2D_coord(coord):
+        new_points = coord.points.transpose()
+        new_bounds = coord.bounds.transpose((1, 0, 2))[:, :, (0, 3, 2, 1)]
+        new_coord = AuxCoord(
+            new_points,
+            bounds=new_bounds,
+            standard_name=coord.standard_name,
+            units=coord.units,
+        )
+        return new_coord
+
+    cube_transposed.transpose()
+
+    new_lat = transpose_2D_coord(cube_transposed.coord("latitude"))
+    new_lon = transpose_2D_coord(cube_transposed.coord("longitude"))
+    cube_transposed.remove_coord("latitude")
+    cube_transposed.remove_coord("longitude")
+    cube_transposed.add_aux_coord(new_lat, (0, 1))
+    cube_transposed.add_aux_coord(new_lon, (0, 1))
+
+    _2D_guess_bounds(cube, extrapolate=True, in_place=True)
+    _2D_guess_bounds(cube_transposed, extrapolate=True, in_place=True)
+
+    cube_transposed.transpose()
+
+    untransposed_lat = transpose_2D_coord(cube_transposed.coord("latitude"))
+    untransposed_lon = transpose_2D_coord(cube_transposed.coord("longitude"))
+
+    assert np.allclose(untransposed_lat.bounds, cube.coord("latitude").bounds)
+    assert np.allclose(untransposed_lon.bounds, cube.coord("longitude").bounds)
+
+
+def test_2D_guess_bounds_slice_equivalence():
+    # Extrapolation should approximate expected values from an extended regular grid when points are
+    # close enough together.
+    shrink_factor = 1000
+    cube = _2d_multicells_testcube()
+    for coord in cube.coords():
+        coord.points = coord.points / shrink_factor
+    sub_cube = cube[1:, 1:-1]
+    _2D_guess_bounds(cube)
+    _2D_guess_bounds(sub_cube)
+    assert np.allclose(
+        cube[1:, 1:-1].coord("latitude").bounds, sub_cube.coord("latitude").bounds
+    )
+    assert np.allclose(
+        cube[1:, 1:-1].coord("longitude").bounds, sub_cube.coord("longitude").bounds
+    )
+
+
+def test_2D_guess_bounds_coord_systems():
+    rotated_cs = RotatedGeogCS(20, 30)
+    mercator_cs = Mercator()
+
+    rotated_cube = _2d_multicells_testcube()
+    rotated_cube.coord("latitude").coord_system = rotated_cs
+    rotated_cube.coord("latitude").standard_name = "grid_latitude"
+    rotated_cube.coord("longitude").coord_system = rotated_cs
+    rotated_cube.coord("longitude").standard_name = "grid_longitude"
+
+    _2D_guess_bounds(rotated_cube, in_place=True)
+
+    assert rotated_cube.coord("grid_latitude").is_contiguous()
+    assert rotated_cube.coord("grid_longitude").is_contiguous()
+
+    mercator_cube = _2d_multicells_testcube()
+    mercator_cube.coord("latitude").coord_system = mercator_cs
+    mercator_cube.coord("longitude").coord_system = mercator_cs
+
+    with pytest.raises(ValueError, match="Coordinate systems are expected geodetic."):
+        _2D_guess_bounds(mercator_cube)
diff --git a/lib/iris/util.py b/lib/iris/util.py
index 551b5aeb68..d2abde989f 100644
--- a/lib/iris/util.py
+++ b/lib/iris/util.py
@@ -37,7 +37,7 @@
 from iris._lazy_data import is_lazy_data, is_lazy_masked_data
 from iris.common import SERVICES
 from iris.common.lenient import _lenient_client
-from iris.coord_systems import GeogCS
+from iris.coord_systems import GeogCS, RotatedGeogCS
 import iris.exceptions
 import iris.warnings
 
@@ -3110,6 +3110,46 @@ def _format(value: Any) -> str:
     return s
 
 
+def guess_2D_bounds(x, y, extrapolate=True, in_place=False):
+    """Guess the bounds of a pair of 2D coords.
+
+    Parameters
+    ----------
+    x : class:`~iris.coords.AuxCoord`
+        A "longitude" or "grid_longitude" coordinate.
+    y : class:`~iris.coords.AuxCoord`
+        A "latitude" or "grid_latitude" coordinate.
+    extrapolate : bool, default=True
+        If True, extend the edge bounds beyond the limits of the edge points.
+    in_place : bool, default=False
+        If True, modify the coordinate arguments in place.
+    """
+    from iris.analysis._grid_angles import _2D_guess_bounds_in_place
+
+    assert len(x.shape) == len(y.shape) == 2
+    assert x.standard_name in ("longitude", "grid_longitude")
+    assert y.standard_name in ("latitude", "grid_latitude")
+
+    if x.units != "degrees" or y.units != "degrees":
+        msg = "Coordinate units are expected to be degrees."
+        raise ValueError(msg)
+    if not all(
+        isinstance(coord.coord_system, GeogCS | RotatedGeogCS | None)
+        for coord in [x, y]
+    ):
+        msg = "Coordinate systems are expected geodetic."
+        raise ValueError(msg)
+
+    if in_place:
+        new_x = x
+        new_y = y
+    else:
+        new_x = x.copy()
+        new_y = y.copy()
+    _2D_guess_bounds_in_place(new_x, new_y, extrapolate=extrapolate)
+    return new_x, new_y
+
+
 @dataclass
 class CMLSettings(threading.local):
     """Settings for controlling the behaviour and formatting of the CML output.

From 2523f0874d64ad035af7bc91de2619e9b33986a0 Mon Sep 17 00:00:00 2001
From: "stephen.worsley" <stephen.worsley@metoffice.gov.uk>
Date: Thu, 23 Apr 2026 16:34:46 +0100
Subject: [PATCH 09/14] move guess_2D_bounds to iris.analysis.cartography

---
 lib/iris/analysis/_grid_angles.py             | 39 +++++++++++++++++
 lib/iris/analysis/cartography.py              |  3 +-
 .../cartography/test_gridcell_angles.py       |  1 -
 .../cartography}/test_guess_2D_bounds.py      |  4 +-
 lib/iris/util.py                              | 42 +------------------
 5 files changed, 44 insertions(+), 45 deletions(-)
 rename lib/iris/tests/unit/{util => analysis/cartography}/test_guess_2D_bounds.py (97%)

diff --git a/lib/iris/analysis/_grid_angles.py b/lib/iris/analysis/_grid_angles.py
index f49d6234ba..11964fad3b 100644
--- a/lib/iris/analysis/_grid_angles.py
+++ b/lib/iris/analysis/_grid_angles.py
@@ -14,6 +14,7 @@
 import numpy as np
 
 import iris
+from iris.coord_systems import GeogCS, RotatedGeogCS
 
 
 def _3d_xyz_from_latlon(lon, lat):
@@ -537,3 +538,41 @@ def _2D_guess_bounds_in_place(lons, lats, extrapolate=True):
         ],
         axis=2,
     )
+
+
+def guess_2D_bounds(x, y, extrapolate=True, in_place=False):
+    """Guess the bounds of a pair of 2D coords.
+
+    Parameters
+    ----------
+    x : class:`~iris.coords.AuxCoord`
+        A "longitude" or "grid_longitude" coordinate.
+    y : class:`~iris.coords.AuxCoord`
+        A "latitude" or "grid_latitude" coordinate.
+    extrapolate : bool, default=True
+        If True, extend the edge bounds beyond the limits of the edge points.
+    in_place : bool, default=False
+        If True, modify the coordinate arguments in place.
+    """
+    assert len(x.shape) == len(y.shape) == 2
+    assert x.standard_name in ("longitude", "grid_longitude")
+    assert y.standard_name in ("latitude", "grid_latitude")
+
+    if x.units != "degrees" or y.units != "degrees":
+        msg = "Coordinate units are expected to be degrees."
+        raise ValueError(msg)
+    if not all(
+        isinstance(coord.coord_system, GeogCS | RotatedGeogCS | None)
+        for coord in [x, y]
+    ):
+        msg = "Coordinate systems are expected geodetic."
+        raise ValueError(msg)
+
+    if in_place:
+        new_x = x
+        new_y = y
+    else:
+        new_x = x.copy()
+        new_y = y.copy()
+    _2D_guess_bounds_in_place(new_x, new_y, extrapolate=extrapolate)
+    return new_x, new_y
diff --git a/lib/iris/analysis/cartography.py b/lib/iris/analysis/cartography.py
index e88711d51f..c3f77cbb11 100644
--- a/lib/iris/analysis/cartography.py
+++ b/lib/iris/analysis/cartography.py
@@ -27,7 +27,7 @@
 from iris.util import _meshgrid
 import iris.warnings
 
-from ._grid_angles import gridcell_angles, rotate_grid_vectors
+from ._grid_angles import gridcell_angles, guess_2D_bounds, rotate_grid_vectors
 
 # List of contents to control Sphinx autodocs.
 # Unfortunately essential to get docs for the grid_angles functions.
@@ -39,6 +39,7 @@
     "get_xy_contiguous_bounded_grids",
     "get_xy_grids",
     "gridcell_angles",
+    "guess_2D_bounds",
     "project",
     "rotate_grid_vectors",
     "rotate_pole",
diff --git a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
index 4446fdb92a..32af84317c 100644
--- a/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
+++ b/lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py
@@ -11,7 +11,6 @@
 import numpy as np
 import pytest
 
-from iris.analysis._grid_angles import _2D_guess_bounds
 from iris.analysis.cartography import gridcell_angles
 from iris.coords import AuxCoord
 from iris.cube import Cube
diff --git a/lib/iris/tests/unit/util/test_guess_2D_bounds.py b/lib/iris/tests/unit/analysis/cartography/test_guess_2D_bounds.py
similarity index 97%
rename from lib/iris/tests/unit/util/test_guess_2D_bounds.py
rename to lib/iris/tests/unit/analysis/cartography/test_guess_2D_bounds.py
index 909d5694fd..b700bad25a 100644
--- a/lib/iris/tests/unit/util/test_guess_2D_bounds.py
+++ b/lib/iris/tests/unit/analysis/cartography/test_guess_2D_bounds.py
@@ -2,7 +2,7 @@
 #
 # This file is part of Iris and is released under the BSD license.
 # See LICENSE in the root of the repository for full licensing details.
-"""Unit tests for the function :func:`iris.util.guess_2D_bounds`."""
+"""Unit tests for the function :func:`iris.analysis.cartography.guess_2D_bounds`."""
 
 import numpy as np
 import pytest
@@ -13,7 +13,7 @@
 from iris.tests.unit.analysis.cartography.test_gridcell_angles import (
     _2d_multicells_testcube,
 )
-from iris.util import guess_2D_bounds
+from iris.analysis.cartography import guess_2D_bounds
 
 
 def _2D_guess_bounds(cube, extrapolate=True, in_place=False):
diff --git a/lib/iris/util.py b/lib/iris/util.py
index d2abde989f..551b5aeb68 100644
--- a/lib/iris/util.py
+++ b/lib/iris/util.py
@@ -37,7 +37,7 @@
 from iris._lazy_data import is_lazy_data, is_lazy_masked_data
 from iris.common import SERVICES
 from iris.common.lenient import _lenient_client
-from iris.coord_systems import GeogCS, RotatedGeogCS
+from iris.coord_systems import GeogCS
 import iris.exceptions
 import iris.warnings
 
@@ -3110,46 +3110,6 @@ def _format(value: Any) -> str:
     return s
 
 
-def guess_2D_bounds(x, y, extrapolate=True, in_place=False):
-    """Guess the bounds of a pair of 2D coords.
-
-    Parameters
-    ----------
-    x : class:`~iris.coords.AuxCoord`
-        A "longitude" or "grid_longitude" coordinate.
-    y : class:`~iris.coords.AuxCoord`
-        A "latitude" or "grid_latitude" coordinate.
-    extrapolate : bool, default=True
-        If True, extend the edge bounds beyond the limits of the edge points.
-    in_place : bool, default=False
-        If True, modify the coordinate arguments in place.
-    """
-    from iris.analysis._grid_angles import _2D_guess_bounds_in_place
-
-    assert len(x.shape) == len(y.shape) == 2
-    assert x.standard_name in ("longitude", "grid_longitude")
-    assert y.standard_name in ("latitude", "grid_latitude")
-
-    if x.units != "degrees" or y.units != "degrees":
-        msg = "Coordinate units are expected to be degrees."
-        raise ValueError(msg)
-    if not all(
-        isinstance(coord.coord_system, GeogCS | RotatedGeogCS | None)
-        for coord in [x, y]
-    ):
-        msg = "Coordinate systems are expected geodetic."
-        raise ValueError(msg)
-
-    if in_place:
-        new_x = x
-        new_y = y
-    else:
-        new_x = x.copy()
-        new_y = y.copy()
-    _2D_guess_bounds_in_place(new_x, new_y, extrapolate=extrapolate)
-    return new_x, new_y
-
-
 @dataclass
 class CMLSettings(threading.local):
     """Settings for controlling the behaviour and formatting of the CML output.

From 6b8261fabb510f65297fec7ba966404a2c22bea3 Mon Sep 17 00:00:00 2001
From: "stephen.worsley" <stephen.worsley@metoffice.gov.uk>
Date: Thu, 23 Apr 2026 16:46:40 +0100
Subject: [PATCH 10/14] sort imports

---
 .../tests/unit/analysis/cartography/test_guess_2D_bounds.py    | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/lib/iris/tests/unit/analysis/cartography/test_guess_2D_bounds.py b/lib/iris/tests/unit/analysis/cartography/test_guess_2D_bounds.py
index b700bad25a..1300cc41eb 100644
--- a/lib/iris/tests/unit/analysis/cartography/test_guess_2D_bounds.py
+++ b/lib/iris/tests/unit/analysis/cartography/test_guess_2D_bounds.py
@@ -7,13 +7,12 @@
 import numpy as np
 import pytest
 
-from iris.analysis.cartography import _transform_xy
+from iris.analysis.cartography import _transform_xy, guess_2D_bounds
 from iris.coord_systems import Mercator, RotatedGeogCS
 from iris.coords import AuxCoord
 from iris.tests.unit.analysis.cartography.test_gridcell_angles import (
     _2d_multicells_testcube,
 )
-from iris.analysis.cartography import guess_2D_bounds
 
 
 def _2D_guess_bounds(cube, extrapolate=True, in_place=False):

From 88976b5783cbc79843ae82097c56e489dde5cffb Mon Sep 17 00:00:00 2001
From: "stephen.worsley" <stephen.worsley@metoffice.gov.uk>
Date: Fri, 1 May 2026 16:57:28 +0100
Subject: [PATCH 11/14] add tests and documentation

---
 lib/iris/analysis/_grid_angles.py             | 55 +++++++++++++++++--
 .../cartography/test_guess_2D_bounds.py       | 47 ++++++++++++++++
 2 files changed, 97 insertions(+), 5 deletions(-)

diff --git a/lib/iris/analysis/_grid_angles.py b/lib/iris/analysis/_grid_angles.py
index 11964fad3b..b850c3a16b 100644
--- a/lib/iris/analysis/_grid_angles.py
+++ b/lib/iris/analysis/_grid_angles.py
@@ -543,20 +543,65 @@ def _2D_guess_bounds_in_place(lons, lats, extrapolate=True):
 def guess_2D_bounds(x, y, extrapolate=True, in_place=False):
     """Guess the bounds of a pair of 2D coords.
 
+    Bounds on a 2D coordinate are structured so that each point of the coordinate
+    has 4 associated bounds, with each pair of neighbouring points sharing 2 of their
+    bounds in common. We can categorise these bounds by how many points they are
+    associated with:
+    - Internal bounds, belonging to 4 surrounding points.
+    - Edge bounds, belonging only to 2 points on the edge of the grid.
+    - Corner bounds, belonging only to 1 point on the corner of the grid.
+    Since each of these categories of bounds has a different number of points associated
+    with them, they must also be calculated differently:
+    - Each internal bound is calculated from the 4 surrounding points. The surrounding
+      points are first transformed from 2D lat-lon space to 3D space on the surface of a
+      unit sphere. The average of these points in 3D space is then taken. This point is then
+      projected onto that unit sphere and transformed back to 2D lat-lon space. Note that in
+      the edge case where the 3D average is precisely the origin, it will not be possible to
+      project onto the unit sphere and an error will be raised.
+    - The calculation for edge bounds depends on if the `extrapolate` keyword is True or
+      False. When `extrapolate` is False, the edge bound is calculated as the midpoint of
+      the 2 neighbouring bounds. This calculation is done, as above, by converting into 3D
+      space, averaging, and then converting back to lat-lon space. When `extrapolate` is True,
+      edge bounds are calculated using the neighbouring internal bound. The edge bound and the
+      neighbouring internal bound will be precisely the 2 bounds which the edge points share
+      in common. The extrapolated edge bound is calculated by rotating the internal bound 180
+      degrees about the midpoint between the 2 neighbouring points calculated previously. This
+      rotation is done in 3D space about the axis defined by the line which passes through the
+      midpoint and the origin. This is equivalent to defining the edge bound so that the
+      midpoint between the two bounds is the same as the midpoint between the two points.
+    - The calculation for the corner bounds similarly depends on `extrapolate`. When
+      `extrapolate` is False, each corner bound is precisely equal to its associated corner
+      point. When `extrapolate` is True, we calculate the corner bound by taking the only
+      internal bound associated with the corner point and, as above, rotating it by 180
+      degrees about the corner point.
+
+
     Parameters
     ----------
     x : class:`~iris.coords.AuxCoord`
-        A "longitude" or "grid_longitude" coordinate.
+        A "longitude" or "grid_longitude" coordinate. Coordinate must be 2D.
     y : class:`~iris.coords.AuxCoord`
-        A "latitude" or "grid_latitude" coordinate.
+        A "latitude" or "grid_latitude" coordinate. Coordinate must be 2D.
     extrapolate : bool, default=True
         If True, extend the edge bounds beyond the limits of the edge points.
     in_place : bool, default=False
         If True, modify the coordinate arguments in place.
     """
-    assert len(x.shape) == len(y.shape) == 2
-    assert x.standard_name in ("longitude", "grid_longitude")
-    assert y.standard_name in ("latitude", "grid_latitude")
+    if x.shape != y.shape:
+        msg = "Coordinates do not have the same shape."
+        raise ValueError(msg)
+    if len(x.shape) != 2:
+        msg = "Coordinates are not 2D."
+        raise ValueError(msg)
+    if x.shape[0] < 2 or x.shape[1] < 2:
+        msg = "Coordinates must have length at least 2 in each dimension."
+        raise ValueError(msg)
+    if x.standard_name not in ("longitude", "grid_longitude"):
+        msg = "X coordinate is not 'longitude' or 'grid_longitude'."
+        raise ValueError(msg)
+    if y.standard_name not in ("latitude", "grid_latitude"):
+        msg = "Y coordinate is not 'latitude' or 'grid_latitude'."
+        raise ValueError(msg)
 
     if x.units != "degrees" or y.units != "degrees":
         msg = "Coordinate units are expected to be degrees."
diff --git a/lib/iris/tests/unit/analysis/cartography/test_guess_2D_bounds.py b/lib/iris/tests/unit/analysis/cartography/test_guess_2D_bounds.py
index 1300cc41eb..45cde05a3c 100644
--- a/lib/iris/tests/unit/analysis/cartography/test_guess_2D_bounds.py
+++ b/lib/iris/tests/unit/analysis/cartography/test_guess_2D_bounds.py
@@ -157,3 +157,50 @@ def test_2D_guess_bounds_coord_systems():
 
     with pytest.raises(ValueError, match="Coordinate systems are expected geodetic."):
         _2D_guess_bounds(mercator_cube)
+
+
+def test_invalid_coords_1D():
+    lat_1D = AuxCoord(np.arange(5), standard_name="latitude")
+    lon_1D = AuxCoord(np.arange(5), standard_name="longitude")
+    with pytest.raises(ValueError, match="Coordinates are not 2D."):
+        guess_2D_bounds(lon_1D, lat_1D)
+
+
+def test_invalid_coords_shape():
+    lat = AuxCoord(np.ones([2, 3]), standard_name="latitude")
+    lon = AuxCoord(np.ones([3, 2]), standard_name="longitude")
+    with pytest.raises(ValueError, match="Coordinates do not have the same shape."):
+        guess_2D_bounds(lon, lat)
+
+
+def test_invalid_coords_size():
+    expected_msg = "Coordinates must have length at least 2 in each dimension."
+
+    lat = AuxCoord(np.ones([1, 3]), standard_name="latitude")
+    lon = AuxCoord(np.ones([1, 3]), standard_name="longitude")
+    with pytest.raises(ValueError, match=expected_msg):
+        guess_2D_bounds(lon, lat)
+
+    lat = AuxCoord(np.ones([3, 1]), standard_name="latitude")
+    lon = AuxCoord(np.ones([3, 1]), standard_name="longitude")
+    with pytest.raises(ValueError, match=expected_msg):
+        guess_2D_bounds(lon, lat)
+
+
+def test_invalid_coords_name():
+    lat_valid = AuxCoord(np.ones([2, 3]), standard_name="latitude")
+    lon_valid = AuxCoord(np.ones([2, 3]), standard_name="longitude")
+    lat_invalid = AuxCoord(np.ones([2, 3]), long_name="lat")
+    lon_invalid = AuxCoord(np.ones([2, 3]), long_name="long")
+    with pytest.raises(
+        ValueError, match="X coordinate is not 'longitude' or 'grid_longitude'."
+    ):
+        guess_2D_bounds(lon_invalid, lat_valid)
+    with pytest.raises(
+        ValueError, match="X coordinate is not 'longitude' or 'grid_longitude'."
+    ):
+        guess_2D_bounds(lat_valid, lon_valid)
+    with pytest.raises(
+        ValueError, match="Y coordinate is not 'latitude' or 'grid_latitude'."
+    ):
+        guess_2D_bounds(lon_valid, lat_invalid)

From ac7e08970ac179bd6046552d767602fcdd90f716 Mon Sep 17 00:00:00 2001
From: "stephen.worsley" <stephen.worsley@metoffice.gov.uk>
Date: Thu, 7 May 2026 11:13:23 +0100
Subject: [PATCH 12/14] add comments

---
 lib/iris/analysis/_grid_angles.py | 76 +++++++++++++++++++++++++++----
 1 file changed, 68 insertions(+), 8 deletions(-)

diff --git a/lib/iris/analysis/_grid_angles.py b/lib/iris/analysis/_grid_angles.py
index b850c3a16b..3c3522294d 100644
--- a/lib/iris/analysis/_grid_angles.py
+++ b/lib/iris/analysis/_grid_angles.py
@@ -457,30 +457,66 @@ def rotate_grid_vectors(u_cube, v_cube, grid_angles_cube=None, grid_angles_kwarg
 
 
 def _vectorised_matmul(mats, vecs):
-    return np.einsum("ijk,ji->ki", mats, vecs)
+    # We interpret the 3D array `mats` as if it were a list of matrices varying over
+    # its last dimension so that mats[:,:,i] represents a single matrix. We consider
+    # the 2D array `vecs` to represent a list of vectors varying over its last dimension
+    # so that vecs[:,i] represents a single vector. The output will be such that for:
+    # result[:,i] == mats[:,:,i] @ vecs[:,i].
+    return np.einsum("jki,ji->ki", mats, vecs)
 
 
 def _generate_180_mats_from_uvecs(uvecs):
-    mats = np.einsum("ji,ki->ijk", uvecs, uvecs) * 2
-    np.einsum("ijj->ij", mats)[:] -= 1
+    # Generates a 3D array representing a list of matrices which can be used by the
+    # function _vectorised_matmul. Both the input `uvecs` and the `result` vary in their
+    # last dimension so that uvecs[:,i] is a length 3 vector corresponding to the
+    # rotation matrix result[:,:,i], where this is a rotation of 180 degrees about the
+    # axis passing through uvecs[:,i] and the origin.
+    # If the vector uvecs[:,i] = (x,y,z), the equivalent rotation matrix result[:,:,i]
+    # will be:
+    # | 2x^2 - 1   2xy      2xz    |
+    # |   2xy    2y^2 - 1   2yz    |
+    # |   2xz      2yz    2z^2 - 1 |
+    mats = np.einsum("ji,ki->jki", uvecs, uvecs) * 2
+
+    # At this point the matrix mats[:,:,i] will be:
+    # |   2x^2     2xy      2xz    |
+    # |   2xy      2y^2     2yz    |
+    # |   2xz      2yz      2z^2   |
+    # to achieve the desired result, we take one from the diagonal.
+    np.einsum("jji->ji", mats)[:] -= 1
     return mats
 
 
 def _2D_guess_bounds_first_pass(array):
-    # average and normalise, boundary buffer represents edges and corners
+    # Average and normalise sets of neighbouring points. This calculation actually
+    # normalises the sum rather than the average, since this is equivalent.
+    # Coners of the resulting array will be the sum of only a single corner of the
+    # source array. Edges of the resulting array will be the midpoint between edge
+    # points in the source. Internal points in the resulting array will correspond
+    # to internal bounds in the final result of `guess_2D_bounds`.
+    # Note: if `extrapolate` is False, this array will contain all the bounds in the
+    # final result of `guess_2D_bounds`.
     result_array = np.zeros((array.shape[0], array.shape[1] + 1, array.shape[2] + 1))
     pads = ((0, 1), (1, 0))
     for pad_i in pads:
         for pad_j in pads:
             result_array += np.pad(array, ((0, 0), pad_i, pad_j))
 
-    # normalise
+    # Normalise
     result_array /= np.linalg.norm(result_array, ord=2, axis=0)[np.newaxis, ...]
     return result_array
 
 
 def _2D_gb_buffer_outer(array_shape):
-    # return appropriate numpy slice for outer halo
+    # Return appropriate numpy slice for outer halo.
+    # This slice preserves the first dimension, which captures the 3D vector points and
+    # lists a series of indices in the next 2 dimensions.
+    # Each index in this slice corresponds to a corner or edge bound.
+
+    # This halo starts at the index [:, 0, 0], the next set of indices increase in the
+    # second dimension until the index [:, -1, 0], then the last dimension increases
+    # until it reaches the index [:, -1, -1], the remaining indices follow the edge of
+    # the array anit-clockwise until it reaches the last index at [:, 0, 1].
     _, x, y = array_shape
     x_i = list(range(x)) + ([x - 1] * (y - 2)) + list(range(x))[::-1] + ([0] * (y - 2))
     y_i = (
@@ -490,7 +526,23 @@ def _2D_gb_buffer_outer(array_shape):
 
 
 def _2D_gb_buffer_inner(array_shape):
-    # return appropriate numpy slice for inner halo
+    # Return appropriate numpy slice for inner halo.
+    # This slice preserves the first dimension, which captures the 3D vector points and
+    # lists a series of indices in the next 2 dimensions.
+    # Each index in this slice corresponds to an internal bound neighbouring a corner or
+    # edge bound.
+    # For every index in the outer halo, this gives the nearest index not in the outer halo.
+    # Note: the internal bounds which are nearest to the corner bounds ar each the nearest
+    # bound for at least two other edge or corner bounds. Therefore, these indices will
+    # occur at least 3 times in the returned slice.
+
+    # This halo starts at the index [:, 1, 1], the next index is also [:, 1, 1]. The
+    # next set of indices increase in the second dimension until the index [:, -2, 1],
+    # this index is repeated twice. In the edge case where the index [:, 1, 1] is
+    # equivalent to [:, -2, 1] the first two copies of the index will be followed
+    # immediately by two more copies. Then the last dimension increases until it reaches
+    # the index [:, -2, -2], repeating twice again, the remaining indices follow this
+    # pattern until it reaches the last index at [:, 1, 1].
     _, x, y = array_shape
     x_i = (
         [1]
@@ -508,18 +560,26 @@ def _2D_gb_buffer_inner(array_shape):
 def _2D_guess_bounds_in_place(lons, lats, extrapolate=True):
     lon_array = lons.points
     lat_array = lats.points
+    # Convert from lat-lon to 3D space.
     xyz_array = _3d_xyz_from_latlon(lon_array, lat_array)
 
+    # Create an array with internal-bounds, edge-midpoints, and corner-points.
     result_xyz = _2D_guess_bounds_first_pass(xyz_array)
     if extrapolate:
+        # Generate slice of edge-midpoints/edge-bounds and corner-points/corner-bounds.
         outer_inds = _2D_gb_buffer_outer(result_xyz.shape)
+        # Generate slice of internal-bounds which neighbour the above bounds.
         inner_inds = _2D_gb_buffer_inner(result_xyz.shape)
+        # Generate rotation matrices about corner-points and edge-midpoints
         mats = _generate_180_mats_from_uvecs(result_xyz[outer_inds])
+        # Rotate internal-bounds about their corresponding corner-point/edge-midpoint.
+        # Replace the corner-point/edge-midpoint with this result.
         result_xyz[outer_inds] = _vectorised_matmul(mats, result_xyz[inner_inds])
 
+    # Convert back from 3D to lat-lon.
     result_lon_bounds, result_lat_bounds = _latlon_from_xyz(result_xyz)
 
-    # add these bounds cf style
+    # Reformat these bounds as CF style bounds.
     lons.bounds = np.stack(
         [
             result_lon_bounds[:-1, :-1],

From 922a01e249b2c2e2a4f6758b4c438f6908d97969 Mon Sep 17 00:00:00 2001
From: "stephen.worsley" <stephen.worsley@metoffice.gov.uk>
Date: Mon, 11 May 2026 12:07:40 +0100
Subject: [PATCH 13/14] fix typo

---
 lib/iris/analysis/_grid_angles.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/iris/analysis/_grid_angles.py b/lib/iris/analysis/_grid_angles.py
index 3c3522294d..9e9c02df15 100644
--- a/lib/iris/analysis/_grid_angles.py
+++ b/lib/iris/analysis/_grid_angles.py
@@ -490,7 +490,7 @@ def _generate_180_mats_from_uvecs(uvecs):
 def _2D_guess_bounds_first_pass(array):
     # Average and normalise sets of neighbouring points. This calculation actually
     # normalises the sum rather than the average, since this is equivalent.
-    # Coners of the resulting array will be the sum of only a single corner of the
+    # Corners of the resulting array will be the sum of only a single corner of the
     # source array. Edges of the resulting array will be the midpoint between edge
     # points in the source. Internal points in the resulting array will correspond
     # to internal bounds in the final result of `guess_2D_bounds`.

From b6fbd515850752e04420602962019d9f18e85662 Mon Sep 17 00:00:00 2001
From: "stephen.worsley" <stephen.worsley@metoffice.gov.uk>
Date: Mon, 11 May 2026 13:21:06 +0100
Subject: [PATCH 14/14] fix docstring bullets

---
 lib/iris/analysis/_grid_angles.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/lib/iris/analysis/_grid_angles.py b/lib/iris/analysis/_grid_angles.py
index 9e9c02df15..f77dd8cbfe 100644
--- a/lib/iris/analysis/_grid_angles.py
+++ b/lib/iris/analysis/_grid_angles.py
@@ -607,11 +607,14 @@ def guess_2D_bounds(x, y, extrapolate=True, in_place=False):
     has 4 associated bounds, with each pair of neighbouring points sharing 2 of their
     bounds in common. We can categorise these bounds by how many points they are
     associated with:
+
     - Internal bounds, belonging to 4 surrounding points.
     - Edge bounds, belonging only to 2 points on the edge of the grid.
     - Corner bounds, belonging only to 1 point on the corner of the grid.
+
     Since each of these categories of bounds has a different number of points associated
     with them, they must also be calculated differently:
+
     - Each internal bound is calculated from the 4 surrounding points. The surrounding
       points are first transformed from 2D lat-lon space to 3D space on the surface of a
       unit sphere. The average of these points in 3D space is then taken. This point is then