Create nodes at lines intersections

pyproject.toml

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

src/snkit/network.py

@@ -290,12 +290,64 @@ 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()
+            # excluding itself
+            # note that __eq__ is used on purpose instead of equals()
+            # this is stricter: for geometries constructed in the same way
+            # WARNING: this excludes sensical self-intersections, to be solved.
+            if edge != hit:
+                intersection = edge.geometry.intersection(hit)
+
+            # then extract the intersection points
+            hits_points = intersection_endpoints(intersection, hits_points)
+
+        # 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 +599,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)
@@ -590,12 +647,42 @@ def line_endpoints(line):
     return start, end
 
 
+def intersection_endpoints(geom, output=[]):
+    """Return the points from an intersection geometry
+
+    It extracts the starting and ending points of intersection
+    geometries recursively and appends them to `output`.
+    This doesn't handle polygons or collections of polygons.
+    """
+    geom_type = geom.geom_type
+    if geom.is_empty:
+        pass
+    elif geom_type == "Point":
+        output.append(geom)
+    elif geom_type == "LineString":
+        start = Point(geom.coords[0])
+        end = Point(geom.coords[-1])
+        output.append(start)
+        output.append(end)
+    # recursively for collections of geometries
+    # note that there is no shared inheritance relationship
+    elif (
+        geom_type == "MultiPoint"
+        or geom_type == "MultiLineString"
+        or geom_type == "GeometryCollection"
+    ):
+        for geom_ in geom.geoms:
+            output = intersection_endpoints(geom_, output)
+
+    return output
+
+
 def split_edge_at_points(edge, points, tolerance=1e-9):
     """Split edge at point/multipoint"""
     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
@@ -696,14 +783,20 @@ def to_networkx(network, directed=False, weight_col=None):
         G.add_nodes_from(network.nodes.id.to_list())
 
         # add nodal positions from geom
-        for node_id, x, y in zip(network.nodes.id, network.nodes.geometry.x, network.nodes.geometry.y):
+        for node_id, x, y in zip(
+            network.nodes.id, network.nodes.geometry.x, network.nodes.geometry.y
+        ):
             G.nodes[node_id]["pos"] = (x, y)
 
         # get edges from network data
         if weight_col is None:
             # default to geometry length
             edges_as_list = list(
-                zip(network.edges.from_id, network.edges.to_id, network.edges.geometry.length)
+                zip(
+                    network.edges.from_id,
+                    network.edges.to_id,
+                    network.edges.geometry.length,
+                )
             )
         else:
             edges_as_list = list(