Index MeshCoord using new MeshXY or _MeshIndexSet

lib/iris/analysis/_interpolation.py

@@ -13,6 +13,7 @@
 import numpy.ma as ma
 
 from iris.coords import AuxCoord, DimCoord
+from iris.exceptions import MonotonicityError
 import iris.util
 
 _DEFAULT_DTYPE = np.float16
@@ -456,7 +457,7 @@ def _validate(self):
                 # Check monotonic.
                 if not iris.util.monotonic(coord.points, strict=True):
                     msg = "Cannot interpolate over the non-monotonic coordinate {}."
-                    raise ValueError(msg.format(coord.name()))
+                    raise MonotonicityError(msg.format(coord.name()))
 
     def _interpolated_dtype(self, dtype):
         """Determine the minimum base dtype required by the underlying interpolator."""

lib/iris/coords.py

@@ -2224,7 +2224,7 @@ def _guess_bounds(self, bound_position=0.5, monthly=False, yearly=False):
         # XXX Consider moving into DimCoord
         # ensure we have monotonic points
         if not self.is_monotonic():
-            raise ValueError(
+            raise iris.exceptions.MonotonicityError(
                 "Need monotonic points to generate bounds for %s" % self.name()
             )
 
@@ -2771,7 +2771,9 @@ def _new_points_requirements(self, points):
             raise TypeError(emsg.format(self.name(), self.__class__.__name__))
         if points.size > 1 and not iris.util.monotonic(points, strict=True):
             emsg = "The {!r} {} points array must be strictly monotonic."
-            raise ValueError(emsg.format(self.name(), self.__class__.__name__))
+            raise iris.exceptions.MonotonicityError(
+                emsg.format(self.name(), self.__class__.__name__)
+            )
 
     @property
     def _values(self):
@@ -2848,7 +2850,7 @@ def _new_bounds_requirements(self, bounds):
                     )
                     if not monotonic:
                         emsg = "The {!r} {} bounds array must be strictly monotonic."
-                        raise ValueError(
+                        raise iris.exceptions.MonotonicityError(
                             emsg.format(self.name(), self.__class__.__name__)
                         )
                     directions.add(direction)
@@ -2858,7 +2860,9 @@ def _new_bounds_requirements(self, bounds):
                         "The direction of monotonicity for {!r} {} must "
                         "be consistent across all bounds."
                     )
-                    raise ValueError(emsg.format(self.name(), self.__class__.__name__))
+                    raise iris.exceptions.MonotonicityError(
+                        emsg.format(self.name(), self.__class__.__name__)
+                    )
 
                 if n_bounds == 2:
                     # Make ordering of bounds consistent with coord's direction

lib/iris/cube.py

@@ -5,6 +5,8 @@
 
 """Classes for representing multi-dimensional data with metadata."""
 
+from __future__ import annotations
+
 from collections import OrderedDict
 import copy
 from copy import deepcopy
@@ -15,6 +17,7 @@
     Container,
     Iterable,
     Iterator,
+    Literal,
     Mapping,
     MutableMapping,
     Optional,
@@ -2364,7 +2367,7 @@ def _any_meshcoord(self):
         return result
 
     @property
-    def mesh(self):
+    def mesh(self) -> "iris.mesh.MeshXY":
         r"""Return the unstructured :class:`~iris.mesh.MeshXY` associated with the cube.
 
         Return the unstructured :class:`~iris.mesh.MeshXY`
@@ -2386,7 +2389,7 @@ def mesh(self):
         return result
 
     @property
-    def location(self):
+    def location(self) -> Literal["face", "edge", "node", None]:
         r"""Return the mesh "location" of the cube data.
 
         Return the mesh "location" of the cube data, if the cube has any
@@ -2890,8 +2893,7 @@ def new_ancillary_variable_dims(av_):
             coord_keys = tuple([full_slice[dim] for dim in self.coord_dims(coord)])
             try:
                 new_coord = coord[coord_keys]
-            except ValueError:
-                # TODO make this except more specific to catch monotonic error
+            except iris.exceptions.MonotonicityError:
                 # Attempt to slice it by converting to AuxCoord first
                 new_coord = iris.coords.AuxCoord.from_coord(coord)[coord_keys]
             cube.add_aux_coord(new_coord, new_coord_dims(coord))
@@ -2910,8 +2912,7 @@ def new_ancillary_variable_dims(av_):
                     cube.add_aux_coord(new_coord, new_dims)
                 else:
                     cube.add_dim_coord(new_coord, new_dims)
-            except ValueError:
-                # TODO make this except more specific to catch monotonic error
+            except iris.exceptions.MonotonicityError:
                 # Attempt to slice it by converting to AuxCoord first
                 new_coord = iris.coords.AuxCoord.from_coord(coord)[coord_keys]
                 cube.add_aux_coord(new_coord, new_dims)

lib/iris/exceptions.py

@@ -161,3 +161,9 @@ class CannotAddError(ValueError):
     """Raised when an object (e.g. coord) cannot be added to a :class:`~iris.cube.Cube`."""
 
     pass
+
+
+class MonotonicityError(ValueError):
+    """Raised when a coordinate values are not monotonic."""
+
+    pass

lib/iris/experimental/__init__.py

@@ -8,3 +8,68 @@
 codebase. The code is expected to graduate, eventually, to "full status".
 
 """
+
+# import enum
+# import threading
+#
+#
+# class MeshCoordIndexing(threading.local):
+#     class Settings(enum.Enum):
+#         """Setting for what is returned when a :class:`~iris.mesh.MeshCoord` is indexed.
+#
+#         Use via the run-time setting :data:`~iris.experimental.MESH_COORD_INDEXING`.
+#         """
+#         # TODO: thread safety?
+#         # TODO: docstring example
+#
+#         AUX_COORD = enum.auto()
+#         """Convert the ``MeshCoord`` to a :class:`~iris.coords.AuxCoord`."""
+#
+#         NEW_MESH = enum.auto()
+#         """
+#         Index the ``MeshCoord`` :attr:`~iris.mesh.MeshCoord.mesh`, returning a new
+#         :class:`~iris.mesh.MeshXY` instance.
+#         """
+#
+#         MESH_INDEX_SET = enum.auto()
+#         """
+#         Index the ``MeshCoord`` :attr:`~iris.mesh.MeshCoord.mesh`, returning a
+#         :class:`iris.mesh.components._MeshIndexSet` instance. This is a read-only
+#         'view' onto the original :class:`~iris.mesh.MeshXY` instance, and is
+#         currently experimental, subject to future iterative changes.
+#         """
+#
+#     def __init__(self):
+#         self._setting = MeshCoordIndexing.Settings.AUX_COORD
+#
+#     @property
+#     def setting(self):
+#         return self._setting
+#
+#     @setting.setter
+#     def setting(self, value):
+#         self._setting = value
+#
+#     # class Setting(threading.local):
+#     #     # setting: "MeshCoordIndexing"
+#     #
+#     #     def __init__(self):
+#     #         self.set(MeshCoordIndexing.AUX_COORD)
+#     #
+#     #     def set(self, value):
+#     #         if isinstance(value, MeshCoordIndexing):
+#     #             self.setting = value
+#     #         else:
+#     #             message = (
+#     #                 "MeshCoordIndexing.Setting.set() must be called with a "
+#     #                 f"MeshCoordIndexing enum value. Instead, got: {value}"
+#     #             )
+#     #             raise ValueError(message)
+#
+#
+# MESH_COORD_INDEXING = MeshCoordIndexing()
+# """
+# Run-time setting for alternative modes of indexing a
+# :class:`~iris.mesh.MeshCoord`. See
+# :class:`~iris.experimental.MeshCoordIndexing` for more.
+# """

lib/iris/experimental/mesh_coord_indexing.py

@@ -0,0 +1,135 @@
+# 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.
+"""Experimental module for alternative modes of indexing a :class:`~iris.mesh.MeshCoord`.
+
+:class:`Options` describes the available indexing modes.
+
+Select the desired option using the run-time setting :data:`SETTING`.
+
+Examples
+--------
+.. testsetup::
+
+    import iris
+    my_mesh_cube = iris.load_cube(iris.sample_data_path("mesh_C4_synthetic_float.nc"))
+
+    # Remove non-compliant content.
+    my_mesh = my_mesh_cube.mesh
+    wanted_roles = ["edge_node_connectivity", "face_node_connectivity"]
+    for conn in my_mesh.all_connectivities:
+        if conn is not None and conn.cf_role not in wanted_roles:
+            my_mesh.remove_connectivity(conn)
+
+Here is a simple :class:`~iris.cube.Cube` with :class:`~iris.mesh.MeshCoord` s:
+
+>>> print(my_mesh)
+synthetic / (1)                     (-- : 96)
+    Mesh coordinates:
+        latitude                        x
+        longitude                       x
+    Mesh:
+        name                        Topology data of 2D unstructured mesh
+        location                    face
+    Attributes:
+        NCO                         'netCDF Operators version 4.7.5 (Homepage = http://nco.sf.net, Code = h ...'
+        history                     'Mon Apr 12 01:44:41 2021: ncap2 -s synthetic=float(synthetic) mesh_C4_synthetic.nc ...'
+        nco_openmp_thread_number    1
+>>> print(my_mesh_cube.aux_coords)
+(<MeshCoord: latitude / (degrees)  mesh(Topology data of 2D unstructured mesh) location(face)  [...]+bounds  shape(96,)>, <MeshCoord: longitude / (degrees)  mesh(Topology data of 2D unstructured mesh) location(face)  [...]+bounds  shape(96,)>)
+
+Here is the default indexing behaviour:
+
+>>> from iris.experimental import mesh_coord_indexing
+>>> print(mesh_coord_indexing.SETTING.value)
+Options.AUX_COORD
+>>> indexed_cube = my_mesh_cube[:36]
+>>> print(indexed_cube.aux_coords)
+(<AuxCoord: latitude / (degrees)  [29.281, 33.301, ..., 33.301, 29.281]+bounds  shape(36,)>, <AuxCoord: longitude / (degrees)  [325.894, 348.621, ..., 191.379, 214.106]+bounds  shape(36,)>)
+
+Set the indexing mode to return a new mesh:
+
+>>> mesh_coord_indexing.SETTING.value = mesh_coord_indexing.Options.NEW_MESH
+>>> indexed_cube = my_mesh_cube[:36]
+>>> print(indexed_cube.aux_coords)
+(<MeshCoord: latitude / (degrees)  mesh(Topology data of 2D unstructured mesh) location(face)  [...]+bounds  shape(36,)>, <MeshCoord: longitude / (degrees)  mesh(Topology data of 2D unstructured mesh) location(face)  [...]+bounds  shape(36,)>)
+
+Or set via a context manager:
+
+>>> with mesh_coord_indexing.SETTING.context(mesh_coord_indexing.Options.AUX_COORD):
+...    indexed_cube = my_mesh_cube[:36]
+...    print(indexed_cube.aux_coords)
+(<AuxCoord: latitude / (degrees)  [29.281, 33.301, ..., 33.301, 29.281]+bounds  shape(36,)>, <AuxCoord: longitude / (degrees)  [325.894, 348.621, ..., 191.379, 214.106]+bounds  shape(36,)>)
+
+"""
+
+from contextlib import contextmanager
+import enum
+import threading
+
+
+class Options(enum.Enum):
+    """Options for what is returned when a :class:`~iris.mesh.MeshCoord` is indexed.
+
+    See the module docstring for usage instructions:
+    :mod:`~iris.experimental.mesh_coord_indexing`.
+    """
+
+    AUX_COORD = enum.auto()
+    """The default. Convert the ``MeshCoord`` to a
+    :class:`~iris.coords.AuxCoord` and index that AuxCoord.
+    """
+
+    NEW_MESH = enum.auto()
+    """Index the :attr:`~iris.mesh.MeshCoord.mesh` of the ``MeshCoord`` to
+    produce a new :class:`~iris.mesh.MeshXY` instance, then return a new
+    :class:`~iris.mesh.MeshCoord` instance based on that new mesh.
+    """
+
+    MESH_INDEX_SET = enum.auto()
+    """**EXPERIMENTAL.** Produce a :class:`iris.mesh.components._MeshIndexSet`
+    instance that references the original :class:`~iris.mesh.MeshXY` instance,
+    then return a new :class:`~iris.mesh.MeshCoord` instance based on that new
+    index set. :class:`~Iris.mesh.components._MeshIndexSet` is a read-only
+    indexed 'view' onto its original :class:`~iris.mesh.MeshXY`; behaviour of
+    this class may change from release to release while the design is
+    finalised.
+    """
+
+
+class _Setting(threading.local):
+    """Setting for what is returned when a :class:`~iris.mesh.MeshCoord` is indexed.
+
+    See the module docstring for usage instructions:
+    :mod:`~iris.experimental.mesh_coord_indexing`.
+    """
+
+    def __init__(self):
+        self._value = Options.AUX_COORD
+
+    @property
+    def value(self):
+        return self._value
+
+    @value.setter
+    def value(self, value):
+        self._value = Options(value)
+
+    @contextmanager
+    def context(self, value):
+        new_value = Options(value)
+        original_value = self._value
+        try:
+            self._value = new_value
+            yield
+        finally:
+            self._value = original_value
+
+
+SETTING = _Setting()
+"""
+Run-time setting for alternative modes of indexing a
+:class:`~iris.mesh.MeshCoord`. See the module docstring for usage
+instructions: :mod:`~iris.experimental.mesh_coord_indexing`.
+"""

lib/iris/experimental/raster.py

@@ -19,6 +19,7 @@
 import iris
 from iris._deprecation import warn_deprecated
 import iris.coord_systems
+from iris.exceptions import MonotonicityError
 
 wmsg = (
     "iris.experimental.raster is deprecated since version 3.2, and will be "
@@ -174,7 +175,7 @@ def export_geotiff(cube, fname):
             raise ValueError(msg)
 
     if coord_x.points[0] > coord_x.points[-1]:
-        raise ValueError(
+        raise MonotonicityError(
             "Coordinate {!r} x-points must be monotonically increasing.".format(name)
         )
 

lib/iris/fileformats/netcdf/saver.py

@@ -945,7 +945,6 @@ def _add_aux_coords(self, cube, cf_var_cube, dimension_names):
             MeshEdgeCoords,
             MeshFaceCoords,
             MeshNodeCoords,
-            MeshXY,
         )
 
         # Exclude any mesh coords, which are bundled in with the aux-coords.
@@ -954,12 +953,11 @@ def _add_aux_coords(self, cube, cf_var_cube, dimension_names):
         ]
 
         # Include any relevant mesh location coordinates.
-        mesh: MeshXY | None = getattr(cube, "mesh")
-        mesh_location: str | None = getattr(cube, "location")
+        mesh = getattr(cube, "mesh")
+        mesh_location = getattr(cube, "location")
+        coords_types = typing.Union[MeshNodeCoords, MeshEdgeCoords, MeshFaceCoords]
         if mesh and mesh_location:
-            location_coords: MeshNodeCoords | MeshEdgeCoords | MeshFaceCoords = getattr(
-                mesh, f"{mesh_location}_coords"
-            )
+            location_coords: coords_types = getattr(mesh, f"{mesh_location}_coords")  # type: ignore[annotation-unchecked]
             coords_to_add.extend(list(location_coords))
 
         return self._add_inner_related_vars(