AbstractUgrid.intersect_linestring cannot intersect with multiple disconnected linestrings

[JoerivanEngelen]

I was a bit confused with the AbstractUgrid.intersect_linestring method. The docstring mentions it is to "Intersect linestrings" (plural). Looking at the code and the method name, it is to intersect with a single linestring.

I was trying to intersect multiple cross-sections linestrings with a network, and ran into the problem that it returned too much intersections. This is caused by:

xugrid/xugrid/ugrid/ugridbase.py

Line 1222 in 3286ff3

edges = np.stack((xy[:-1], xy[1:]), axis=1)

which connects all seperate linestring with each other.

This seems to work to create edges where linestrings are not connected:

    def intersect_linestring(
        self,
        obj: Union[xr.DataArray, xr.Dataset],
        linestring: "shapely.geometry.LineString",  # type: ignore # noqa
    ) -> Union[xr.DataArray, xr.Dataset]:
...
        xys, linesection_ids = shapely.get_coordinates([linestrings], return_index=True)
        _, split_ids = np.unique(linesection_ids, return_index=True)
        split_xys = np.split(xys, split_ids[1:])
        
        edges = np.concatenate([np.stack((xy[:-1], xy[1:]), axis=1) for xy in split_xys], axis=0)
        edge_line_ids = np.delete(linesection_ids, split_ids)
        edge_index, core_index, intersections = self.intersect_edges(edges)
        # Linestring ids, which can be added to coords
        line_ids = edge_line_ids[edge_index]
        ...
        coords, core_index = get_sorted_section_coords(
            s, intersection_for_coord, line_ids, dim, core_index, self.name
        )
        ...

I also thought it is useful to keep the linestring ids and add them to coords:

def get_sorted_section_coords(
    s: FloatArray, xy: FloatArray, line_ids: IntArray, dim: str, index: IntArray, name: str
):
    order = np.argsort(s)
    coords = {
        f"{name}_x": (dim, xy[order, 0]),
        f"{name}_y": (dim, xy[order, 1]),
        f"{name}_s": (dim, s[order]),
        f"{name}_line_id": (dim, line_ids[order]),
    }
    return coords, index[order]

@Huite Let me know if this is within the scope of intersect_linestring or requires a separate method intersect_linestrings.

[Huite]

Uh, oh right. I guess I had a single linestring in mind, because of the signature, the stack logic, and the user guide part: https://deltares.github.io/xugrid/examples/selection.html#line-selection

I'm not entirely sure whether we should create a new method intersect_linestrings (multiple) or whether we should check the type; seems easy to call np.atleast_1d and then assume it's multiple linestrings all the time?

This function is probably what we should use for the linestring to edges conversion in case of multiple linestrings:

xugrid/xugrid/conversion.py

Line 70 in 3286ff3

def linestrings_to_edges(edges: LineArray) -> Tuple[FloatArray, FloatArray, IntArray]:

It should be more efficient than numpy split since everything's vectorized. The version there returns x and y coordinates via a unique call though, I think you only want this logic:

    xy, index = shapely.get_coordinates(edges, return_index=True)
    linear_index = np.arange(index.size)
    segments = np.column_stack([linear_index[:-1], linear_index[1:]])
    segments = segments[np.diff(index) == 0]

You could probably pull that into a separate function?

[JoerivanEngelen]

Thanks! I was struggling to find a vectorized solution quickly, so this certainly helps. The vectorized version of the required logic would then be:

    def intersect_linestring(
        self,
        obj: Union[xr.DataArray, xr.Dataset],
        linestring: "shapely.geometry.LineString",  # type: ignore # noqa
    ) -> Union[xr.DataArray, xr.Dataset]:
...
        xys, linesection_ids = shapely.get_coordinates([linestrings], return_index=True)
        linear_index = np.arange(linesection_ids.size)
        segments = np.column_stack([linear_index[:-1], linear_index[1:]])
        is_same_section = np.diff(linesection_ids) == 0
        segments = segments[is_same_section]
        edges = xys[segments]
        edge_line_ids = linesection_ids[1:][is_same_section]
        edge_index, core_index, intersections = self.intersect_edges(edges)
        # Linestring ids, which can be added to coords
        line_ids = edge_line_ids[edge_index]
        ...
        coords, core_index = get_sorted_section_coords(
            s, intersection_for_coord, line_ids, dim, core_index, self.name
        )
        ...

It is better to refactor this code more to reduce duplication with linestrings_to_edges.