Merge branch 'main' into dask-in-cfdm

Author: davidhassell

Changelog.rst

@@ -3,6 +3,9 @@ version NEXTVERSION
 
 **2024-12-??**
 
+* Allow ``'nearest_dtos'`` 2-d regridding to work with discrete
+  sampling geometry source grids
+  (https://github.com/NCAS-CMS/cf-python/issues/832)
 * New method: `cf.Field.filled`
   (https://github.com/NCAS-CMS/cf-python/issues/811)
 * New method: `cf.Field.is_discrete_axis`
@@ -35,13 +38,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/data/dask_regrid.py

@@ -506,17 +506,20 @@ def _regrid(
             # Note: It is much more efficient to access
             #       'weights.indptr', 'weights.indices', and
             #       'weights.data' directly, rather than iterating
-            #       over rows of 'weights' and using 'weights.getrow'.
+            #       over rows of 'weights' and using
+            #       'weights.getrow'. Also, 'np.count_nonzero' is much
+            #       faster than 'np.any' and 'np.all'.
             count_nonzero = np.count_nonzero
             indptr = weights.indptr.tolist()
             indices = weights.indices
             data = weights.data
             for j, (i0, i1) in enumerate(zip(indptr[:-1], indptr[1:])):
                 mask = src_mask[indices[i0:i1]]
-                if not count_nonzero(mask):
+                n_masked = count_nonzero(mask)
+                if not n_masked:
                     continue
 
-                if mask.all():
+                if n_masked == mask.size:
                     dst_mask[j] = True
                     continue
 
@@ -528,8 +531,8 @@ def _regrid(
 
             del indptr
 
-        elif method in ("linear", "bilinear", "nearest_dtos"):
-            # 2) Linear and nearest neighbour methods:
+        elif method in ("linear", "bilinear"):
+            # 2) Linear methods:
             #
             # Mask out any row j that contains at least one positive
             # (i.e. greater than or equal to 'min_weight') w_ji that
@@ -545,7 +548,9 @@ def _regrid(
             # Note: It is much more efficient to access
             #       'weights.indptr', 'weights.indices', and
             #       'weights.data' directly, rather than iterating
-            #       over rows of 'weights' and using 'weights.getrow'.
+            #       over rows of 'weights' and using
+            #       'weights.getrow'. Also, 'np.count_nonzero' is much
+            #       faster than 'np.any' and 'np.all'.
             count_nonzero = np.count_nonzero
             where = np.where
             indptr = weights.indptr.tolist()
@@ -561,12 +566,45 @@ def _regrid(
 
             del indptr, pos_data
 
+        elif method == "nearest_dtos":
+            # 3) Nearest neighbour dtos method:
+            #
+            # Set to 0 any weight that corresponds to a masked source
+            # grid cell.
+            #
+            # Mask out any row j for which all source grid cells are
+            # masked.
+            dst_size = weights.shape[0]
+            if dst_mask is None:
+                dst_mask = np.zeros((dst_size,), dtype=bool)
+            else:
+                dst_mask = dst_mask.copy()
+
+            # Note: It is much more efficient to access
+            #       'weights.indptr', 'weights.indices', and
+            #       'weights.data' directly, rather than iterating
+            #       over rows of 'weights' and using
+            #       'weights.getrow'. Also, 'np.count_nonzero' is much
+            #       faster than 'np.any' and 'np.all'.
+            count_nonzero = np.count_nonzero
+            indptr = weights.indptr.tolist()
+            indices = weights.indices
+            for j, (i0, i1) in enumerate(zip(indptr[:-1], indptr[1:])):
+                mask = src_mask[indices[i0:i1]]
+                n_masked = count_nonzero(mask)
+                if n_masked == mask.size:
+                    dst_mask[j] = True
+                elif n_masked:
+                    weights.data[np.arange(i0, i1)[mask]] = 0
+
+            del indptr
+
         elif method in (
             "patch",
             "conservative_2nd",
             "nearest_stod",
         ):
-            # 3) Patch recovery and second-order conservative methods:
+            # 4) Patch recovery and second-order conservative methods:
             #
             # A reference source data mask has already been
             # incorporated into the weights matrix, and 'a' is assumed

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/docstring/docstring.py

@@ -161,7 +161,9 @@
                   mapped to the closest destination point. A
                   destination point can be mapped to multiple source
                   points. Some destination points may not be
-                  mapped. Useful for regridding of categorical data.
+                  mapped. Each regridded value is the sum of its
+                  contributing source elements. Useful for binning or
+                  for categorical data.
 
                 * `None`: This is the default and can only be used
                   when *dst* is a `RegridOperator`.""",