Enable line and point selection on multi-topology datasets.

Prefix name on x, y, s coordinates to avoid collisions.

Add tests for intersect_line, intersect_linestring.

docs/changelog.rst

@@ -17,6 +17,13 @@ Added
   :meth:`xugrid.UgridDataArrayAccessor.intersect_line`, and
   :meth:`xugrid.UgridDataArrayAccessor.intersect_linestring` have been added to
   intersect line and linestrings and extract the associated face data.
+
+Changed
+~~~~~~~
+
+- Selection operations along a line, or at point locations, will now prefix the
+  name of the grid in the x and y coordinates. This avoids name collisions when
+  multiple topologies are present in a dataset.
 
 [0.6.4] 2023-08-22
 ------------------

examples/selection.py

@@ -95,7 +95,7 @@
 def show_point_selection(uda, da):
     _, ax = plt.subplots()
     uda.ugrid.plot(ax=ax, vmin=-20, vmax=90, cmap="terrain")
-    ax.scatter(da["x"], da["y"], color="red")
+    ax.scatter(da["mesh2d_x"], da["mesh2d_y"], color="red")
     ax.set_aspect(1.0)
 
 
@@ -164,7 +164,7 @@ def show_point_selection(uda, da):
 def show_line_selection(uda, da, line_x=None, line_y=None):
     _, (ax0, ax1) = plt.subplots(ncols=2, figsize=(10, 5))
     uda.ugrid.plot(ax=ax0, vmin=-20, vmax=90, cmap="terrain")
-    da.plot(ax=ax1, x="s")
+    da.plot(ax=ax1, x="mesh2d_s")
     if line_x is None:
         ax0.axhline(line_y, color="red")
     elif line_y is None:
@@ -207,6 +207,7 @@ def show_line_selection(uda, da, line_x=None, line_y=None):
 # This will work for any type of shapely line:
 
 ring = shapely.geometry.Point(155_000.0, 463_000).buffer(50_000.0).exterior
+da = uda.ugrid.intersect_linestring(ring)
 show_line_selection(uda, da, *shapely.get_coordinates(ring).T)
 
 # %%
@@ -233,3 +234,5 @@ def show_line_selection(uda, da, line_x=None, line_y=None):
 # In general, index selection should only be performed on the "core" dimension
 # of the UGRID topology. This is the edge dimension for 1D topologies, and the
 # face dimension for 2D topologies.
+
+# %%

tests/test_ugrid2d.py

@@ -608,9 +608,9 @@ def test_sel_points(self):
         expected = xr.DataArray(
             data=[0, 3],
             coords={
-                "index": (dim, [0, 1]),
-                "x": (dim, x),
-                "y": (dim, y),
+                f"{NAME}_index": (dim, [0, 1]),
+                f"{NAME}_x": (dim, x),
+                f"{NAME}_y": (dim, y),
             },
             dims=[dim],
         )
@@ -620,7 +620,7 @@ 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)
-        assert np.array_equal(actual["index"], [1, 3])
+        assert np.array_equal(actual[f"{NAME}_index"], [1, 3])
 
     def test_validate_indexer(self):
         with pytest.raises(ValueError, match="slice stop should be larger than"):
@@ -701,9 +701,9 @@ def check_output(actual):
             expected = xr.DataArray(
                 data=[0],
                 coords={
-                    "index": (dim, [0]),
-                    "x": (dim, [0.5]),
-                    "y": (dim, [0.5]),
+                    f"{NAME}_index": (dim, [0]),
+                    f"{NAME}_x": (dim, [0.5]),
+                    f"{NAME}_y": (dim, [0.5]),
                 },
                 dims=[dim],
             )
@@ -725,17 +725,19 @@ def check_output(actual):
             expected = xr.DataArray(
                 data=[0, 0, 1, 2, 2, 3],
                 coords={
-                    "x": (dim, [0.4, 0.8, 1.2, 0.4, 0.8, 1.2]),
-                    "y": (dim, [0.5, 0.5, 0.5, 1.1, 1.1, 1.1]),
+                    f"{NAME}_x": (dim, [0.4, 0.8, 1.2, 0.4, 0.8, 1.2]),
+                    f"{NAME}_y": (dim, [0.5, 0.5, 0.5, 1.1, 1.1, 1.1]),
                 },
                 dims=[dim],
             )
             # This fails for some reason:
             # assert expected.equals(actual)
             assert np.array_equal(expected.values, actual.values)
             assert expected.dims == actual.dims
-            assert np.allclose(expected["y"].values, actual["y"].values)
-            assert np.allclose(expected["x"].values, actual["x"].values)
+            x = f"{NAME}_x"
+            y = f"{NAME}_y"
+            assert np.allclose(expected[y].values, actual[y].values)
+            assert np.allclose(expected[x].values, actual[x].values)
 
         x = [0.4, 0.8, 1.2]
         y = [0.5, 1.1]
@@ -760,9 +762,9 @@ def test_sel__edges_from_slice(self):
         expected = xr.DataArray(
             data=[0, 1],
             coords={
-                "x": (dim, [0.5, 1.5]),
-                "y": (dim, [0.5, 0.5]),
-                "s": (dim, [0.5, 1.5]),
+                f"{NAME}_x": (dim, [0.5, 1.5]),
+                f"{NAME}_y": (dim, [0.5, 0.5]),
+                f"{NAME}_s": (dim, [0.5, 1.5]),
             },
             dims=[dim],
         )
@@ -773,14 +775,56 @@ def test_sel__edges_from_slice(self):
         expected = xr.DataArray(
             data=[0, 2],
             coords={
-                "x": (dim, [0.5, 0.5]),
-                "y": (dim, [0.5, 1.25]),
-                "s": (dim, [0.5, 1.25]),
+                f"{NAME}_x": (dim, [0.5, 0.5]),
+                f"{NAME}_y": (dim, [0.5, 1.25]),
+                f"{NAME}_s": (dim, [0.5, 1.25]),
             },
             dims=[dim],
         )
         assert expected.equals(actual)
 
+    def test_intersect_line_error(self):
+        with pytest.raises(ValueError, match="Start and end coordinate pairs"):
+            self.grid.intersect_line(
+                obj=None, start=(0.0, 0.0, 0.0), end=(1.0, 1.0, 1.0)
+            )
+
+    def test_intersect_line(self):
+        grid = self.grid
+        obj = xr.DataArray([0, 1, 2, 3], dims=[grid.face_dimension])
+
+        p0 = (0.0, 0.0)
+        p1 = (2.0, 2.0)
+        actual = grid.intersect_line(obj, start=p0, end=p1)
+        sqrt2 = np.sqrt(2.0)
+        assert isinstance(actual, xr.DataArray)
+        assert actual.dims == (grid.face_dimension,)
+        assert np.array_equal(actual.to_numpy(), [0, 3])
+        assert np.allclose(actual[f"{NAME}_x"], [0.5, 1.25])
+        assert np.allclose(actual[f"{NAME}_y"], [0.5, 1.25])
+        assert np.allclose(actual[f"{NAME}_s"], [0.5 * sqrt2, 1.25 * sqrt2])
+
+        actual = grid.intersect_line(obj, start=p1, end=p0)
+        assert np.array_equal(actual.to_numpy(), [3, 0])
+
+    def test_intersect_linestring(self):
+        grid = self.grid
+        obj = xr.DataArray([0, 1, 2, 3], dims=[grid.face_dimension])
+        linestring = shapely.geometry.LineString(
+            [
+                [0.5, 0.5],
+                [1.5, 0.5],
+                [1.5, 1.5],
+            ]
+        )
+        actual = grid.intersect_linestring(obj, linestring)
+        assert isinstance(actual, xr.DataArray)
+        assert actual.dims == (grid.face_dimension,)
+        assert np.array_equal(actual.to_numpy(), [0, 1, 1, 3])
+        assert np.allclose(actual[f"{NAME}_x"], [0.75, 1.25, 1.5, 1.5])
+        assert np.allclose(actual[f"{NAME}_y"], [0.5, 0.5, 0.75, 1.25])
+        assert np.allclose(actual[f"{NAME}_s"], [0.25, 0.75, 1.25, 1.75])
+
 
 def test_topology_subset():
     grid = grid2d()

tests/test_ugrid_dataset.py

@@ -244,6 +244,32 @@ def test_sel(self):
         assert actual.shape == (1,)
         assert actual.ugrid.grid.n_face == 1
 
+    def test_intersect_line(self):
+        p0 = (0.0, 0.0)
+        p1 = (2.0, 2.0)
+        actual = self.uda.ugrid.intersect_line(start=p0, end=p1)
+        sqrt2 = np.sqrt(2.0)
+        assert isinstance(actual, xr.DataArray)
+        assert actual.dims == ("mesh2d_nFaces",)
+        assert np.allclose(actual["mesh2d_x"], [0.5, 1.25])
+        assert np.allclose(actual["mesh2d_y"], [0.5, 1.25])
+        assert np.allclose(actual["mesh2d_s"], [0.5 * sqrt2, 1.25 * sqrt2])
+
+    def test_intersect_linestring(self):
+        linestring = shapely.geometry.LineString(
+            [
+                [0.5, 0.5],
+                [1.5, 0.5],
+                [1.5, 1.5],
+            ]
+        )
+        actual = self.uda.ugrid.intersect_linestring(linestring)
+        assert isinstance(actual, xr.DataArray)
+        assert actual.dims == ("mesh2d_nFaces",)
+        assert np.allclose(actual["mesh2d_x"], [0.75, 1.25, 1.5, 1.5])
+        assert np.allclose(actual["mesh2d_y"], [0.5, 0.5, 0.75, 1.25])
+        assert np.allclose(actual["mesh2d_s"], [0.25, 0.75, 1.25, 1.75])
+
     def test_rasterize(self):
         actual = self.uda.ugrid.rasterize(resolution=0.5)
         x = [0.25, 0.75, 1.25, 1.75]
@@ -573,6 +599,36 @@ def test_rasterize(self):
         assert np.allclose(actual["x"], x)
         assert np.allclose(actual["y"], y)
 
+    def test_intersect_line(self):
+        p0 = (0.0, 0.0)
+        p1 = (2.0, 2.0)
+        actual = self.uds.ugrid.intersect_line(start=p0, end=p1)
+        sqrt2 = np.sqrt(2.0)
+        assert isinstance(actual, xr.Dataset)
+        assert actual.dims == {"mesh2d_nFaces": 2}
+        assert np.allclose(actual["mesh2d_x"], [0.5, 1.25])
+        assert np.allclose(actual["mesh2d_y"], [0.5, 1.25])
+        assert np.allclose(actual["mesh2d_s"], [0.5 * sqrt2, 1.25 * sqrt2])
+        assert "a" in actual
+        assert "b" in actual
+
+    def test_intersect_linestring(self):
+        linestring = shapely.geometry.LineString(
+            [
+                [0.5, 0.5],
+                [1.5, 0.5],
+                [1.5, 1.5],
+            ]
+        )
+        actual = self.uds.ugrid.intersect_linestring(linestring)
+        assert isinstance(actual, xr.Dataset)
+        assert actual.dims == {"mesh2d_nFaces": 4}
+        assert np.allclose(actual["mesh2d_x"], [0.75, 1.25, 1.5, 1.5])
+        assert np.allclose(actual["mesh2d_y"], [0.5, 0.5, 0.75, 1.25])
+        assert np.allclose(actual["mesh2d_s"], [0.25, 0.75, 1.25, 1.75])
+        assert "a" in actual
+        assert "b" in actual
+
     def test_partitioning(self):
         partitions = self.uds.ugrid.partition(n_part=2)
         assert len(partitions) == 2

xugrid/core/accessorbase.py

@@ -32,6 +32,18 @@ def set_crs():
     def to_crs():
         """ """
 
+    @abc.abstractmethod
+    def sel_points():
+        """ """
+
+    @abc.abstractmethod
+    def intersect_line():
+        """ """
+
+    @abc.abstractmethod
+    def intersect_linestring():
+        """ """
+
     @abc.abstractproperty
     def bounds():
         """ """

xugrid/core/dataset_accessor.py

@@ -232,6 +232,7 @@ def sel_points(self, x, y):
         Returns
         -------
         points: Union[xr.DataArray, xr.Dataset]
+            The name of the topology is prefixed in the x, y coordinates.
         """
         result = self.obj
         for grid in self.grids:
@@ -281,6 +282,52 @@ def rasterize_like(self, other: Union[xr.DataArray, xr.Dataset]) -> xr.Dataset:
             datasets.append(self._raster(xx, yy, index))
         return xr.merge(datasets)
 
+    def intersect_line(
+        self, start: Sequence[float], end: Sequence[float]
+    ) -> xr.Dataset:
+        """
+        Intersect a line with the grid of this data, and fetch the values of
+        the intersected faces.
+
+        Parameters
+        ----------
+        obj: xr.DataArray or xr.Dataset
+        start: sequence of two floats
+            coordinate pair (x, y), designating the start point of the line.
+        end: sequence of two floats
+            coordinate pair (x, y), designating the end point of the line.
+
+        Returns
+        -------
+        intersection: xr.Dataset
+            The name of the topology is prefixed in the x, y and s
+            (spatium=distance) coordinates.
+        """
+        obj = self.obj
+        for grid in self.grids:
+            obj = grid.intersect_line(obj, start, end)
+        return obj
+
+    def intersect_linestring(self, linestring) -> xr.Dataset:
+        """
+        Intersect the grid along a collection of linestrings. Returns a new Dataset
+        with the values for each intersected segment.
+
+        Parameters
+        ----------
+        linestring: shapely.LineString
+
+        Returns
+        -------
+        intersection: xr.Dataset
+            The name of the topology is prefixed in the x, y and s
+            (spatium=distance) coordinates.
+        """
+        obj = self.obj
+        for grid in self.grids:
+            obj = grid.intersect_linestring(obj, linestring)
+        return obj
+
     def to_dataset(self, optional_attributes: bool = False):
         """
         Converts this UgridDataArray or UgridDataset into a standard

xugrid/ugrid/ugrid1d.py

@@ -356,9 +356,6 @@ def to_shapely(self, dim):
                 " Ugrid1d topology."
             )
 
-    def sel_points(obj, x, y):
-        raise NotImplementedError("Cannot select points in a UGRID 1D topology")
-
     def isel(self, indexers=None, return_index=False, **indexers_kwargs):
         """
         Select based on node or edge.
@@ -523,6 +520,15 @@ def clip_box(
     ):
         return self.sel(x=slice(xmin, xmax), y=slice(ymin, ymax))
 
+    def sel_points(obj, x, y):
+        return obj
+
+    def intersect_line(self, obj, start, stop):
+        return obj
+
+    def intersect_linestring(self, obj, linestring):
+        return obj
+
     def topological_sort_by_dfs(self) -> IntArray:
         """
         Returns an array of vertices in topological order.