Add inject_points for projecting external points onto the network

docs/source/api.rst

@@ -41,6 +41,7 @@ A subset of functions dealing with network nodes:
    consolidate_nodes
    remove_interstitial_nodes
    induce_nodes
+   inject_points
 
 Face artifact detection
 -----------------------

neatnet/__init__.py

@@ -8,6 +8,7 @@
     consolidate_nodes,
     fix_topology,
     induce_nodes,
+    inject_points,
     remove_interstitial_nodes,
     split,
 )

neatnet/nodes.py

@@ -323,6 +323,10 @@ def induce_nodes(streets: gpd.GeoDataFrame, *, eps: float = 1e-4) -> gpd.GeoData
     -------
     geopandas.GeoDataFrame
         Updated ``streets`` with (potentially) added nodes.
+
+    See Also
+    --------
+    inject_points : add nodes from external Point features (POIs).
     """
 
     sindex_kws = {"predicate": "dwithin", "distance": 1e-4}
@@ -343,6 +347,78 @@ def induce_nodes(streets: gpd.GeoDataFrame, *, eps: float = 1e-4) -> gpd.GeoData
     return split(nodes_to_induce.geometry, streets, streets.crs, eps=eps)
 
 
+def inject_points(
+    streets: gpd.GeoDataFrame,
+    points: gpd.GeoDataFrame,
+    *,
+    snap_radius: float | None = None,
+    eps: float = 1e-4,
+) -> gpd.GeoDataFrame:
+    """Project external points onto the nearest LineString and split there.
+
+    For every input point within ``snap_radius`` of any LineString in
+    ``streets``, the point is projected onto the nearest LineString and the
+    target line is split at that projected location. Points outside the
+    radius are ignored. When the CRS differs, ``points`` is reprojected to
+    match ``streets`` before processing.
+
+    Parameters
+    ----------
+    streets : geopandas.GeoDataFrame
+        LineString network. CRS must be set.
+    points : geopandas.GeoDataFrame
+        Point features to project onto ``streets``. CRS must be set; will
+        be reprojected to match ``streets`` if different.
+    snap_radius : float | None = None
+        Maximum projection distance, in ``streets`` CRS units. Points
+        beyond this distance from any LineString are dropped from the
+        result. ``None`` keeps every point regardless of distance.
+    eps : float = 1e-4
+        Tolerance epsilon passed to :func:`split` for the actual snap.
+
+    Returns
+    -------
+    geopandas.GeoDataFrame
+        ``streets`` with the relevant LineStrings split at projection
+        locations.
+
+    See Also
+    --------
+    induce_nodes : symmetric case where the new node comes from a LineString endpoint.
+    split : underlying split-at-points primitive.
+    """
+    if streets.crs is None or points.crs is None:
+        raise ValueError("Both streets and points must have a CRS set.")
+    if points.crs != streets.crs:
+        points = points.to_crs(streets.crs)
+
+    line_geoms = streets.geometry.values
+    point_geoms = points.geometry.values
+    if len(point_geoms) == 0 or len(line_geoms) == 0:
+        return streets.copy()
+
+    # Nearest line per input point (sindex-accelerated)
+    input_idx, line_idx = streets.sindex.nearest(point_geoms, return_all=False)
+    distances = shapely.distance(line_geoms[line_idx], point_geoms[input_idx])
+
+    if snap_radius is not None:
+        keep = distances <= snap_radius
+        input_idx = input_idx[keep]
+        line_idx = line_idx[keep]
+
+    if len(input_idx) == 0:
+        return streets.copy()
+
+    # Project each accepted point onto its nearest LineString (foot of perpendicular)
+    accepted_lines = line_geoms[line_idx]
+    accepted_points = point_geoms[input_idx]
+    proj_dists = shapely.line_locate_point(accepted_lines, accepted_points)
+    projected = shapely.line_interpolate_point(accepted_lines, proj_dists)
+
+    projected_series = gpd.GeoSeries(projected, crs=streets.crs)
+    return split(projected_series, streets, streets.crs, eps=eps)
+
+
 def _identify_degree_mismatch(
     edges: gpd.GeoDataFrame, sindex_kws: dict
 ) -> gpd.GeoSeries:

neatnet/tests/test_nodes.py

@@ -1316,3 +1316,109 @@ def test_fill_attrs():
     )
 
     assert_frame_equal(known, observed)
+
+
+# ----- inject_points -----
+
+
+_inj_line_h = shapely.LineString([(0, 0), (100, 0)])
+_inj_line_v = shapely.LineString([(50, 0), (50, 100)])
+
+
+def _inj_network():
+    return geopandas.GeoDataFrame(
+        {"road": ["H", "V"]},
+        geometry=[_inj_line_h, _inj_line_v],
+        crs="EPSG:3035",
+    )
+
+
+@pytest.mark.parametrize(
+    ("point", "snap_radius", "expected_h_segments"),
+    [
+        (shapely.Point(40, 0), 1, 2),
+        (shapely.Point(40, 1), 5, 2),
+    ],
+)
+def test_inject_points_projects_to_target_line(point, snap_radius, expected_h_segments):
+    """Nearby points split the target line at the projected location."""
+    streets = _inj_network()
+    points = geopandas.GeoDataFrame(geometry=[point], crs="EPSG:3035")
+    augmented = neatnet.inject_points(streets, points, snap_radius=snap_radius)
+    assert (augmented["road"] == "H").sum() == expected_h_segments
+    assert (augmented["road"] == "V").sum() == 1
+    h_rows = augmented[augmented["road"] == "H"]
+    for geom in h_rows.geometry:
+        for _, y in geom.coords:
+            assert y == 0.0, "y coordinate drifted, the split should stay on the line"
+
+
+def test_inject_points_snap_radius_rejects_far_points():
+    """Points beyond snap_radius do not change the topology."""
+    streets = _inj_network()
+    points = geopandas.GeoDataFrame(
+        geometry=[shapely.Point(40, 50)],
+        crs="EPSG:3035",
+    )
+    augmented = neatnet.inject_points(streets, points, snap_radius=5)
+    geopandas.testing.assert_geodataframe_equal(
+        augmented.sort_values("road").reset_index(drop=True),
+        streets.sort_values("road").reset_index(drop=True),
+    )
+
+
+def test_inject_points_none_snap_radius_keeps_all():
+    """snap_radius=None splits regardless of distance."""
+    streets = _inj_network()
+    # Point is below line H at distance 100; projection (40, 0) is interior of H.
+    # snap_radius=None must accept it; snap_radius=10 would reject.
+    points = geopandas.GeoDataFrame(
+        geometry=[shapely.Point(40, -100)],
+        crs="EPSG:3035",
+    )
+    augmented = neatnet.inject_points(streets, points, snap_radius=None)
+    assert len(augmented) > len(streets)
+
+
+def test_inject_points_reprojects_points_to_streets_crs():
+    """Points in a different CRS are reprojected to match streets."""
+    streets_3035 = _inj_network()
+    point_in_4326 = geopandas.GeoSeries(
+        [shapely.Point(40, 0)], crs="EPSG:3035"
+    ).to_crs("EPSG:4326").iloc[0]
+    points_4326 = geopandas.GeoDataFrame(
+        geometry=[point_in_4326], crs="EPSG:4326"
+    )
+    augmented = neatnet.inject_points(streets_3035, points_4326, snap_radius=10)
+    assert (augmented["road"] == "H").sum() == 2
+
+
+def test_inject_points_requires_crs():
+    """Missing CRS on either input raises ValueError."""
+    streets = _inj_network()
+    points = geopandas.GeoDataFrame(geometry=[shapely.Point(40, 0)])
+    with pytest.raises(ValueError, match="CRS"):
+        neatnet.inject_points(streets, points)
+
+
+def test_inject_points_empty_points_returns_unchanged():
+    """Empty input points returns streets unchanged."""
+    streets = _inj_network()
+    empty_points = geopandas.GeoDataFrame(geometry=[], crs="EPSG:3035")
+    augmented = neatnet.inject_points(streets, empty_points, snap_radius=5)
+    assert len(augmented) == len(streets)
+
+
+def test_inject_points_multiple_pois_on_one_line():
+    """Two POIs on the same LineString split it into three segments."""
+    streets = geopandas.GeoDataFrame(
+        {"road": ["canal"]},
+        geometry=[shapely.LineString([(0, 0), (100, 0)])],
+        crs="EPSG:3035",
+    )
+    pois = geopandas.GeoDataFrame(
+        geometry=[shapely.Point(30, 1), shapely.Point(70, 1)],
+        crs="EPSG:3035",
+    )
+    augmented = neatnet.inject_points(streets, pois, snap_radius=5)
+    assert len(augmented) == 3