Merge pull request #830 from davidhassell/fix-cyclic-subspace

Fix bug where `cf.Field.subspace` doesn't always correctly handle   global cyclic subspaces

Author: davidhassell

Changelog.rst

@@ -33,13 +33,15 @@ version NEXTVERSION
 * Fix bug where `cf.normalize_slice` doesn't correctly
   handle certain cyclic slices
   (https://github.com/NCAS-CMS/cf-python/issues/774)
+* Fix bug where `cf.Field.subspace` doesn't always correctly handle
+  global or near-global cyclic subspaces
+  (https://github.com/NCAS-CMS/cf-python/issues/828)
 * New dependency: ``h5netcdf>=1.3.0``
 * New dependency: ``h5py>=3.10.0``
 * New dependency: ``s3fs>=2024.2.0``
 * Changed dependency: ``1.11.2.0<=cfdm<1.11.3.0``
 * Changed dependency: ``cfunits>=3.3.7``
 
-
 ----
 
 version 3.16.2

cf/cfimplementation.py

@@ -26,7 +26,6 @@
     TiePointIndex,
 )
 from .data import Data
-
 from .data.array import (
     BoundsFromNodesArray,
     CellConnectivityArray,

cf/data/data.py

@@ -43,7 +43,6 @@
 from ..units import Units
 from .collapse import Collapse
 from .creation import generate_axis_identifiers, to_dask
-
 from .dask_utils import (
     _da_ma_allclose,
     cf_asanyarray,

cf/dimensioncoordinate.py

@@ -8,7 +8,11 @@
     _inplace_enabled,
     _inplace_enabled_define_and_cleanup,
 )
-from .functions import _DEPRECATION_ERROR_ATTRIBUTE, _DEPRECATION_ERROR_KWARGS
+from .functions import (
+    _DEPRECATION_ERROR_ATTRIBUTE,
+    _DEPRECATION_ERROR_KWARGS,
+    bounds_combination_mode,
+)
 from .timeduration import TimeDuration
 from .units import Units
 
@@ -246,6 +250,154 @@ def increasing(self):
         """
         return self.direction()
 
+    @_inplace_enabled(default=False)
+    def anchor(self, value, cell=False, parameters=None, inplace=False):
+        """Anchor the coordinate values.
+
+        By default, the coordinate values are transformed so that the
+        first coordinate is the closest to *value* from above (below)
+        for increasing (decreasing) coordinates.
+
+        If the *cell* parameter is True, then the coordinate values
+        are transformed so that the first cell either contains
+        *value*; or is the closest to cell to *value* from above
+        (below) for increasing (decreasing) coordinates.
+
+        .. versionadded:: NEXTVERSION
+
+        .. seealso:: `period`, `roll`
+
+        :Parameters:
+
+            value: scalar array_like
+                Anchor the coordinate values for the selected cyclic
+                axis to the *value*. May be any numeric scalar object
+                that can be converted to a `Data` object (which
+                includes `numpy` and `Data` objects). If *value* has
+                units then they must be compatible with those of the
+                coordinates, otherwise it is assumed to have the same
+                units as the coordinates.
+
+                The coordinate values are transformed so the first
+                corodinate is the closest to *value* from above (for
+                increasing coordinates), or the closest to *value* from
+                above (for decreasing coordinates)
+
+                  * Increasing coordinates with positive period, P,
+                    are transformed so that *value* lies in the
+                    half-open range (L-P, F], where F and L are the
+                    transformed first and last coordinate values,
+                    respectively.
+
+            ..
+
+                  * Decreasing coordinates with positive period, P,
+                    are transformed so that *value* lies in the
+                    half-open range (L+P, F], where F and L are the
+                    transformed first and last coordinate values,
+                    respectively.
+
+                *Parameter example:*
+                  If the original coordinates are ``0, 5, ..., 355``
+                  (evenly spaced) and the period is ``360`` then
+                  ``value=0`` implies transformed coordinates of ``0,
+                  5, ..., 355``; ``value=-12`` implies transformed
+                  coordinates of ``-10, -5, ..., 345``; ``value=380``
+                  implies transformed coordinates of ``380, 385, ...,
+                  735``.
+
+                *Parameter example:*
+                  If the original coordinates are ``355, 350, ..., 0``
+                  (evenly spaced) and the period is ``360`` then
+                  ``value=355`` implies transformed coordinates of
+                  ``355, 350, ..., 0``; ``value=0`` implies
+                  transformed coordinates of ``0, -5, ..., -355``;
+                  ``value=392`` implies transformed coordinates of
+                  ``390, 385, ..., 35``.
+
+            cell: `bool`, optional
+                If True, then the coordinate values are transformed so
+                that the first cell either contains *value*, or is the
+                closest to cell to *value* from above (below) for
+                increasing (decreasing) coordinates.
+
+                If False (the default) then the coordinate values are
+                transformed so that the first coordinate is the closest
+                to *value* from above (below) for increasing
+                (decreasing) coordinates.
+
+            parameters: `dict`, optional
+                If a `dict` is provided then it will be updated
+                in-place with parameters which describe thethe
+                anchoring process.
+
+            {{inplace: `bool`, optional}}
+
+        :Returns:
+
+            `{{class}}` or `None`
+                The anchored dimension coordinates, or `None` if the
+                operation was in-place.
+
+        """
+        d = _inplace_enabled_define_and_cleanup(self)
+
+        period = d.period()
+        if period is None:
+            raise ValueError(f"Cyclic {d!r} has no period")
+
+        value = d._Data.asdata(value)
+        if not value.Units:
+            value = value.override_units(d.Units)
+        elif not value.Units.equivalent(d.Units):
+            raise ValueError(
+                f"Anchor value has incompatible units: {value.Units!r}"
+            )
+
+        if cell:
+            c = d.upper_bounds.persist()
+        else:
+            d.persist(inplace=True)
+            c = d.get_data(_fill_value=False)
+
+        if d.increasing:
+            # Adjust value so it's in the range [c[0], c[0]+period)
+            n = ((c[0] - value) / period).ceil()
+            value1 = value + n * period
+            shift = c.size - np.argmax((c - value1 >= 0).array)
+            d.roll(0, shift, inplace=True)
+            if cell:
+                d0 = d[0].upper_bounds
+            else:
+                d0 = d.get_data(_fill_value=False)[0]
+
+            n = ((value - d0) / period).ceil()
+        else:
+            # Adjust value so it's in the range (c[0]-period, c[0]]
+            n = ((c[0] - value) / period).floor()
+            value1 = value + n * period
+            shift = c.size - np.argmax((value1 - c >= 0).array)
+            d.roll(0, shift, inplace=True)
+            if cell:
+                d0 = d[0].upper_bounds
+            else:
+                d0 = d.get_data(_fill_value=False)[0]
+
+            n = ((value - d0) / period).floor()
+
+        n.persist(inplace=True)
+        if n:
+            nperiod = n * period
+            with bounds_combination_mode("OR"):
+                d += nperiod
+        else:
+            nperiod = 0
+
+        if parameters is not None:
+            parameters.update({"shift": shift, "nperiod": nperiod})
+
+        return d
+
     def direction(self):
         """Return True if the dimension coordinate values are
         increasing, otherwise return False.

cf/mixin/fielddomain.py

@@ -18,7 +18,7 @@
     bounds_combination_mode,
     normalize_slice,
 )
-from ..query import Query
+from ..query import Query, wi
 from ..units import Units
 
 logger = logging.getLogger(__name__)
@@ -440,50 +440,41 @@ def _indices(self, config, data_axes, ancillary_mask, kwargs):
                     if debug:
                         logger.debug("  1-d CASE 2:")  # pragma: no cover
 
+                    size = item.size
                     if item.increasing:
-                        anchor0 = value.value[0]
-                        anchor1 = value.value[1]
+                        anchor = value.value[0]
                     else:
-                        anchor0 = value.value[1]
-                        anchor1 = value.value[0]
-
-                    a = self.anchor(axis, anchor0, dry_run=True)["roll"]
-                    b = self.flip(axis).anchor(axis, anchor1, dry_run=True)[
-                        "roll"
-                    ]
-
-                    size = item.size
-                    if abs(anchor1 - anchor0) >= item.period():
-                        if value.operator == "wo":
-                            set_start_stop = 0
-                        else:
-                            set_start_stop = -a
-
-                        start = set_start_stop
-                        stop = set_start_stop
-                    elif a + b == size:
-                        b = self.anchor(axis, anchor1, dry_run=True)["roll"]
-                        if (b == a and value.operator == "wo") or not (
-                            b == a or value.operator == "wo"
-                        ):
-                            set_start_stop = -a
-                        else:
-                            set_start_stop = 0
+                        anchor = value.value[1]
+
+                    item = item.persist()
+                    parameters = {}
+                    item = item.anchor(anchor, parameters=parameters)
+                    n = np.roll(np.arange(size), parameters["shift"], 0)
+                    if value.operator == "wi":
+                        n = n[item == value]
+                        if not n.size:
+                            raise ValueError(
+                                f"No indices found from: {identity}={value!r}"
+                            )
 
-                        start = set_start_stop
-                        stop = set_start_stop
+                        start = n[0]
+                        stop = n[-1] + 1
                     else:
-                        if value.operator == "wo":
-                            start = b - size
-                            stop = -a + size
-                        else:
-                            start = -a
-                            stop = b - size
+                        # "wo" operator
+                        n = n[item == wi(*value.value)]
+                        if n.size == size:
+                            raise ValueError(
+                                f"No indices found from: {identity}={value!r}"
+                            )
 
-                    if start == stop == 0:
-                        raise ValueError(
-                            f"No indices found from: {identity}={value!r}"
-                        )
+                        if n.size:
+                            start = n[-1] + 1
+                            stop = start - n.size
+                        else:
+                            start = size - parameters["shift"]
+                            stop = start + size
+                            if stop > size:
+                                stop -= size
 
                     index = slice(start, stop, 1)
 
@@ -1287,77 +1278,34 @@ def anchor(
                 self, "anchor", kwargs
             )  # pragma: no cover
 
-        da_key, axis = self.domain_axis(axis, item=True)
+        axis = self.domain_axis(axis, key=True)
 
         if dry_run:
             f = self
         else:
             f = _inplace_enabled_define_and_cleanup(self)
 
-        dim = f.dimension_coordinate(filter_by_axis=(da_key,), default=None)
+        dim = f.dimension_coordinate(filter_by_axis=(axis,), default=None)
         if dim is None:
             raise ValueError(
                 "Can't shift non-cyclic "
-                f"{f.constructs.domain_axis_identity(da_key)!r} axis"
+                f"{f.constructs.domain_axis_identity(axis)!r} axis"
             )
 
-        period = dim.period()
-        if period is None:
-            raise ValueError(f"Cyclic {dim.identity()!r} axis has no period")
-
-        value = f._Data.asdata(value)
-        if not value.Units:
-            value = value.override_units(dim.Units)
-        elif not value.Units.equivalent(dim.Units):
-            raise ValueError(
-                f"Anchor value has incompatible units: {value.Units!r}"
-            )
-
-        axis_size = axis.get_size()
-
-        if axis_size <= 1:
-            # Don't need to roll a size one axis
-            if dry_run:
-                return {"axis": da_key, "roll": 0, "nperiod": 0}
-
-            return f
-
-        c = dim.get_data(_fill_value=False)
-
-        if dim.increasing:
-            # Adjust value so it's in the range [c[0], c[0]+period)
-            n = ((c[0] - value) / period).ceil()
-            value1 = value + n * period
-
-            shift = axis_size - np.argmax((c - value1 >= 0).array)
-            if not dry_run:
-                f.roll(da_key, shift, inplace=True)
-
-            # Re-get dim
-            dim = f.dimension_coordinate(filter_by_axis=(da_key,))
-            # TODO CHECK n for dry run or not
-            n = ((value - dim.data[0]) / period).ceil()
-        else:
-            # Adjust value so it's in the range (c[0]-period, c[0]]
-            n = ((c[0] - value) / period).floor()
-            value1 = value + n * period
-
-            shift = axis_size - np.argmax((value1 - c >= 0).array)
-
-            if not dry_run:
-                f.roll(da_key, shift, inplace=True)
-
-            # Re-get dim
-            dim = f.dimension_coordinate(filter_by_axis=(da_key,))
-            # TODO CHECK n for dry run or not
-            n = ((value - dim.data[0]) / period).floor()
+        parameters = {"axis": axis}
+        dim = dim.anchor(value, parameters=parameters)
 
         if dry_run:
-            return {"axis": da_key, "roll": shift, "nperiod": n * period}
+            return parameters
+
+        f.roll(axis, parameters["shift"], inplace=True)
 
-        if n:
+        if parameters["nperiod"]:
+            # Get the rolled dimension coordinate and adjust the
+            # values by the non-zero integer multiple of 'period'
+            dim = f.dimension_coordinate(filter_by_axis=(axis,))
             with bounds_combination_mode("OR"):
-                dim += n * period
+                dim += parameters["nperiod"]
 
         return f
 

cf/regrid/regrid.py

@@ -2465,7 +2465,6 @@ def create_esmpy_weights(
             from netCDF4 import Dataset
 
             from .. import __version__
-
             from ..data.array.locks import netcdf_lock
 
             if (