Merge pull request #271 from Deltares/sel-points-bounds

Do not drop out-of-bounds point silently in sel_points

docs/changelog.rst

@@ -30,6 +30,14 @@ Changed
 
 - Custom reduction functions provide to the overlap regridders no longer require
   an ``indices`` argument.
+- :meth:`xugrid.Ugrid2d.sel_points`,
+  :meth:`xugrid.UgridDataArrayAccessor.sel_points` and
+  :meth:`xugrid.UgridDatasetAccessor.sel_points` now take an ``out_of_bounds``
+  and ``fill_value`` argument to determine what to with points that do not fall
+  inside of any grid feature. Previously, the method silently dropped these
+  points. The method now takes a ``fill_value`` argument to assign to
+  out-of-bounds points. It gives a warning return uses ``fill_value=np.nan`` by
+  default. To enable the old behavior, set ``out_of_bounds="drop"``.
 
 [0.11.0] 2024-08-05
 -------------------

tests/test_ugrid1d.py

@@ -342,6 +342,20 @@ def test_sel():
     assert np.array_equal(new_obj.values, [0])
 
 
+def test_sel_points():
+    grid = grid1d()
+    obj = xr.DataArray(
+        data=[0, 1],
+        dims=[grid.edge_dimension],
+    )
+    # For now, this function does nothing so it'll work for multi-topology
+    # UgridDatasets.
+    actual = grid.sel_points(
+        obj=obj, x=None, y=None, out_of_bounds=None, fill_value=None
+    )
+    assert actual.identical(obj)
+
+
 def test_topology_subset():
     grid = grid1d()
     edge_indices = np.array([1])

tests/test_ugrid2d.py

@@ -670,6 +670,8 @@ def test_sel_points(self):
         x = [0.5, 1.5]
         y = [0.5, 1.25]
 
+        with pytest.raises(ValueError, match="out_of_bounds must be one of"):
+            self.grid.sel_points(obj=self.obj, x=x, y=y, out_of_bounds="nothing")
         with pytest.raises(ValueError, match="shape of x does not match shape of y"):
             self.grid.sel_points(obj=self.obj, x=[0.5, 1.5], y=[0.5])
         with pytest.raises(ValueError, match="x and y must be 1d"):
@@ -693,9 +695,29 @@ def test_sel_points(self):
     def test_sel_points_out_of_bounds(self):
         x = [-10.0, 0.5, -20.0, 1.5, -30.0]
         y = [-10.0, 0.5, -20.0, 1.25, -30.0]
-        actual = self.grid.sel_points(obj=self.obj, x=x, y=y)
+
+        with pytest.raises(
+            ValueError, match="Not all points are located inside of the grid"
+        ):
+            self.grid.sel_points(obj=self.obj, x=x, y=y, out_of_bounds="raise")
+
+        actual = self.grid.sel_points(obj=self.obj, x=x, y=y, out_of_bounds="drop")
         assert np.array_equal(actual[f"{NAME}_index"], [1, 3])
 
+        with pytest.warns(
+            UserWarning, match="Not all points are located inside of the grid"
+        ):
+            actual = self.grid.sel_points(obj=self.obj, x=x, y=y, out_of_bounds="warn")
+            assert np.allclose(actual, [np.nan, 0, np.nan, 3, np.nan], equal_nan=True)
+
+        actual = self.grid.sel_points(obj=self.obj, x=x, y=y, out_of_bounds="ignore")
+        assert np.allclose(actual, [np.nan, 0, np.nan, 3, np.nan], equal_nan=True)
+
+        actual = self.grid.sel_points(
+            obj=self.obj, x=x, y=y, out_of_bounds="ignore", fill_value=-1
+        )
+        assert np.allclose(actual, [-1, 0, -1, 3, -1])
+
     def test_validate_indexer(self):
         with pytest.raises(ValueError, match="slice stop should be larger than"):
             self.grid._validate_indexer(slice(2, 0))

xugrid/core/accessorbase.py

@@ -38,7 +38,7 @@ def sel(self):
         pass
 
     @abc.abstractmethod
-    def sel_points(self):
+    def sel_points(self, x, y, out_of_bounds, fill_value):
         pass
 
     @abc.abstractmethod

xugrid/core/dataarray_accessor.py

@@ -1,5 +1,6 @@
 from typing import Dict, List, Sequence, Tuple, Union
 
+import numpy as np
 import scipy.sparse
 import xarray as xr
 
@@ -175,7 +176,7 @@ def sel(self, x=None, y=None):
         else:
             return result
 
-    def sel_points(self, x, y):
+    def sel_points(self, x, y, out_of_bounds="warn", fill_value=np.nan):
         """
         Select points in the unstructured grid.
 
@@ -186,12 +187,23 @@ def sel_points(self, x, y):
         ----------
         x: ndarray of floats with shape ``(n_points,)``
         y: ndarray of floats with shape ``(n_points,)``
+        out_of_bounds: str, default: "warn"
+            What to do when points are located outside of any feature:
+
+            * raise: raise a ValueError.
+            * ignore: return ``fill_value`` for the out of bounds points.
+            * warn: give a warning and return NaN for the out of bounds points.
+            * drop: drop the out of bounds points. They may be identified
+              via the ``index`` coordinate of the returned selection.
+        fill_value: scalar, DataArray, Dataset, or callable, optional, default: np.nan
+            Value to assign to out-of-bounds points if out_of_bounds is warn
+            or ignore. Forwarded to xarray's ``.where()`` method.
 
         Returns
         -------
         points: Union[xr.DataArray, xr.Dataset]
         """
-        return self.grid.sel_points(self.obj, x, y)
+        return self.grid.sel_points(self.obj, x, y, out_of_bounds, fill_value)
 
     def rasterize(self, resolution: float) -> xr.DataArray:
         """

xugrid/core/dataset_accessor.py

@@ -229,7 +229,7 @@ def sel(self, x=None, y=None):
         else:
             return result
 
-    def sel_points(self, x, y):
+    def sel_points(self, x, y, out_of_bounds="warn", fill_value=np.nan):
         """
         Select points in the unstructured grid.
 
@@ -240,6 +240,17 @@ def sel_points(self, x, y):
         ----------
         x: ndarray of floats with shape ``(n_points,)``
         y: ndarray of floats with shape ``(n_points,)``
+        out_of_bounds: str, default ``"warn"``
+            What to do when points are located outside of any feature:
+
+            * raise: raise a ValueError.
+            * ignore: return ``fill_value`` for the out of bounds points.
+            * warn: give a warning and return NaN for the out of bounds points.
+            * drop: drop the out of bounds points. They may be identified
+              via the ``index`` coordinate of the returned selection.
+        fill_value: scalar, DataArray, Dataset, or callable, optional, default: np.nan
+            Value to assign to out-of-bounds points if out_of_bounds is warn
+            or ignore. Forwarded to xarray's ``.where()`` method.
 
         Returns
         -------
@@ -248,7 +259,7 @@ def sel_points(self, x, y):
         """
         result = self.obj
         for grid in self.grids:
-            result = grid.sel_points(result, x, y)
+            result = grid.sel_points(result, x, y, out_of_bounds, fill_value)
         return result
 
     def rasterize(self, resolution: float) -> xr.Dataset:

xugrid/ugrid/ugrid1d.py

@@ -574,7 +574,7 @@ def clip_box(
     ):
         return self.sel(x=slice(xmin, xmax), y=slice(ymin, ymax))
 
-    def sel_points(obj, x, y):
+    def sel_points(self, obj, x, y, out_of_bounds, fill_value):
         return obj
 
     def intersect_line(self, obj, start, stop):

xugrid/ugrid/ugrid2d.py

@@ -1,5 +1,6 @@
 from __future__ import annotations
 
+import warnings
 from itertools import chain
 from typing import Any, Dict, Optional, Sequence, Tuple, Union
 
@@ -1292,24 +1293,43 @@ def _sel_xline(
         ystop = numeric_bound(y.stop, ymax)
         return self._sel_line(obj, start=(x, ystart), end=(x, ystop))
 
-    def sel_points(self, obj, x: FloatArray, y: FloatArray):
+    def sel_points(
+        self, obj, x: FloatArray, y: FloatArray, out_of_bounds="warn", fill_value=np.nan
+    ):
         """
         Select points in the unstructured grid.
 
-        Out-of-bounds points are ignored. They may be identified via the
-        ``index`` coordinate of the returned selection.
 
         Parameters
         ----------
         x: 1d array of floats with shape ``(n_points,)``
         y: 1d array of floats with shape ``(n_points,)``
         obj: xr.DataArray or xr.Dataset
+        out_of_bounds: str, default ``"warn"``
+            What to do when points are located outside of any feature:
+
+            * raise: raise a ValueError.
+            * ignore: return ``fill_value`` for the out of bounds points.
+            * warn: give a warning and return NaN for the out of bounds points.
+            * drop: drop the out of bounds points. They may be identified
+              via the ``index`` coordinate of the returned selection.
+        fill_value: scalar, DataArray, Dataset, or callable, optional, default: np.nan
+            Value to assign to out-of-bounds points if out_of_bounds is warn
+            or ignore. Forwarded to xarray's ``.where()`` method.
 
         Returns
         -------
         selection: xr.DataArray or xr.Dataset
             The name of the topology is prefixed in the x, y coordinates.
         """
+        options = ("warn", "raise", "ignore", "drop")
+        if out_of_bounds not in options:
+            str_options = ", ".join(options)
+            raise ValueError(
+                f"out_of_bounds must be one of {str_options}, "
+                f"received: {out_of_bounds}"
+            )
+
         x = np.atleast_1d(x)
         y = np.atleast_1d(y)
         if x.shape != y.shape:
@@ -1319,14 +1339,34 @@ def sel_points(self, obj, x: FloatArray, y: FloatArray):
         dim = self.face_dimension
         xy = np.column_stack([x, y])
         index = self.locate_points(xy)
+
+        keep = slice(None, None)  # keep all by default
+        condition = None
         valid = index != -1
-        index = index[valid]
+        if not valid.all():
+            msg = "Not all points are located inside of the grid."
+            if out_of_bounds == "raise":
+                raise ValueError(msg)
+            elif out_of_bounds in ("warn", "ignore"):
+                if out_of_bounds == "warn":
+                    warnings.warn(msg)
+                condition = xr.DataArray(valid, dims=(dim,))
+            elif out_of_bounds == "drop":
+                index = index[valid]
+                keep = valid
+
+        # Create the selection DataArray or Dataset
         coords = {
-            f"{self.name}_index": (dim, np.arange(len(valid))[valid]),
-            f"{self.name}_x": (dim, xy[valid, 0]),
-            f"{self.name}_y": (dim, xy[valid, 1]),
+            f"{self.name}_index": (dim, np.arange(len(xy))[keep]),
+            f"{self.name}_x": (dim, xy[keep, 0]),
+            f"{self.name}_y": (dim, xy[keep, 1]),
         }
-        return obj.isel({dim: index}).assign_coords(coords)
+        selection = obj.isel({dim: index}).assign_coords(coords)
+
+        # Set values to fill_value for out-of-bounds
+        if condition is not None:
+            selection = selection.where(condition, other=fill_value)
+        return selection
 
     def intersect_line(self, obj, start: Sequence[float], end: Sequence[float]):
         """

xugrid/ugrid/ugridbase.py

@@ -129,7 +129,7 @@ def clip_box(self):
         pass
 
     @abc.abstractmethod
-    def sel_points(self):
+    def sel_points(self, obj, x, y, out_bounds, fill_value):
         pass
 
     @abc.abstractmethod