Create nodes at lines intersections

pyproject.toml

@@ -0,0 +1,4 @@
+[tool.pytest.ini_options]
+pythonpath = [
+    "src",
+]

src/snkit/network.py

@@ -290,12 +290,77 @@ def split_edges_at_nodes(network, tolerance=1e-9):
 
     # combine dfs
     edges = pandas.concat(split_edges, axis=0)
-    # reset index and drop
     edges = edges.reset_index().drop("index", axis=1)
-    # return new network with split edges
+
     return Network(nodes=network.nodes, edges=edges)
 
 
+def split_edges_at_intersections(network, tolerance=1e-9):
+    """Split network edges where they intersect line geometries"""
+    split_edges = []
+    split_points = []
+    for edge in tqdm(
+        network.edges.itertuples(index=False), desc='split', total=len(network.edges)
+    ):
+        # note: the symmetry of intersection is not exploited here.
+        # (If A intersects B, then B intersects A)
+        # since edges are not modified within the loop, this has just
+        # potential performance consequences.
+
+        hits = edges_intersecting(edge.geometry, network.edges, tolerance)
+
+        hits_points = []
+        for hit in hits.geometry:
+            # first extract the actual intersections from the hits
+            # these are *new geometrical objects* not in the sindex
+            # therefore they cannot be returned by internal _intersects*
+            intersection = Point()
+            # restrict intersection to crossing
+            # this excludes self overlap and end-points intersections
+            # WARNING: how are self-crossing handled? (e.g. loop-like)
+            if edge.geometry.crosses(hit):
+                intersection = edge.geometry.intersection(hit)
+
+            # then keep track of the intersection points
+            geom_type = intersection.geom_type
+            if intersection.is_empty:
+                continue
+            elif geom_type == 'Point':
+                hits_points.append(intersection)
+            elif geom_type == 'MultiPoint':
+                for point in intersection.geoms:
+                    hits_points.append(point)
+            elif geom_type == 'MultiLineString':
+                # when lines almost overlap for a stretch
+                for line in intersection.geoms:
+                    start = Point(line.coords[0])
+                    end = Point(line.coords[-1])
+                    hits_points.append(start)
+                    hits_points.append(end)
+
+        # store the split edges and intersection points
+        split_points.extend(hits_points)
+        hits_points = MultiPoint(hits_points)
+        edges = split_edge_at_points(edge, hits_points, tolerance)
+        split_edges.append(edges)
+
+    # add the (potentially) split edges
+    edges = pandas.concat(split_edges, axis=0)
+    edges = edges.reset_index().drop("index", axis=1)
+
+    # combine the original nodes with the new intersection nodes
+    # dropping the duplicates.
+    # note: there are at least duplicates from above since intersections
+    # are checked twice
+    # note: intersection nodes are appended, and if any duplicates, the
+    # original counterparts are kept.
+    nodes = GeoDataFrame(geometry=split_points)
+    nodes = pandas.concat([network.nodes, nodes], axis=0).drop_duplicates()
+    nodes = nodes.reset_index().drop("index", axis=1)
+
+    return Network(nodes=nodes, edges=edges)
+
+
 def link_nodes_to_edges_within(network, distance, condition=None, tolerance=1e-9):
     """Link nodes to all edges within some distance"""
     new_node_geoms = []
@@ -547,6 +612,11 @@ def nodes_intersecting(line, nodes, tolerance=1e-9):
     return intersects(line, nodes, tolerance)
 
 
+def edges_intersecting(line, edges, tolerance=1e-9):
+    """Find edges intersecting line"""
+    return intersects(line, edges, tolerance)
+
+
 def intersects(geom, gdf, tolerance=1e-9):
     """Find the subset of a GeoDataFrame intersecting with a shapely geometry"""
     return _intersects(geom, gdf, tolerance)
@@ -595,7 +665,7 @@ def split_edge_at_points(edge, points, tolerance=1e-9):
     try:
         segments = split_line(edge.geometry, points, tolerance)
     except ValueError:
-        # if splitting fails, e.g. becuase points is empty GeometryCollection
+        # if splitting fails, e.g. because points is empty GeometryCollection
         segments = [edge.geometry]
     edges = GeoDataFrame([edge] * len(segments))
     edges.geometry = segments

tests/test_init.py

@@ -140,6 +140,49 @@ def split_with_ids():
     return snkit.Network(edges=edges, nodes=nodes)
 
 
+@fixture
+def unsplit_intersection():
+    """Edges intersection, both edges not split
+       b
+       |
+    c--|--d
+       |
+       a
+    """
+    a = Point((1, 0))
+    b = Point((1, 2))
+    c = Point((0, 1))
+    d = Point((2, 1))
+    nodes = GeoDataFrame(data={"geometry": [a, b, c, d]})
+    ab = LineString([a, b])
+    cd = LineString([c, d])
+    edges = GeoDataFrame(data={"geometry": [ab, cd]})
+    return snkit.Network(edges=edges, nodes=nodes)
+
+
+@fixture
+def split_intersection():
+    """Edges intersection, both edges split
+       b
+       |
+    c--x--d
+       |
+       a
+    """
+    a = Point((1, 0))
+    b = Point((1, 2))
+    c = Point((0, 1))
+    d = Point((2, 1))
+    x = Point((1, 1))
+    nodes = GeoDataFrame(data={"geometry": [a, b, c, d, x]})
+    ax = LineString([a, x])
+    xb = LineString([x, b])
+    cx = LineString([c, x])
+    xd = LineString([x, d])
+    edges = GeoDataFrame(data={"geometry": [ax, xb, cx, xd]})
+    return snkit.Network(edges=edges, nodes=nodes)
+
+
 @fixture
 def gap():
     """T-junction with nodes, edges not quite intersecting:
@@ -229,6 +272,13 @@ def test_split_at_nodes(unsplit, split):
     assert_frame_equal(split.edges, actual.edges)
 
 
+def test_split_at_intersections(unsplit_intersection, split_intersection):
+    """Should split edges at edges intersections"""
+    actual = snkit.network.split_edges_at_intersections(unsplit_intersection)
+    assert_frame_equal(split_intersection.edges, actual.edges)
+    assert_frame_equal(split_intersection.nodes, actual.nodes)
+
+
 def test_split_line():
     """Definitively split line at a point"""
     line = LineString(