From 0126dc5c54e09f53cdb41568e519e4ff45757019 Mon Sep 17 00:00:00 2001 From: "stephen.worsley" Date: Mon, 29 Sep 2025 22:07:24 +0100 Subject: [PATCH 1/2] 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 2/2] [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() -