diff --git a/lib/iris/analysis/_interpolation.py b/lib/iris/analysis/_interpolation.py index 6904c5ae4f..1aafe4e041 100644 --- a/lib/iris/analysis/_interpolation.py +++ b/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.""" diff --git a/lib/iris/coords.py b/lib/iris/coords.py index d2f5b05f89..a7373aca01 100644 --- a/lib/iris/coords.py +++ b/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 diff --git a/lib/iris/cube.py b/lib/iris/cube.py index 54e086937d..2f6a460683 100644 --- a/lib/iris/cube.py +++ b/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) diff --git a/lib/iris/exceptions.py b/lib/iris/exceptions.py index d6d2084d3c..cb3f426bbc 100644 --- a/lib/iris/exceptions.py +++ b/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 diff --git a/lib/iris/experimental/__init__.py b/lib/iris/experimental/__init__.py index eea4259355..dd4da95505 100644 --- a/lib/iris/experimental/__init__.py +++ b/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. +# """ diff --git a/lib/iris/experimental/mesh_coord_indexing.py b/lib/iris/experimental/mesh_coord_indexing.py new file mode 100644 index 0000000000..0e19d04fb3 --- /dev/null +++ b/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) +(, ) + +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) +(, ) + +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) +(, ) + +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) +(, ) + +""" + +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`. +""" diff --git a/lib/iris/experimental/raster.py b/lib/iris/experimental/raster.py index 52ef2f651b..b4975d7a8c 100644 --- a/lib/iris/experimental/raster.py +++ b/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) ) diff --git a/lib/iris/fileformats/netcdf/saver.py b/lib/iris/fileformats/netcdf/saver.py index cfc69143ae..481386749f 100644 --- a/lib/iris/fileformats/netcdf/saver.py +++ b/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( diff --git a/lib/iris/mesh/components.py b/lib/iris/mesh/components.py index a5936388f8..0a2dc6ede3 100644 --- a/lib/iris/mesh/components.py +++ b/lib/iris/mesh/components.py @@ -13,14 +13,17 @@ from collections import namedtuple from collections.abc import Container from contextlib import contextmanager -from typing import Iterable +from copy import copy, deepcopy +from typing import Iterable, Literal import warnings from cf_units import Unit from dask import array as da import numpy as np +from numpy.typing import ArrayLike from iris.common.metadata import ConnectivityMetadata, MeshCoordMetadata, MeshMetadata +import iris.util from .. import _lazy_data as _lazy from ..common import CFVariableMixin, metadata_filter, metadata_manager_factory @@ -597,6 +600,7 @@ class Mesh(CFVariableMixin, ABC): - Move whatever is appropriate from :class:`MeshXY` into this class, leaving behind only those elements specific to the assumption of X and Y node coordinates. + - Set :class:`_MeshIndexSet` to subclass `Mesh` instead of `MeshXY`. - Remove the docstring warning, the NotImplementedError, and the uses of ABC/abstractmethod. - Add a cross-reference in the docstring for :class:`MeshXY`. @@ -727,11 +731,11 @@ def normalise(element, axis): raise ValueError(emsg) if self.topology_dimension == 1: - self._coord_manager = _Mesh1DCoordinateManager(**kwargs) - self._connectivity_manager = _Mesh1DConnectivityManager(*connectivities) + self._coord_man = _Mesh1DCoordinateManager(**kwargs) + self._connectivity_man = _Mesh1DConnectivityManager(*connectivities) elif self.topology_dimension == 2: - self._coord_manager = _Mesh2DCoordinateManager(**kwargs) - self._connectivity_manager = _Mesh2DConnectivityManager(*connectivities) + self._coord_man = _Mesh2DCoordinateManager(**kwargs) + self._connectivity_man = _Mesh2DConnectivityManager(*connectivities) else: emsg = f"Unsupported 'topology_dimension', got {topology_dimension!r}." raise NotImplementedError(emsg) @@ -956,8 +960,8 @@ def __hash__(self): def __getstate__(self): return ( self._metadata_manager, - self._coord_manager, - self._connectivity_manager, + self._coord_man, + self._connectivity_man, ) def __ne__(self, other): @@ -1002,11 +1006,11 @@ def _summary_oneline(self): mesh_name = None if mesh_name: # Use a more human-readable form - mesh_string = f"" + mesh_string = f"<{self.__class__.__name__}: '{mesh_name}'>" else: # Mimic the generic object.__str__ style. mesh_id = id(self) - mesh_string = f"" + mesh_string = f"<{self.__class__.__name__} object at {hex(mesh_id)}>" return mesh_string @@ -1020,7 +1024,7 @@ def line(text, i_indent=0): indent = indent_str * i_indent lines.append(f"{indent}{text}") - line(f"MeshXY : '{self.name()}'") + line(f"{self.__class__.__name__} : '{self.name()}'") line(f"topology_dimension: {self.topology_dimension}", 1) for element in ("node", "edge", "face"): if element == "node": @@ -1101,8 +1105,8 @@ def line(text, i_indent=0): def __setstate__(self, state): metadata_manager, coord_manager, connectivity_manager = state self._metadata_manager = metadata_manager - self._coord_manager = coord_manager - self._connectivity_manager = connectivity_manager + self._coord_man = coord_manager + self._connectivity_man = connectivity_manager def _set_dimension_names(self, node, edge, face, reset=False): args = (node, edge, face) @@ -1133,6 +1137,16 @@ def _set_dimension_names(self, node, edge, face, reset=False): return result + @property + def _connectivity_manager(self): + # Delivered as a property to enable more sophisticated overriding. + return self._connectivity_man + + @property + def _coord_manager(self): + # Delivered as a property to enable more sophisticated overriding. + return self._coord_man + @property def all_connectivities(self): """All the :class:`~iris.mesh.Connectivity` instances of the :class:`MeshXY`.""" @@ -1956,6 +1970,24 @@ def dimension_names(self, node=None, edge=None, face=None): """ return self._set_dimension_names(node, edge, face, reset=False) + def is_view_of(self, other: "MeshXY") -> bool: + """Whether this instance is either itself, or a view of the given :class:`MeshXY`. + + Parameters + ---------- + other : iris.mesh.components.MeshXY + The :class:`~MeshXY` to compare against. + + Returns + ------- + bool + """ + # The Mesh parent class does not implement any operations for returning + # alternative views of its content, so the only possible 'view' of + # itself IS itself. + # Subclasses may be more sophisticated in this regard. + return other is self + @property def cf_role(self): """The UGRID ``cf_role`` attribute of the :class:`MeshXY`.""" @@ -1973,6 +2005,278 @@ def topology_dimension(self): return self._metadata_manager.topology_dimension +class _MeshIndexSet(MeshXY): + """A container representing the UGRID ``cf_role``: ``location_index_set``. + + A container representing the UGRID [1]_ ``cf_role``: + ``location_index_set``. Achieved by referencing an original :class:`MeshXY` + instance (:attr:`super_mesh`), together with a specific :attr:`location` + (``node``/``edge``/``face``) and a set of :attr:`indices`. This is strictly + a view onto the original :class:`MeshXY` - it does not store its own + coordinates or connectivities. + + Instances should be created via the :meth:`from_mesh` class method, rather + than directly via ``__init__()``. + + Warnings + -------- + This class is not yet stable, hence being marked as private. + + References + ---------- + .. [1] The UGRID Conventions, https://ugrid-conventions.github.io/ugrid-conventions/ + """ + + # TODO: make the metadata manager a read-only view of the original Mesh. + # Or maybe disconnect, since I/O will allow different long_names etc. + # TODO: implement I/O (iris#6123). + # TODO: cf_role + + super_mesh: MeshXY | None = None + """The :class:`MeshXY` instance that this is a view onto.""" + + location: Literal["node", "edge", "face", None] = None + """The location on :attr:`super_mesh` that is being indexed.""" + + indices: ArrayLike | None = None + """Which elements of the :attr:`location` of :attr:`super_mesh` belong.""" + + class _IndexedMembers(dict): + # Supports the read-only nature of _connectivity_manager and _coord_manager. + # Replaces the `_members` attribute of the managers. + + def _readonly(self, *args, **kwargs): + message = ( + "Modification of _MeshIndexSet is forbidden - this is only " + f"a view onto an original MeshXY: id={self.mesh_id}." + ) + raise RuntimeError(message) + + __setitem__ = _readonly + __delitem__ = _readonly + pop = _readonly # type: ignore[assignment] + popitem = _readonly + clear = _readonly + update = _readonly # type: ignore[assignment] + setdefault = _readonly # type: ignore[assignment] + + mesh_id: int + + def __init__(self, seq, **kwargs): + self.mesh_id = kwargs.pop("mesh_id") + super().__init__(seq, **kwargs) + + def __getstate__(self): + self._validate() + state_inherited = super().__getstate__() + return ( + self.super_mesh, + self.location, + self.indices, + state_inherited, + ) + + def __setstate__(self, state): + self.super_mesh, self.location, self.indices, state_inherited = state + super().__setstate__(state_inherited) + self._validate() + + @classmethod + def from_mesh( + cls, mesh: MeshXY, location: Literal["node", "edge", "face"], indices: ArrayLike + ) -> "_MeshIndexSet": + """Create a new instance by indexing a :class:`MeshXY` instance. + + Parameters + ---------- + mesh : MeshXY + The :class:`MeshXY` to index. The created instance will be a view + onto `mesh`. + location : {"node", "edge", "face"} + The location on `mesh` that is being indexed. + indices : ArrayLike + Which elements of the `location` of `mesh` belong. + + Returns + ------- + _MeshIndexSet + + """ + if not hasattr(mesh, "topology_dimension"): + message = ( + "_MeshIndexSet.from_mesh() requires a MeshXY instance, got: " + f"{type(mesh).__name__}" + ) + raise ValueError(message) + else: + instance_copied = copy(mesh) + instance_copied.__class__ = cls + instance: _MeshIndexSet = instance_copied # type: ignore[assignment] + instance.super_mesh = mesh + instance.location = location + instance.indices = indices + return instance + + def _validate(self) -> None: + # Instances are expected to be created via the from_mesh() method. + # If __init__() is used instead, the class will not be fully set up. + # Actually disabling __init__() is considered an anti-pattern. + if any(v is None for v in [self.super_mesh, self.location, self.indices]): + message = ( + "_MeshIndexSet instance not fully set up. Note that " + "_MeshIndexSet is designed to be created via the from_mesh() " + "method, not via __init__()." + ) + raise RuntimeError(message) + + # The 3 managers are part of MeshXY.__getstate__(), so the copy() op + # in from_mesh() should produce identical references. + if not all( + getattr(self, man) is getattr(self.super_mesh, man) + for man in ["_metadata_manager", "_coord_man", "_connectivity_man"] + ): + message = ( + f"_MeshIndexSet (id={id(self)}) is not a valid view of " + f"its original MeshXY (id={id(self.super_mesh)})." + ) + raise RuntimeError(message) + + def _calculate_node_indices(self): + # Use self.location and self.indices to work out the indices to use + # when indexing the nodes of self.super_mesh. + self._validate() + if self.location == "node": + result = self.indices + elif self.location in ["edge", "face"]: + (connectivity,) = [ + c + for c in self.super_mesh.all_connectivities + if ( + c is not None + and c.location == self.location + and c.connected == "node" + ) + ] + # Doesn't matter if connectivity is transposed or not in this case. + conn_indices = connectivity.indices[self.indices] + node_set = np.unique(conn_indices) + if iris.util.is_masked(node_set): + node_set_unmasked = node_set.compressed() + else: + node_set_unmasked = node_set + result = node_set_unmasked + else: + result = None + # TODO: should this be validated earlier? + # Maybe even with an Enum? + message = ( + f"Expected location to be one of `node`, `edge` or `face`, " + f"got `{self.location}`" + ) + raise NotImplementedError(message) + + return result + + def _calculate_edge_indices(self): + # Use self.location and self.indices to work out the indices to use + # when indexing the edges of self.super_mesh. + self._validate() + if self.location == "edge": + result = self.indices + else: + result = None + return result + + def _calculate_face_indices(self): + # Use self.location and self.indices to work out the indices to use + # when indexing the faces of self.super_mesh. + self._validate() + if self.location == "face": + result = self.indices + else: + result = None + return result + + @property + def _coord_manager(self): + # Intended to be a 'view' on the original, and Meshes are mutable, so + # must re-index every time it is accessed. + # TODO: only re-index when the managers have changed (also being + # considered in iris#4757). + return self._coord_man.indexed( + self._calculate_node_indices(), + self._calculate_edge_indices(), + self._calculate_face_indices(), + mesh_id=id(self.super_mesh), + ) + + @property + def _connectivity_manager(self): + # Intended to be a 'view' on the original, and Meshes are mutable, so + # must re-index every time it is accessed. + # TODO: only re-index when the managers have changed (also being + # considered in iris#4757). + return self._connectivity_man.indexed( + self._calculate_node_indices(), + self._calculate_edge_indices(), + self._calculate_face_indices(), + mesh_id=id(self.super_mesh), + ) + + def is_view_of(self, other: MeshXY) -> bool: + """Whether this instance is either itself, or a view of the given :class:`MeshXY`. + + Parameters + ---------- + other : MeshXY + The :class:`MeshXY` to compare against. + + Returns + ------- + bool + """ + return other is self.super_mesh or other is self + + def convert_to_mesh(self) -> None: + """Convert this instance **in-place** to be a :class:`MeshXY`. + + This will no longer be a view onto another :class:`MeshXY`. All + connectivities and coordinates will be deep-copied. + """ + final_coords = self._coord_manager + final_coords.force_mutability() + final_coords = deepcopy(final_coords) + + final_connectivities = self._connectivity_manager + final_connectivities.force_mutability() + final_connectivities = deepcopy(final_connectivities) + + self._metadata_manager = deepcopy(self._metadata_manager) + self._coord_man = final_coords + self._connectivity_man = final_connectivities + + self.__class__ = MeshXY # type: ignore[assignment] + del self.super_mesh + del self.location + del self.indices + + def as_mesh(self) -> MeshXY: + """Return a :class:`MeshXY` representation of this instance. + + The returned instance will no longer be a view onto another + :class:`MeshXY`. All connectivities and coordinates will be deep-copied. + + Returns + ------- + MeshXY + """ + result = copy(self) + # Using `convert_to_mesh()` avoids modifying private attributes of + # `result`. + result.convert_to_mesh() + return result + + class _Mesh1DCoordinateManager: """TBD: require clarity on coord_systems validation. @@ -2098,7 +2402,7 @@ def all_members(self): return Mesh1DCoords(**self._members) @property - def edge_coords(self): + def edge_coords(self) -> MeshEdgeCoords: return MeshEdgeCoords(edge_x=self.edge_x, edge_y=self.edge_y) @property @@ -2118,7 +2422,7 @@ def edge_y(self, coord): self._setter(element="edge", axis="y", coord=coord, shape=self._edge_shape) @property - def node_coords(self): + def node_coords(self) -> MeshNodeCoords: return MeshNodeCoords(node_x=self.node_x, node_y=self.node_y) @property @@ -2257,6 +2561,87 @@ def populated_coords(coords_tuple): result_dict = {k: v for k, v in self._members.items() if id(v) in result_ids} return result_dict + def index( + self, + node_indices: ArrayLike, + edge_indices: ArrayLike, + face_indices: ArrayLike, + mesh_id: int, + ) -> None: + """Permanently index the members of this coordinate manager **in-place**. + + Members become read-only - this method is intended to support the creation + of indexed views of original meshes. See :meth:`force_mutability` for + undoing the read-only nature. + + Parameters + ---------- + node_indices, edge_indices, face_indices : ArrayLike + The indices to use when indexing member node, edge or face + coordinates respectively. + mesh_id : int + The ID of the mesh that these indices refer to - used to + produce a meaningful read-only error. + """ + indices_dict = { + "node": node_indices, + "edge": edge_indices, + "face": face_indices, + } + members_indexed = {} + for key, coord in self._members.items(): + indexing = None + indexed = None + if coord is not None: + indexing = indices_dict[key.split("_")[0]] + if indexing is not None: + indexed = coord[indexing] + members_indexed[key] = indexed + + self._members = _MeshIndexSet._IndexedMembers( + members_indexed, + mesh_id=mesh_id, + ) + + def indexed( + self, + node_indices: ArrayLike, + edge_indices: ArrayLike, + face_indices: ArrayLike, + mesh_id: int, + ) -> "_Mesh1DCoordinateManager": + """Return an indexed copy of this coordinate manager. + + Members are read-only - this method is intended to support the creation + of indexed views of original meshes. See :meth:`force_mutability` for + undoing the read-only nature. + + Parameters + ---------- + node_indices, edge_indices, face_indices : ArrayLike + The indices to use when indexing member node, edge or face + coordinates respectively. + mesh_id : int + The ID of the mesh that these indices refer to - used to + produce a meaningful read-only error. + + Returns + ------- + _Mesh1DCoordinateManager + """ + result = copy(self) + # Using `index()` avoids modifying private attributes of `result`. + result.index(node_indices, edge_indices, face_indices, mesh_id) + return result + + def force_mutability(self) -> None: + """Permanently allow modification of this instance's members **in-place**. + + Intended to support the conversion of mesh views to be independent meshes. + """ + if hasattr(self._members, "mesh_id"): + self._members = dict(self._members) + def remove( self, item=None, @@ -2310,7 +2695,7 @@ def all_members(self): return Mesh2DCoords(**self._members) @property - def face_coords(self): + def face_coords(self) -> MeshFaceCoords: return MeshFaceCoords(face_x=self.face_x, face_y=self.face_y) @property @@ -2552,6 +2937,96 @@ def element_filter(instances, loc_arg, loc_name): result_dict = {k: v for k, v in self._members.items() if id(v) in result_ids} return result_dict + def index( + self, + node_indices: ArrayLike, + edge_indices: ArrayLike, + face_indices: ArrayLike, + mesh_id: int, + ) -> None: + """Permanently index the members of this connectivity manager **in-place**. + + Members become read-only - this method is intended to support the creation + of indexed views of original meshes. See :meth:`force_mutability` for + undoing the read-only nature. + + Parameters + ---------- + node_indices, edge_indices, face_indices : ArrayLike + The indices to use when indexing member connectivities with node, + edge or face :attr:`~Connectivity.location` respectively. + mesh_id : int + The ID of the mesh that these indices refer to - used to + produce a meaningful read-only error. + """ + indices_dict = { + "node": node_indices, + "edge": edge_indices, + "face": face_indices, + } + node_index_mapping = { + old_index: new_index for new_index, old_index in enumerate(node_indices) + } + members_indexed = {} + for key, connectivity in self._members.items(): + indexing = None + indexed = None + if connectivity is not None: + indexing = indices_dict[connectivity.location] + if indexing is not None: + new_values = connectivity.indices_by_location()[indexing] + new_values = np.vectorize(node_index_mapping.get)(new_values) + if connectivity.location_axis == 1: + new_values = new_values.T + if connectivity.start_index == 1: + new_values = new_values + 1 + indexed = connectivity.copy(new_values) + members_indexed[key] = indexed + + self._members = _MeshIndexSet._IndexedMembers( + members_indexed, + mesh_id=mesh_id, + ) + + def indexed( + self, + node_indices: ArrayLike, + edge_indices: ArrayLike, + face_indices: ArrayLike, + mesh_id: int, + ) -> "_MeshConnectivityManagerBase": + """Return an indexed copy of this connectivity manager. + + Members are read-only - this method is intended to support the creation + of indexed views of original meshes. See :meth:`force_mutability` for + undoing the read-only nature. + + Parameters + ---------- + node_indices, edge_indices, face_indices : ArrayLike + The indices to use when indexing member connectivities with node, + edge or face :attr:`~Connectivity.location` respectively. + mesh_id : int + The ID of the mesh that these indices refer to - used to + produce a meaningful read-only error. + + Returns + ------- + _MeshConnectivityManagerBase + """ + result = copy(self) + # Using `index()` avoids modifying private attributes of `result`. + result.index(node_indices, edge_indices, face_indices, mesh_id) + return result + + def force_mutability(self) -> None: + """Permanently allow modification of this instance's members **in-place**. + + Intended to support the conversion of mesh views to be independent meshes. + """ + if hasattr(self._members, "mesh_id"): + self._members = dict(self._members) + def remove( self, item=None, @@ -2852,17 +3327,55 @@ def climatological(self, value): raise ValueError(msg) def __getitem__(self, keys): - # Disallow any sub-indexing, permitting *only* "self[:,]". - # We *don't* intend here to support indexing as such : the exception is - # just sufficient to enable cube slicing, when it does not affect the - # mesh dimension. This works because Cube.__getitem__ passes us keys - # "normalised" with iris.util._build_full_slice_given_keys. - if keys != (slice(None),): - msg = "Cannot index a MeshCoord." - raise ValueError(msg) + # TODO: validate keys as 1-dimensional? + # TODO: handle problems caused by asking for 1 index + # E.g. coord[0:1] works fine, coord[0] causes errors difficult to + # understand. + from ..experimental.mesh_coord_indexing import SETTING, Options + + def get_index_set(): + return _MeshIndexSet.from_mesh(self.mesh, self.location, keys) + + match SETTING.value: + case Options.AUX_COORD: + result = AuxCoord.from_coord(self)[keys] + + case Options.NEW_MESH: + new_mesh = get_index_set().as_mesh() + result = self.__class__( + mesh=new_mesh, + location=self.location, + axis=self.axis, + ) - # Translate "self[:,]" as "self.copy()". - return self.copy() + case Options.MESH_INDEX_SET: + index_set = get_index_set() + result = self.__class__( + mesh=index_set, + location=self.location, + axis=self.axis, + ) + + case _: + message = ( + "Unsupported mesh_coord_indexing setting. Expected one of: " + f"{Options.__members__}, got: {SETTING.value}." + ) + raise NotImplementedError(message) + + return result + + # # Disallow any sub-indexing, permitting *only* "self[:,]". + # # We *don't* intend here to support indexing as such : the exception is + # # just sufficient to enable cube slicing, when it does not affect the + # # mesh dimension. This works because Cube.__getitem__ passes us keys + # # "normalised" with iris.util._build_full_slice_given_keys. + # if keys != (slice(None),): + # msg = "Cannot index a MeshCoord." + # raise ValueError(msg) + # + # # Translate "self[:,]" as "self.copy()". + # return self.copy() def collapsed(self, dims_to_collapse=None): """Return a copy of this coordinate, which has been collapsed along the specified dimensions. diff --git a/lib/iris/pandas.py b/lib/iris/pandas.py index dd0f9fec89..800193f51a 100644 --- a/lib/iris/pandas.py +++ b/lib/iris/pandas.py @@ -28,6 +28,7 @@ from iris._deprecation import warn_deprecated from iris.coords import AncillaryVariable, AuxCoord, CellMeasure, DimCoord from iris.cube import Cube, CubeList +from iris.exceptions import MonotonicityError from iris.warnings import IrisIgnoringWarning @@ -383,7 +384,7 @@ def as_cubes( "Pandas index is not monotonic. Consider using the " "sort_index() method before passing in." ) - raise ValueError(message) + raise MonotonicityError(message) cube_shape = getattr(pandas_index, "levshape", (pandas_index.nunique(),)) n_rows = len(pandas_structure) diff --git a/lib/iris/tests/unit/cube/test_Cube.py b/lib/iris/tests/unit/cube/test_Cube.py index 878a793448..6c661635a5 100644 --- a/lib/iris/tests/unit/cube/test_Cube.py +++ b/lib/iris/tests/unit/cube/test_Cube.py @@ -2589,6 +2589,39 @@ def test(self): self.assertEqual((aux / cube).units, "m.s-1") +class Test__getitem_DimCoord: + @pytest.fixture(autouse=True) + def _setup(self): + cube = Cube(np.arange(6).reshape(2, 3)) + x_coord = DimCoord(points=np.array([2, 3, 4]), long_name="x") + cube.add_dim_coord(x_coord, 1) + x_coord_aux = DimCoord(points=np.array([2, 3, 4]), long_name="x_aux") + cube.add_aux_coord(x_coord_aux, 1) + y_coord = DimCoord(points=np.array([5, 6]), long_name="y") + cube.add_dim_coord(y_coord, 0) + z_coord = AuxCoord(points=np.arange(6).reshape(2, 3), long_name="z") + cube.add_aux_coord(z_coord, [0, 1]) + self.cube = cube + + def test_dim_coords_2d(self): + result = self.cube[0:2, 0:2] + assert len(result.dim_coords) == 2 + for ix, size in enumerate(result.shape): + assert result.coord(dim_coords=True, dimensions=ix).shape == (size,) + assert result.coord("x_aux").shape == (result.shape[0],) + + def test_dim_coord_1d(self): + result = self.cube[0, 0:2] + assert len(result.dim_coords) == 1 + assert result.dim_coords[0].shape == result.shape + assert result.coord("x_aux").shape == (result.shape[0],) + + def test_demote(self): + result = self.cube[:, [0, 2, 1]] + assert len(result.dim_coords) == 1 + assert result.coords(dimensions=1, dim_coords=True) == [] + + class Test__getitem_CellMeasure(tests.IrisTest): def setUp(self): cube = Cube(np.arange(6).reshape(2, 3))