_MeshIndexSet first pass (#7149)

docs/src/conf.py

@@ -252,7 +252,7 @@ def _dotv(version):
     "members": True,
     "member-order": "alphabetical",
     "undoc-members": True,
-    "private-members": False,
+    "private-members": "_MeshIndexSet",
     "special-members": False,
     # Enums are most valuable when documented as concisely as possible.
     "inherited-members": "Enum,IntEnum,ReprEnum,StrEnum",

lib/iris/analysis/_interpolation.py

@@ -14,6 +14,7 @@
 
 from iris._lazy_data import map_complete_blocks
 from iris.coords import AuxCoord, DimCoord
+from iris.exceptions import MonotonicityError
 import iris.util
 
 _DEFAULT_DTYPE = np.float16
@@ -478,7 +479,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 _points(self, sample_points, data, data_dims=None):
         """Interpolate at the specified points.

lib/iris/common/metadata.py

@@ -1594,6 +1594,114 @@ def equal(self, other, lenient=None):
         return super().equal(other, lenient=lenient)
 
 
+class MeshIndexSetMetadata(BaseMetadata):
+    """Metadata container for a :class:`~iris.mesh.components._MeshIndexSet`."""
+
+    _members = ("mesh", "location", "start_index")
+
+    __slots__ = ()
+
+    @wraps(BaseMetadata.__eq__, assigned=("__doc__",), updated=())
+    @lenient_service
+    def __eq__(self, other):
+        return super().__eq__(other)
+
+    def _combine_lenient(self, other):
+        """Perform lenient combination of metadata members for _MeshIndexSet.
+
+        Parameters
+        ----------
+        other : MeshIndexSetMetadata
+            The other metadata participating in the lenient combination.
+
+        Returns
+        -------
+        A list of combined metadata member values.
+
+        """
+
+        # It is actually "strict" : return None except where members are equal.
+        def func(field):
+            left = getattr(self, field)
+            right = getattr(other, field)
+            return left if left == right else None
+
+        # Note that, we use "_members" not "_fields".
+        values = [func(field) for field in self._members]
+        # Perform lenient combination of the other parent members.
+        result = super()._combine_lenient(other)
+        result.extend(values)
+
+        return result
+
+    def _compare_lenient(self, other):
+        """Perform lenient equality of metadata members for _MeshIndexSet.
+
+        Parameters
+        ----------
+        other : MeshIndexSetMetadata
+            The other metadata participating in the lenient comparison.
+
+        Returns
+        -------
+        bool
+
+        """
+        # Perform "strict" comparison for the _MeshIndexSet specific members
+        # 'mesh', 'location', 'start_index' : for equality, they must all match.
+        result = all(
+            [getattr(self, field) == getattr(other, field) for field in self._members]
+        )
+        if result:
+            # Perform lenient comparison of the other parent members.
+            result = super()._compare_lenient(other)
+
+        return result
+
+    def _difference_lenient(self, other):
+        """Perform lenient difference of metadata members for _MeshIndexSet.
+
+        Parameters
+        ----------
+        other : MeshIndexSetMetadata
+            The other metadata participating in the lenient difference.
+
+        Returns
+        -------
+        A list of different metadata member values.
+
+        """
+
+        # Perform "strict" difference for location / axis.
+        def func(field):
+            left = getattr(self, field)
+            right = getattr(other, field)
+            return None if left == right else (left, right)
+
+        # Note that, we use "_members" not "_fields".
+        values = [func(field) for field in self._members]
+        # Perform lenient difference of the other parent members.
+        result = super()._difference_lenient(other)
+        result.extend(values)
+
+        return result
+
+    @wraps(BaseMetadata.combine, assigned=("__doc__",), updated=())
+    @lenient_service
+    def combine(self, other, lenient=None):
+        return super().combine(other, lenient=lenient)
+
+    @wraps(BaseMetadata.difference, assigned=("__doc__",), updated=())
+    @lenient_service
+    def difference(self, other, lenient=None):
+        return super().difference(other, lenient=lenient)
+
+    @wraps(BaseMetadata.equal, assigned=("__doc__",), updated=())
+    @lenient_service
+    def equal(self, other, lenient=None):
+        return super().equal(other, lenient=lenient)
+
+
 class MeshCoordMetadata(BaseMetadata):
     """Metadata container for a :class:`~iris.coords.MeshCoord`."""
 

lib/iris/coords.py

@@ -277,7 +277,18 @@ def _sanitise_array(self, src, ndmin):
     @property
     def _values(self):
         """The _DimensionalMetadata values as a NumPy array."""
-        return self._values_dm.data.view()
+
+        def data_id():
+            return id(self._values_dm.core_data())
+
+        original_id = data_id()
+        result = self._values_dm.data.view()
+        if data_id() != original_id:
+            # Realisation has occurred - potential effect on deferred mesh computations
+            #  (MeshCoord, _MeshIndexSet).
+            for timestamp in self._mesh_timestamps:
+                timestamp.update()
+        return result
 
     @_values.setter
     def _values(self, values):
@@ -2425,7 +2436,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()
             )
 
@@ -2999,7 +3010,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):
@@ -3079,7 +3092,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)
@@ -3089,7 +3102,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

@@ -26,7 +26,7 @@
 from functools import partial, reduce
 import itertools
 import operator
-from typing import TYPE_CHECKING, Any, Optional, TypeAlias, TypeGuard
+from typing import TYPE_CHECKING, Any, Literal, Optional, TypeAlias, TypeGuard
 import warnings
 from xml.dom.minidom import Document
 
@@ -3289,8 +3289,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]
             aux_coords.append((new_coord, new_coord_dims(coord)))
@@ -3315,8 +3314,7 @@ def new_ancillary_variable_dims(av_):
                     else:
                         dim_coords.append((new_coord, new_dims))
                         shape += new_coord.core_points().shape
-                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]
                     aux_coords.append((new_coord, new_dims))

lib/iris/exceptions.py

@@ -185,3 +185,9 @@ class CFParseError(IrisError):
     """Raised when a string associated with a CF defined syntax could not be parsed."""
 
     pass
+
+
+class MonotonicityError(ValueError):
+    """Raised when a coordinate values are not monotonic."""
+
+    pass

lib/iris/experimental/mesh_coord_indexing.py

@@ -0,0 +1,151 @@
+# 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`.
+
+.. z_reference:: iris.experimental.mesh_coord_indexing
+   :tags: topic_mesh
+
+   API reference
+
+:class:`Options` describes the available indexing modes.
+
+Select the desired option using the run-time setting :data:`SETTING`.
+
+Examples
+--------
+.. testsetup::
+
+    import iris
+    from iris.experimental.mesh_coord_indexing import SETTING
+    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_connectivities(conn)
+
+    # Capture original state.
+    original_setting = SETTING.value
+
+.. testcleanup::
+
+   # Restore original state.
+   SETTING.value = original_setting
+
+Here is a simple :class:`~iris.cube.Cube` with :class:`~iris.mesh.MeshCoord` s:
+
+>>> print(my_mesh_cube)
+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(<MeshXY object at ...>) location(face)  [...]+bounds  shape(36,)>, <MeshCoord: longitude / (degrees)  mesh(<MeshXY object at ...>) 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
+
+# TODO: update the full documentation
+
+
+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`.
+"""