dev

Author: davidhassell

Changelog.rst

@@ -3,6 +3,9 @@ version NEXTVERSION
 
 **2024-??-??**
 
+* Allowed ``'nearest_dtos'`` 2-d regridding to work with discrete
+  sampling geometry source grids ()
+  (https://github.com/NCAS-CMS/cf-python/issues/811)
 * New method: `cf.Field.filled`
   (https://github.com/NCAS-CMS/cf-python/issues/811)
 * New method: `cf.Field.is_discrete_axis`

cf/cfimplementation.py

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

cf/data/dask_regrid.py

@@ -514,10 +514,11 @@ def _regrid(
             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
 
@@ -529,8 +530,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
@@ -562,12 +563,43 @@ 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'.
+            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/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/docstring/docstring.py

@@ -161,7 +161,11 @@
                   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 the
+                  contributing source elements. Masked source grid
+                  cells are mapped, but do not contribute to the
+                  regridded values. Useful for binning or for
+                  categorical data.
 
                 * `None`: This is the default and can only be used
                   when *dst* is a `RegridOperator`.""",

cf/regrid/regrid.py

@@ -523,11 +523,10 @@ def regrid(
                     "are a UGRID mesh"
                 )
 
-            if src_grid.is_locstream or dst_grid.is_locstream:
+            if dst_grid.is_locstream:
                 raise ValueError(
-                    f"{method!r} regridding is (at the moment) only available "
-                    "when neither the source and destination grids are "
-                    "DSG featureTypes."
+                    f"{method!r} regridding is (at the moment) not available "
+                    "when the destination grid is a DSG featureType."
                 )
 
     elif cartesian and (src_grid.is_mesh or dst_grid.is_mesh):
@@ -656,6 +655,7 @@ def regrid(
             dst=dst,
             weights_file=weights_file if from_file else None,
             src_mesh_location=src_grid.mesh_location,
+            src_featureType=src_grid.featureType,
             dst_featureType=dst_grid.featureType,
             src_z=src_grid.z,
             dst_z=dst_grid.z,
@@ -1279,7 +1279,7 @@ def spherical_grid(
 
     # Set cyclicity of X axis
     if mesh_location or featureType:
-        cyclic = None
+        cyclic = False
     elif cyclic is None:
         cyclic = f.iscyclic(x_axis)
     else:
@@ -2281,6 +2281,11 @@ def create_esmpy_locstream(grid, mask=None):
         #       but the esmpy mask requires 0/1 for masked/unmasked
         #       elements.
         mask = np.invert(mask).astype("int32")
+        if mask.size == 1:
+            # Make sure that there's a mask element for each point in
+            # the locstream (rather than a scalar that applies to all
+            # elements).
+            mask = np.full((location_count,), mask, dtype="int32")
     else:
         # No masked points
         mask = np.full((location_count,), 1, dtype="int32")
@@ -2465,7 +2470,6 @@ def create_esmpy_weights(
             from netCDF4 import Dataset
 
             from .. import __version__
-
             from ..data.array.locks import netcdf_lock
 
             if (

cf/test/test_regrid_featureType.py

@@ -28,7 +28,7 @@
 disallowed_methods = (
     "conservative",
     "conservative_2nd",
-    "nearest_dtos",
+    #    "nearest_dtos",
 )
 
 methods = (
@@ -161,6 +161,75 @@ def test_Field_regrid_grid_to_featureType_3d(self):
 
                     y = esmpy_regrid(coord_sys, method, src, dst, **kwargs)
 
+                    self.assertEqual(y.size, a.size)
+                    self.assertTrue(np.allclose(y, a, atol=4e-3, rtol=rtol))
+
+                    if isinstance(a, np.ma.MaskedArray):
+                        self.assertTrue((y.mask == a.mask).all())
+                    else:
+                        self.assertFalse(y.mask.any())
+
+    @unittest.skipUnless(esmpy_imported, "Requires esmpy/ESMF package.")
+    def test_Field_regrid_featureType_to_grid_2d(self):
+        self.assertFalse(cf.regrid_logging())
+
+        src = self.dst_featureType
+        src.del_construct("cellmethod0")
+        src = src[:12]
+        src[...] = 273 + np.arange(12)
+        x = src.coord("X")
+        x[...] = [4, 6, 9, 11, 14, 16, 4, 6, 9, 11, 14, 16]
+        y = src.coord("Y")
+        y[...] = [41, 41, 31, 31, 21, 21, 39, 39, 29, 29, 19, 19]
+
+        dst = self.src_grid.copy()
+        x = dst.coord("X")
+        x[...] = [5, 10, 15, 20]
+        y = dst.coord("Y")
+        y[...] = [10, 20, 30, 40]
+
+        # Mask some destination grid points
+        dst[0, 0, 1, 2] = cf.masked
+
+        y0 = np.ma.array(
+            [[0, 0, 0, 0], [0, 0, 1122, 0], [0, 1114, 0, 0], [1106, 0, 0, 0]],
+            mask=[
+                [True, True, True, True],
+                [True, True, False, True],
+                [True, False, True, True],
+                [False, True, True, True],
+            ],
+        )
+
+        coord_sys = "spherical"
+
+        for src_masked in (False, True):
+            y = y0.copy()
+            if src_masked:
+                src = src.copy()
+                src[6:8] = cf.masked
+                y[3, 0] = 547
+
+            # Loop over whether or not to use the destination grid
+            # masked points
+            for use_dst_mask in (False, True):
+                if use_dst_mask:
+                    y = y.copy()
+                    y[1, 2] = np.ma.masked
+
+                kwargs = {"use_dst_mask": use_dst_mask}
+                method = "nearest_dtos"
+                for return_operator in (False, True):
+                    if return_operator:
+                        r = src.regrids(
+                            dst, method=method, return_operator=True, **kwargs
+                        )
+                        x = src.regrids(r)
+                    else:
+                        x = src.regrids(dst, method=method, **kwargs)
+
+                    a = x.array
+
                     self.assertEqual(y.size, a.size)
                     self.assertTrue(np.allclose(y, a, atol=atol, rtol=rtol))
 
@@ -196,7 +265,6 @@ def test_Field_regrid_grid_to_featureType_2d(self):
                     a = x.array
 
                     y = esmpy_regrid(coord_sys, method, src, dst, **kwargs)
-
                     self.assertEqual(y.size, a.size)
                     self.assertTrue(np.allclose(y, a, atol=atol, rtol=rtol))