From 514d8d9c1f5239eb181b1b16e3c5408759d1feb2 Mon Sep 17 00:00:00 2001 From: jg-codes <53511569+jg-codes@users.noreply.github.com> Date: Mon, 15 Jun 2026 21:15:55 +0200 Subject: [PATCH 1/8] Add inject_points for projecting external points onto the network --- neatnet/__init__.py | 1 + neatnet/nodes.py | 76 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 77 insertions(+) diff --git a/neatnet/__init__.py b/neatnet/__init__.py index ecc7260..c09bd21 100644 --- a/neatnet/__init__.py +++ b/neatnet/__init__.py @@ -8,6 +8,7 @@ consolidate_nodes, fix_topology, induce_nodes, + inject_points, remove_interstitial_nodes, split, ) diff --git a/neatnet/nodes.py b/neatnet/nodes.py index 8ecdf5c..77e8ea4 100644 --- a/neatnet/nodes.py +++ b/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: From 85b5c12cf20473d133aaed1b6d9bd261843a3c97 Mon Sep 17 00:00:00 2001 From: jg-codes <53511569+jg-codes@users.noreply.github.com> Date: Mon, 15 Jun 2026 21:16:02 +0200 Subject: [PATCH 2/8] Add inject_points coverage for projected point snapping --- neatnet/tests/test_nodes.py | 109 ++++++++++++++++++++++++++++++++++++ 1 file changed, 109 insertions(+) diff --git a/neatnet/tests/test_nodes.py b/neatnet/tests/test_nodes.py index d3addde..4f8f0f1 100644 --- a/neatnet/tests/test_nodes.py +++ b/neatnet/tests/test_nodes.py @@ -1316,3 +1316,112 @@ 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 + + + From ee80cdbbba1fadcc4a886685197f58235fd2165c Mon Sep 17 00:00:00 2001 From: jg-codes <53511569+jg-codes@users.noreply.github.com> Date: Mon, 15 Jun 2026 21:16:07 +0200 Subject: [PATCH 3/8] Document inject_points in the API reference --- docs/source/api.rst | 1 + neatnet/tests/test_nodes.py | 3 --- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/docs/source/api.rst b/docs/source/api.rst index 6655ff7..7c8346b 100644 --- a/docs/source/api.rst +++ b/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 ----------------------- diff --git a/neatnet/tests/test_nodes.py b/neatnet/tests/test_nodes.py index 4f8f0f1..12dd456 100644 --- a/neatnet/tests/test_nodes.py +++ b/neatnet/tests/test_nodes.py @@ -1422,6 +1422,3 @@ def test_inject_points_multiple_pois_on_one_line(): ) augmented = neatnet.inject_points(streets, pois, snap_radius=5) assert len(augmented) == 3 - - - From 23678cf608be4d63925117ddb78cc446d6882d89 Mon Sep 17 00:00:00 2001 From: jg-codes <53511569+jg-codes@users.noreply.github.com> Date: Tue, 16 Jun 2026 17:56:56 +0200 Subject: [PATCH 4/8] Address inject_points review comments - Raise on CRS mismatch instead of implicitly reprojecting, matching the behaviour of `neatify` and letting users fix mismatched inputs explicitly (review: check-and-raise rather than silent reprojection). - Replace the separate `shapely.distance` pass with `sindex.nearest(..., return_distance=True)`, which returns the point-to-line distances directly (one geometry pass instead of two). - Run ruff-format on the CRS test (fixes the failing pre-commit.ci check). Co-Authored-By: Claude Opus 4.8 --- neatnet/nodes.py | 19 +++++++++++-------- neatnet/tests/test_nodes.py | 18 +++++++++--------- 2 files changed, 20 insertions(+), 17 deletions(-) diff --git a/neatnet/nodes.py b/neatnet/nodes.py index 77e8ea4..b5c261d 100644 --- a/neatnet/nodes.py +++ b/neatnet/nodes.py @@ -359,16 +359,15 @@ def inject_points( 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. + radius are ignored. 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. + Point features to project onto ``streets``. Must share the same CRS + as ``streets``. snap_radius : float | None = None Maximum projection distance, in ``streets`` CRS units. Points beyond this distance from any LineString are dropped from the @@ -390,16 +389,20 @@ def inject_points( 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) + raise ValueError( + "The input `streets` and `points` data are in " + "different coordinate reference systems. Reproject and rerun." + ) 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]) + # Nearest line per input point, with distance (sindex-accelerated) + (input_idx, line_idx), distances = streets.sindex.nearest( + point_geoms, return_all=False, return_distance=True + ) if snap_radius is not None: keep = distances <= snap_radius diff --git a/neatnet/tests/test_nodes.py b/neatnet/tests/test_nodes.py index 12dd456..31803a8 100644 --- a/neatnet/tests/test_nodes.py +++ b/neatnet/tests/test_nodes.py @@ -1380,17 +1380,17 @@ def test_inject_points_none_snap_radius_keeps_all(): assert len(augmented) > len(streets) -def test_inject_points_reprojects_points_to_streets_crs(): - """Points in a different CRS are reprojected to match streets.""" +def test_inject_points_crs_mismatch_raises(): + """Mismatched CRS raises rather than implicitly reprojecting.""" 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" + point_in_4326 = ( + geopandas.GeoSeries([shapely.Point(40, 0)], crs="EPSG:3035") + .to_crs("EPSG:4326") + .iloc[0] ) - augmented = neatnet.inject_points(streets_3035, points_4326, snap_radius=10) - assert (augmented["road"] == "H").sum() == 2 + points_4326 = geopandas.GeoDataFrame(geometry=[point_in_4326], crs="EPSG:4326") + with pytest.raises(ValueError, match="coordinate reference system"): + neatnet.inject_points(streets_3035, points_4326, snap_radius=10) def test_inject_points_requires_crs(): From f58d5684691333302e1abf4d7da6ea7db84f450d Mon Sep 17 00:00:00 2001 From: jg-codes <53511569+jg-codes@users.noreply.github.com> Date: Tue, 16 Jun 2026 21:13:54 +0200 Subject: [PATCH 5/8] DOC: scope inject_points docstring to near-line use, document limitations Clarify that inject_points projects to the foot of the perpendicular and adds a node *on* the existing geometry (no kink), unlike split() which snaps the line to the point; the two coincide only when the point already lies on the line. Narrow the stated scope to near-line features and document the limitations raised in review: point attributes are not carried over, the projection clamps to endpoints (which can be a no-op), and a genuinely off-network point needs a connector edge rather than inject_points. Co-Authored-By: Claude Opus 4.8 --- neatnet/nodes.py | 50 ++++++++++++++++++++++++++++++++++-------------- 1 file changed, 36 insertions(+), 14 deletions(-) diff --git a/neatnet/nodes.py b/neatnet/nodes.py index b5c261d..c3258a3 100644 --- a/neatnet/nodes.py +++ b/neatnet/nodes.py @@ -354,37 +354,59 @@ def inject_points( 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. + """Project external points onto the nearest line and split it at the foot. + + Each input point within ``snap_radius`` of the network is projected onto + its *nearest* LineString in ``streets`` -- at the foot of the + perpendicular, via linear referencing -- and that line is split at the + projected location. The new node therefore lands exactly *on* the original + geometry, so the line's shape (and hence its length and any along-line + referencing) is left unchanged. Points farther than ``snap_radius`` from + every line are ignored; no node is injected for them. + + This makes ``inject_points`` the geometry-preserving way to turn near-line + features (e.g. a structure surveyed a few metres off the digitised + centreline) into exact topological nodes. It differs from :func:`split`, + which snaps the *line* onto the raw point (via ``shapely.snap``) and so + bends the geometry toward an off-line point; the two coincide only when the + point already lies on the line. Parameters ---------- streets : geopandas.GeoDataFrame LineString network. CRS must be set. points : geopandas.GeoDataFrame - Point features to project onto ``streets``. Must share the same CRS - as ``streets``. + Point features to project onto ``streets``. Must share the same CRS as + ``streets`` (a mismatch raises rather than reprojecting implicitly). + Only the geometries are used; point attributes are not carried over. 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. + Maximum projection distance, in ``streets`` CRS units. Points beyond + this distance from every line are ignored (no node injected). Choose a + value that reflects genuine on-network membership. ``None`` (the + default) injects *every* point, however far off-network it lies. 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. + ``streets`` with the relevant LineStrings split at the projected + locations. Point attributes are not propagated. + + Notes + ----- + The projection clamps to a line's endpoints: a point beyond the end of its + nearest line projects onto that endpoint, which -- if it is an existing + node -- injects no new node. Because the projection is the foot of the + perpendicular, ``inject_points`` is not a substitute for a connector + (spur) edge: it does not link a genuinely off-network point -- one whose + true location is not meant to lie on the network -- back to the line. For + that case, add an explicit connector edge to the projected node instead. See Also -------- + split : split-at-points primitive; snaps the line to the point (kink). 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.") From 56c0a493bc79911c4685d7514575c6e4aca58cbd Mon Sep 17 00:00:00 2001 From: jg-codes <53511569+jg-codes@users.noreply.github.com> Date: Fri, 19 Jun 2026 18:27:49 +0200 Subject: [PATCH 6/8] DOC: contrast split (snaps line through point) with inject_points in docstrings Per review: make the distinction explicit in both functions' docstrings. `split` now states that snapping forces the line through an off-line point (a kink) and points to `inject_points`; `inject_points` notes it keeps the original linework, just subdivided. Co-Authored-By: Claude Opus 4.8 --- neatnet/nodes.py | 26 ++++++++++++++++++++------ 1 file changed, 20 insertions(+), 6 deletions(-) diff --git a/neatnet/nodes.py b/neatnet/nodes.py index c3258a3..7890272 100644 --- a/neatnet/nodes.py +++ b/neatnet/nodes.py @@ -53,7 +53,15 @@ def split( *, eps: float = 1e-4, ) -> gpd.GeoSeries | gpd.GeoDataFrame: - """Split lines on new nodes. + """Split lines at the given points, snapping each line through the point. + + Each point is snapped onto the nearest line within ``eps`` (via + ``shapely.snap``) and that line is split there. Because the snap moves the + line onto the point, a point that does not already lie on the line forces + the line to pass through it, introducing a vertex (a kink) at the point. + Use :func:`inject_points` instead when points lie off the line and you want + to keep the original linework, splitting it at the projected foot rather + than bending it toward the point. Parameters ---------- @@ -64,12 +72,17 @@ def split( crs : str | pyproj.CRS Anything accepted by ``pyproj.CRS``. eps : float = 1e-4 - Tolerance epsilon for point snapping. + Tolerance epsilon for point snapping. Points within ``eps`` of a line + are snapped onto it; the line is moved to pass through them. Returns ------- geopandas.GeoSeries | geopandas.GeoDataFrame Resultant split line geometries. + + See Also + -------- + inject_points : project off-line points onto the nearest line (no kink) and split. """ split_points = gpd.GeoSeries(split_points, crs=crs) for split in split_points.drop_duplicates(): @@ -366,10 +379,11 @@ def inject_points( This makes ``inject_points`` the geometry-preserving way to turn near-line features (e.g. a structure surveyed a few metres off the digitised - centreline) into exact topological nodes. It differs from :func:`split`, - which snaps the *line* onto the raw point (via ``shapely.snap``) and so - bends the geometry toward an off-line point; the two coincide only when the - point already lies on the line. + centreline) into exact topological nodes: the original linework is kept, + just subdivided. It differs from :func:`split`, which snaps the *line* onto + the raw point (via ``shapely.snap``) and so bends the geometry toward an + off-line point. The two coincide only when the point already lies on the + line. Parameters ---------- From a57c0b407a44c27afc4b283c011928da583bf056 Mon Sep 17 00:00:00 2001 From: jg-codes <53511569+jg-codes@users.noreply.github.com> Date: Fri, 19 Jun 2026 18:33:06 +0200 Subject: [PATCH 7/8] DOC: add inject_points vs split demo to the simple preprocessing guide Show on a small example that split snaps the line through an off-line point (a kink) while inject_points projects the point onto the line and subdivides it, keeping the original geometry. Addresses the review request to make the distinction clear in the user guide notebook. Co-Authored-By: Claude Opus 4.8 --- docs/source/simple_preprocessing.ipynb | 70 +++++++++++++++++++++++++- 1 file changed, 69 insertions(+), 1 deletion(-) diff --git a/docs/source/simple_preprocessing.ipynb b/docs/source/simple_preprocessing.ipynb index 619d83b..22a6a6e 100644 --- a/docs/source/simple_preprocessing.ipynb +++ b/docs/source/simple_preprocessing.ipynb @@ -524,8 +524,76 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The figures above are self-explanatory. However, remember that the extended network is not topologically correct and is not suitable for network analysis directly. Use `induce_nodes` to fix it if needed.\n", + "The figures above are self-explanatory. However, remember that the extended network is not topologically correct and is not suitable for network analysis directly. Use `induce_nodes` to fix it if needed.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Inject points\n", + "\n", + "Sometimes you need a node at a specific external location, for example a point of interest sitting near (but not on) the network. `neatnet.split` and `neatnet.inject_points` both add a node, but they treat an off-line point differently.\n", + "\n", + "`split` snaps the line onto the point, so a point that is not already on the line forces the line to bend through it. `inject_points` projects the point onto the nearest line and splits it at that projection, leaving the original linework unchanged. Both require the network and the points to share a CRS." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [ + { + "output_type": "display_data", + "metadata": {}, + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAABHMAAACWCAYAAABdE8RVAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlcelbwAAAAlwSFlzAAASdAAAEnQB3mYfeAAALxRJREFUeJzt3Qm0XEWZwPFKyIIkGHgJECAEhLAMyCKSREEdkJ0gso0aRJBBGYGDw6AOMDqiE2ckLmBUQEDUkZHgsESRALIvhiUIAgYUErYQCVsCISwhCfScfzE33O7XfXt/vbz/75x3ktd93+17627VX31VNSCXy+WCJEmSJEmSOsLAVm+AJEmSJEmSKmcwR5IkSZIkqYMYzJEkSZIkSeogBnMkSZIkSZI6iMEcSZIkSZKkDmIwR5IkSZIkqYMYzJEkSZIkSeogBnMkSZIkSZI6iMEcSZIkSZKkDmIwR5IkSZIkqYN0bDDn2muvDd/73vfC888/31br6kbdUj6l9qNb9q8Rfvvb38ayePnll9tiPe3qkUceifs3e/bs0J/df//9sRz4txFeeumluL4rr7yyKeVdy3pKbVMj9cVnSOo7v/zlL8MZZ5xR8evdug3Vuuuuu8LPfvaz8P3vfz9cdtlloVu88cYb8R5/ySWXNGR9ybPszjvvDN3u4osvjvu6YsWKrvqsVkjq5q+88kqrN0WNlOtQRx99dI7N/9Of/tRW66rUzJkzc9/97ndzixYtyrW7RpZPK/e71H604vi3q49//OOxLB5//PG2WE+jPPTQQ/G8u+eeexqyvksuuSTu37e//e1cp+5DI5xzzjmxHPi3EebOnRvX98lPfrIp5V3LekptUyP1xWdI6jvbb799brXVVqv49Z/97Ge5adOmtXQbWu2QQw6J98HkZ++99851ixdffDHu0+67717R8tSdeN7/4Q9/yHyWTZkyJdftdtlll7ivS5curXtdl156aSzX1157remf1Y6SuvlTTz3V9LJU3xkUOtTee+8d1lprrbDuuuu21boqddFFF4Vf/epXYf/99w89PT2hnTWyfNpxv1tx/NW37rnnnvCVr3wlnHnmmWHHHXese31bbrll+NKXvhQmTpwYOnUfOkmjyrsVx02SEkceeWR45plnehXI1KlTw4IFC8IXv/jFlm1DK11zzTUxE2e99dYLBx54YBg+fHjYaqutQn/117/+NT7vv/rVr4Zddtml1ZvTNc4///zw+9//Pnz2s58N73rXu0J/w7U1bty4sOaaa9a9rv5elu2kY4M5//AP/xB/2m1d3ajby6fb90+Nt+2228ZUVXVWeXvcJLXSv/zLv7T8ALTDNhRKuurSWDF58uRWb47UlQi8qPt0bDCHcU4eeOCB2MKwzjrr5PUHnDt3bjjmmGNiZP+GG24IDz/8cMy82GuvvYpmX5RaV2LZsmXhtttuC/PmzYu/b7PNNuFDH/pQGDhwYMm+sbNmzYp9WgcPHhy23nrr8IEPfCAMGDAgvv+DH/wg/OUvf4n///nPf77qM9dee+1w9NFHV1UO6f0lMnrdddeFxx57LK6T/WWdpbzwwgvhpptuii00I0aMiNH/zTbbrGllXel+05fzlltuCU899VRYY401whZbbBFb0pPya+dzqZpzhmNFBeYzn/lMbI1K3HvvveHGG28M66+/fvj0pz+dt17KbenSpb1a7yo9lul9GjZsWPwb9mn06NHhkEMOySyn6dOnh7/97W/xXGZfapW1Hvb9T3/6U3jttdfC2LFjw0c/+tFeLQiUGWW35557hu233z7cfffdMWuFa23XXXfttd/s81VXXRX/z/6uXLly1Xuf+9zn4vGsFtf2FVdcET7ykY+ECRMmrHqd6/6OO+4Ihx56aNhkk03y/mbx4sVxLADuB/vtt1/ee2+99Vb4wx/+EK/d5cuXh/e85z3hgx/8YDzvGrEP5dZfz7anlTsWaS+++GJs1eHcpaVojz32qKq8b7311jj2DdvDdhXzwx/+MKy22mrh+OOPL7meWrep2vO23s+Q1HgLFy6M4448/fTTYdSoUeG9731vfF6nnXvuueHNN98Mxx13XBwHjuv32WefjfdJnkNDhw6t+PMYr4Zr/6STToq/L1q0KD7XuccyTkc6aE1mClnMjVa4DcX2k3sZ2TI8qzfccMOYwUydIWvcr5tvvjlmF62++uph/Pjx8flcDvW8X//617FMQd2HzwRBHT478fjjj8fnGPdQ6iw8Y4rVw9L78eqrr8Y63vz588NOO+2Ul+VCfZRnCHVOnonUR9797ncX3c4///nP4Y9//GNcdoMNNoj3+FL1a+po7A/1sY033jjWF6vB9ibjBXFups8J6oTUDQtRj6POyTlEfblY9m7y/GM/qX+R/XP77bfH8qSskuwKtptj+dxzz8UMep6XPNeKZWbwefxtoaz6BObMmRM/m3o9dRGuO37n5+CDDw6bbrpp0bLh85LvOVyvlG0lWf783bRp08ITTzwRfz/77LPj9wtwPA877LC6P4tzjXoJ5yn1Ds7/ar67MF4P188///M/x/OXc4jrg+PNZ2dl01R6zJLvAF/4whfy6n/pz0bWftdSlmqyXIcqNc5J0h/wtttuy+2www55/W/XWGON3BVXXFHxunD55ZfnRo8enbcefrbeeuvcnDlzivYhXH/99XstTx/lZPyQoUOH9nqfn80226zqckj297rrrsttvvnmeetba621cr/97W+L/t13vvOd3Oqrr563/IABA3JHHnlk7vXXX29KWVey3z/96U/jdhcus9VWW+XuvffeVcsx5g59Nc8+++ymjZlTy7lUzTnDtvM6+5zGMUg+44033lj1On14Bw0aFPv01nosk326/vrr4/YU9k0vNtbNm2++mTv++OPj6/RpL1xnMdWu55FHHsmNHz++V5lxLvzqV78qOk7LD37wg1Wfk/wMHDgwd8oppxTdlmI/jFXSyLFXvvrVr8bXr7766l5/85e//CW+9+lPfzrvdc7rLbbYote2rb322rkLLrig7n2oZP21bnu1xyJ9jYwYMaLXNX7ttddWPGZOsuwRRxxR9DNmzZrVa11ZY+ZUu03Vnre1foak5vjSl76UGzJkSK/r9yMf+UjuhRdeWLXchhtumBs5cmR8bhbWTzbddNPcn//855rHq+FvS93b0/eDZMw07iHNGjMn2c/bb7+9Vx2G9x588MGinzN16tTcsGHDem3/XnvtlXv++eczt/Gmm24quf/UvUA96HOf+1ys16Tfp075rW99q+R+3HjjjblRo0atWp5nHJ588sncbrvtVvS+/ctf/jJvXYwtsuuuu/Zalv0tVv/kHFlnnXXylt14443jvlQ6Zk5SJy32c8cdd+Q9y/7jP/4jd+KJJ/Za7qijjuq13uRvvvnNb8byTC/PcXrrrbdiGQ0ePDjvPc6VE044IbdixYq89bFflHMxpeoTrINtK9zez3/+87mvf/3r8f+/+93vio5jw3iBhXUZjtnNN99ctkypP5cq04kTJ9b9Wb/4xS9yPT09vda900475R577LGy25f+bMZJGjNmTN561l133XitFKr2mJUaM6ea/a60LNV3OjYzpxyi10uWLIlR4ZEjR8ZoKVF4si+IdGa1MCRoCSdbgajqbrvtFjNEiJYSbX7wwQdjSwX/kgmRRDzprpPL5WJUmYgsn8MyRNeJYhKhJsWVZdmeo446alXEkwhoPalztBixvbTSE6EnCs/20FJO1DtxzjnnhH/913+N/ycivt1228VWC1oD/vu//ztGXRnXptFlXW6/KZ9/+qd/imW88847xy4RZCDQepCU+fve9764LNFn+hOT0XLsscfWXG6N3L9qz5ndd989/g0ZP+nMJFqmaB3ieHLeEGEHn0t5JH9Xz7Fk22k9ogWEc/Lv/u7viu7766+/HqPsv/nNb2I2ECnQpTLSsmSthxaFD3/4w7GlkxZRWvSI8tPSRFmQuUTr3N///d/nrfO73/1ubE0ls4EMEI4JZXT66afHlowke4k+wpQDx4Zjkm6xyspc6yuHH354vFZpGaSlkdaxRx99NLZQsf//+I//WNc+VLL+elV6LEAW3Cc/+cm4P7Q801JKyyDnLPeFSnEdjBkzJrZgnnXWWXmtTOD8T871cmrZpmrP20btt6T6zZgxI86YNGjQoLDPPvvE+xbZJWRgcO9idkue9wkyVahP8bwnW4ZsR57dtF7vu+++8bleKqsjC/UfxvEqlnWb1HfSY6ZNmjQpHHTQQaFZ2M8DDjggPid4dgwZMiRcffXVMVuGLFCeG2lTpkwJX//612M2zsc+9rGYDcB+XH/99fHeRr2JjIFSNtpoo7j/1FnJkqH+lGRycH8HdTwyQ8l2IBOKbBeOE9vyta99Ld53C7uOsR98Nu9RhyOTh4wUMqC4J1Pf5FjyTGSbqTexDWSFc+8G9T6W5RiTRUm9lOwInp9kLpCRwvON5zPIyPj4xz8eMzQ4n6i7kcnDsp/61KcqPgbUEzkXebaRQZPOJiLzIe2nP/1p3Hb2g3rck08+GY8X5xPn9Sc+8Yle6z/vvPPi8WTf+G7AsaOcvvWtb4X//M//jMvwHs8p9pXz/Ec/+lGs15KRUQ/GAGLbqP9xLKmXPPTQQzHLJznepXDec1x5jnKt8Zxl+zheSTZMKZzHnGeUKcc+nYnE+VTPZ1144YXxOxj3ErJYOPZk55NFRjYXx4Hnf6UZfHwmdWaOHceF64dt5rrkvE9vb6OPWSX7XW1Zqg/kOlS5bApaOxcuXLjqdaKTH/7wh+N7v/nNb8quiywCMkZovU5nhCRoceZvzjjjjPj7ypUrY5Sa177yla/0iobed999uSeeeGLV77Rwsywt3vVI9peWiHT0d/ny5bnDDjusV+sOLRxJq8GPf/zjvHXRcpC0UrG9WeVTa1ln7feMGTPie8cee2yv9/7617/mtQo999xzsVWNFoZmZ+ZUun/VnjMg+r7eeuvltfSzzGmnnRYzG2ipSJx00knxvSRCXsuxTPaJ1sRio9mnM2pomdx5551jaxjZP9WoZj1f+MIXSs6IdMMNN8RyoIWvMBuE9RWeX2Sa8B7lks5quvDCC+PrZ555Zq4RGpWZQ3ZSco4tW7as1+wXtCzWsw/VrL+ezJxqjsWhhx66quUwfZ/k3OdaKLxnZZX3qaeeGl+nVaxwv2lNonWZe3O59dSyTdWet7V8hqTmSO536ezExK233pqXUUL9imU/9KEP5ZYsWZKXIUzLO++RNVNPVsyWW24Zsz1KocWcOk+x7W1kZg77csABB+TNUMN+JveoBQsWrHr96aefjvWMbbbZJq+OBO75++23X/wbyrOck08+uejz5+GHH47PFz6ncD1JdjP3+ldeeaXXfuy55555r6c/533ve1+vOhDPxFtuuWXV71/72tfissWyf2bPnh0zonfcccdez4SDDjoo73k7f/78WOeqNDMHlEM6m6hQ8izjOVOYJU59kPd45hT7G35+/etf99p3ssF57+KLL85778orr4yfwzlDVlOtmTmLFy+OZcZ6eEamkaGebFupzJwJEybEdaSf89tuu+2qbJZKkIWeZCIVU+1n8Tp1nA022CCeq4VZM8ccc0z8m8KMr6zP5l7AtZXgWpw0aVJ8j+z2eo5Zucycasq4XFmq71TfxN4hiFQSiU8QMaWVAUQZy7nvvvvickRYidQz3gvZBGeccUb8SfpAkjmRtJoQESfKTWs0n5dG38lmRixpnSDCnWDMCsaLIIJKv2eyhUCUmBYn+g4n40gkaLUgMwbJ+Bx9UdaJZHwNWjXIZimchSY9LgZj29CPmBahZqt0/6o9Z0C/a1r2ibaDaDpoAdxhhx1i61aC94jSk4FT77H893//98xWECLwtAbxGWT20CJYi0rWc/nll8cWMloykjJLyo1xSMjgSpdZglYKWsHSyDKh9Yxy4TPbHS1itLSRrUFLTBoZdmThtPP6qz0WtGbTWst5zDFO3yc333zz8G//9m81DeaXZOEkyACkZZPrNKu1rp5tqua8bfR+S6pPUt8gA7YQGXfFMqW5dtPZN1zj3/nOd6quM9WCbEzqPI3IpCyHFv307DTsJ/f4wjrP7373uzgGGxnIjLmRvgf++Mc/XpXZVOz5Xamk/sp+c1zSyNjhGcO9vjBjCGxHYRY+921ccMEFvepAPBOTTGiQecDzg+do4T2e7C0yw7nXM05msq3U89j3dBYG2UdkMDUDz10yooo9F0vVv8lqLczYYX/IyCCLhKyMNLLBeC0Zx6VW1It5XrIu6r5pZKeThVQuAzidicxxSTKeqvmuUYlKPyvJ4uM5zj0gXe/nX/6u2muAbJv0uEhci1yThfeZZhyzvixjNU7XdrNKdytKJAOlFX6pKYZuCeBLSNaXQtLkkQx0y8Omlm4o9SLFvxAPUgI8pP2znTyQ+XKNUjdNAgXcNJLl+qKsE3Sr+vKXvxxTnwlccLNnv0gf5b1WqXT/qj1nwD4yGCH7yz6S0kgFgW4bvMeDgZRlKgukadKNhQAd6jmW6fTtYqi4kR7MYHmFFYVqlFsPlTC6zFUywwYPpfQX82LnfPI6FTvSP6notTsGjyOFlYolwRW6Z7LdVCoLg8LtuP5qjgUp7nx5IlCZdE9NKzYocRa+RHD+k4ZMMD0JmCfBnUpmbqhlm6o9b7nmG7nfkupDoJeBd7l+6XZCPYOGEZ6xxQYO5V7J9Vvq3pcMBtrp6EZUbP+z6jx0WeOnkjpPtcrVc3idZ0xhPYd6UuHA+ASF6DJFQK5cHQgMFMv9m3ppFrpj0cBIdye61hZ2hWrmPb5Y/ZQAFj+l6t/F9r2S+iSTV1Tz3aBQco28//3vL/o+119W0KNR3zWaUe8nUMVPI66BYvUpvsvxnY56DucxQcNmHLO+LGM1TtcGc5IvvMUkWSpZktliuKjSkfpiXybS6yzXCtwspb6UJa8n21duOwuX74uyLowKM34EX/65qZPlRB9ybkr/+7//W7ZPbTNUun/VnjNIj5vDKPLMVJTMfMB7tMQl0Xc+Kz1eTj3HstgXysIgzEUXXRSmTp0a++HWMhZAJetJyoxW0HJjmxTOCFDpOd8u0jNQFbas0eJBNgnHmocvGWdUCOlbTr/yejRi/aW2vdZjUW75ahCw4V5Bv3UyFBnLhiwYKoyFs9JkqWabaj1vG7nfkmpH9jKZFIxrwRgjNMDwrCLoStYCAeF0dgWNdMUa6nj+co232/OmGfWdUnUexhAie7qUema/rLWeQ1Cq8HixDD+V1NPJpuSHoAgz/2QhcyJZd1/f42upfxer/zX6u0FWnaHUZ5Q7Lo38rtHoej91/lKNWqgkeFhNfSp9z2nX73PqO9YgS0haeBlUMz0tYCnJgG133XVXRQXf6Gm26eJTOA0kGR18ieNhlKTNke6ZTKVbDN3F0ss1WiX7TWtKukWFtEIqCwyyRZCnW84ZEJwiPZOIPmXPdKFJwIYsLyqcBHoI5iCdmtrMY0nXMrbrm9/8ZvyyT5pmLVN4l1sP2WJUuqhAf+Mb3+g1iG25cz7r9fR+N2ta+0JJ5Z9sj0JkyJVCqx6DS/IDWg85D+hux8DCdM+pZx8qWX+t217NseA+lAwSTCtPOpUfpKxXixTgE088MWa4Ecz5n//5n9iaWsnAx7VuU7XnbTP2W1L9CDQkwQa+mNGgRLdH6iDpbtx0J2Jg98IA8f333x+/5NRbZ+qrZ1Qz6jwMvJt0N2u0dD0n6d5eaz2H4A5ZR2SIcC/OCkCxLPUzsm3I3kpPkV4KDSRMisEztHAK62rv8X19PtRSn6TOwEDKhVnTpeoMSRkmwwoUIvu82Rpdrsk1QEZ9pfX+cqg3FZ5vNFIxJAPlnwQpO+H7nPpG146ZUy/S1hgnhUyJc889t+gyBEqStEFSKFn+7rvvLnpBL1iwICxcuHDV71TswcXZCHxpXrRoUa9xUeifype25OJnO8mOYNT+Sy+9NG95Zo1iVhjUmw1QStZ+8/mUUyHGXCFanNyYuuWcSXB8CLyRhZT8DoJwdIkhmEP3K4Ig6VmMmn0s+ZL67W9/O84wQRCp8PxqxHo4L+l+RUWZQEPS97wwhblYt7VLLrkkzhKWRuCP8qLlKd1K0ujrrVxQtzDlnONLP+hi+8Y9oxCz4bEu3k/6KdeyD9Wsv9ptr+VYUOGjqxfr5LxIo9vSf/3Xf4VqsX6yj0iJJ92eoA5BUGZPq0Qt21TteduM/ZZUO+6LXKOFLdlJVm2x+sYpp5wSr/kE9Stm52lEnYn7O+MF0i25U3AP5N7GGDFkOJUqZ+57tUrKlToVM4alMTtmMpZgesanLMnYP8ccc0yv48+xJTiXSGapYjyXpEEtjeOV7hZEFz0CGyeffHLM6knwOYX3/XL6qs6S4Lynnk2ZUt8sPIY8V/nyzj4mqDMQAL3yyivzludcKHwt+QyuMRpcCgM31BeoOzdbo8uVOi0NO4yvRD2oGILA6e9/5Zx22ml51wznUjJrbfo+U8sxa6S+PkdVmpk5JfBlgAGs+EJAiiXT/9F6Q59FItFz5swJs2bNil1/uKGxPF/GSflnkFcuapbnCznLzpw5M2YlJINaJVNBM3BtMt0lrbfpKaqrwRdkxlxJpibnpsj28aXj1FNPXbUcLcLcFGjBJpWYhxUZPbQ+MIBdMphWMshuo2XtN4P90tWI7ji0flFW7Bc3Kr4spfty0vLBVJX8fbOnJm/WOZMgePOTn/wkDszHoIzpgbJ5j6wWcKzSrR99cSypvPI5ZD4wpgDHiBaIRq6HgQFJc2f8AgJTVBBpEeMYJ1OAUqaUZ+GXeNbFIG/JdNjJw5RpE9PZD8l5RzlT+Ug+my/iSaYQU2NS6SKjg0yWWjClKK1VBNd4gJJdxTnMg55BiAsxcB5BOc53AnecF1QoGQOGChHnd3KuVLIP9ay/2m2v9VhwP+J405JLRZgKCf3J+dxy3blK4b7L+cO9l1Y/ptdMTytcTi3bVO1524z9llQbGjt4FnOvo+szmRRkVfBasbEjuDfyhYnum0xFTjYOAwCTgcC9plxXnHK4v9PCznOcwC+fxzTDZCWDeyr1SO6vzZyavBpM180znXEOuZ/RRZxxhfiSx9geZKoTgGHMDr7w1mK77baLZUI9kEYBnjE8s6hP8WzCSSedVLbreLouQv2ILsccd4I7ydTk1NHZh+SeTYYW92dep17GdrAsGZlkt3KP58s82waeP3Rj5u+TMQ4J+NBAwr/VIGuIL+Ksj8xPsn74nWnW04PjNgp1HrLfGaeRMmBad85JnmU846gDkBmVHiaARhS6NJMde8QRR8RsEqYZp8y4nqh/pFFnYRgF6lo0fibHMqkvcB1R52hmxkdSj6LOwLGnbkLZVtr4U4jzjsaYE044IdbDCSqyb9TJuJ+QZUNDJt/JKj1uXC+c95Qv1xJlTHCZe0J6/KZajlk7l6XqkOtQ5aaTZkrkUlNfT5kypaJ1genkmGo6mTIv/TNu3LjcH//4x15TAjK9ZOGyG220Ud60dUytzVR26WWY1rpayf5efvnluXXXXTdvfUOHDo3T/RViCu0vfvGLcbrHwu3cY4898qala3RZZ+339ddfH8upWFlvt912uUcffbTXVMnpab2bNTV5NftXyznD1N3JsWAawzSmyUz+9oc//GHdxzJrn7Le/8lPfhI/gykT//a3vxX923rWw/TpTG9arMyYdpS/K5wOm6lgP/jBD/Za/sgjj8ybjjrxsY99rNeyc+fOzZtms9R9oFCpKa6TfWRKyPTnMC0p020WTu/NlI7jx48vut+jRo3qNd13uX0oVO36q9n2eo7FueeeG6eaTS/LNOKXXXZZVVOTJ/iM9H2lcL8qWU+121TteVvrZ0hqvLPOOiteo8Wu3f333z+3dOnSvKmumYKZ6ZypVxXeR2+77ba6pwW/884747TN6XWn7wcXXnhhfI0pips5NXk1U00ndRDeKyyX5IepjpnavNapyfHSSy/l9t13317rpi5x7LHH5lasWFHxfuCBBx7IbbXVVr3Wx735Rz/6Ud6yjzzySMln6PDhw3Onn3563vKcI8l00ckP9UGmiq5manIcddRRvT7zjjvuyHuWFauDgu8h1LPSyv0N08nzzC62r0y3XjjN+/Lly2M9s3DZz3zmM7lTTz216PFkHcmU1umfAw44IHfCCSfE/990001Fp81OX5OJ888/P75HXaQS1El7enryPnvixIl1fxb18zXXXLNo2VFHKJy2vJjks2fOnJkbMWJE3jo4npdeemndx6zc1OTV7He5slTf6djMHFqRiXwmo2wniJzTWlAsSs/YHbQSF2YqJGmzxQadYhYYWkGIzBNxZln6H5IFUyzjgYwTos2kDDLDFZFK+l7vt99+MXMjHVGlxYKIPcuReVJsKsxKMbgW0W2yO4jqsn7KIp0BkiBbZ9q0aXFbGQCQvph016ElvliqaiPLOmu/yUKhCxKt3LSus10sT0sP2TppRPBZfzUD85baj0buXy3nDPuSdJMjWymNCD+ReFoBk/Tgeo5l1j5lvc805xwnutIwblG5Vshq10NGEcec1jJa85gtiFYc1sF5kUzvmEZrFcvTWkl/dK5fWslKzWDFOcd1SesnqeyUaTKWFP+ntY2U1UpaT7LuGewj5U+aMeukhe3ggw+OmVKcM+mZHCgLWm24dikTzn/2i9ZXppYs3O+sfSim2vVXs+3g2uR1MtBooa30WJDiTrowLZqc93T74nNIUWd9hYMFsh28TnZRMWSsMQ0o+8r/ud8Wk7WeareplvO2ls+Q1Hi0aDPdNVmiZNdw7dKqzP2v2IwuoOWdTEcGk6d7JPUrWsOL3YPJ8OR5XOnr3JPoHk22D5m8tKqn7wfUI7lHFM7QlKXabeB5XCoLk3s8n8+sOoV1ELrhkplAnYesAF4jk4X6C9lFlSCzhwzFwvWDegTPPTIaGV+QDFrqa9TdkuyASvcD1MXI7KGLFt1R6C7HvZjshsKxSqjr8Vzh+Uk3Xu7b1EtZnsybwqnPk+wM6uKUMZk83OPJTKL8qsmSYPp0vk/w+TyPed4nM2Ulz7JSmddkuRfWv8r9DfWfX/ziF3EZrgsyazi3yRRjpqlCfKfhmPPcpxz5nbo6x5L6KPWkwuNJefEeWU2UJ1k41BN4XpJJXmyMl8mTJ8dzqdjgvGSwsL3FZporhmuW64y6FPUhtjGdDV/rZ3H+U/enPFg/z3SOPddwqdm7SqH8mCWLc4jsHuqknEPFZkmr9pgldfPCTLla9rtcWarvDCCiE/o5brzcVOjTyJgnnYQLk4oFAZxigRup29DNiK5155xzTt2p7QkqdlTwWB/rLYcuNgyOyWxQlUx/LUnqTHSf5As/ExRIqh1BU8aDKQzY0AWZBqZkEOn+OLguAVOGYqB8qpkMROrYzJxGoN8mWSAEcoigdlogR1JjkFlBq1J6fKliaMmjVem8886Lv5fKFJEkSdI7yFYio5UsHMbyI2jDoNNkXpFbQAZIfwzkSPXo18GcqVOnxpRQbhzlZmvpK6SrJbPLZGmnQfCkTkfKNhk5BHXLXZ/JLGEMLFcsxVuSJEn56AJEwznd7PlJ8D2M4QLoHiapOv06mMMsMPRrZGaC9JTPrUQfWfqflkM6IsGccmOgSN0mGael0j7SlWDcg0rQR57xlZhVg9mLJEndrdwYLJIqw5hyjJfJmDmM5cfsjwR3GK+IMX36s6xxa6QsjpnTZphum0Fzy2EQPAbvkyRJkiRJ/YvBHEmSJEmSpA4ysNUbIEmSJEmSpMoZzJEkSZIkSeogBnMkSZIkSZI6iMEcSZIkSZKkDmIwR1I+pmCdNCmEAQNCGDjw7X+Tn+T3/fd/ezlJkiTVz/qXpCoZzJGU7/TTQ7jqqrf/n8vlv5f8PnNmCFOnWnKSJEmNYP1LUpWcmlxSfqvQRhuF8NxzvQM5eXeOASGst14I8+eHMHiwJShJklQr61+SamBmjqR3zJoVwrPPZgdywPvPPPP28pIkSaqd9S9JNTCYI+kdixc3d3lJkiTlyS1aZP1LUtUM5kh6R09Pc5eXJElStPjV5eH8Wx8LJ13/VFUlMuX2Z8KVDzwdlq98y5KU+jHHzJH0DvtsS5IkNU0ulwt3Pb44XHTX/HDNnGfC8jffCoPeXBnuOPuzYeRrL2W2tBO6eWHY2mHnY38eVq42KIwcNiQcutOYMHn82LDJqGEeNamfMTNH0jsYzPj44ysbM4flHPxYkiSp4iyc3b9/S/jUeXeGK+5/OgZyMHb0iPDk5M+W/WLG+w8eeHgYNvxd8fdFry4P597yWNj1ezeHw86/02wdqZ8xM0dS7+ycgw56e/pxZq1KB3aS3ydNCmHGDIM5kiRJVWThJIasNjDsu+3oMHnC2DDxPT1hwMqVFde/loWBcX2sd/YT+eMXmq0j9R8GcyQVD+hMnRrCWWe9PWtVYvTotzNyTj7ZQI4kSVKJLJzL7lkQps+eHx574dW89zZdZ1g4bMLYcPCOY0LPsCF117/mPbc0XHTXU+GyexeEJa+vyHtv581GhsMmjg17bT06DBlkhwyp2xjMkVTaihXh5auvDgNefDGsufHGIeyyi0EcSZKkerJwyLTJQlCH6cqZNZTJJiqofy1b8abZOlI/YzBHkiRJkvoyC6eJzNaR+geDOZIyLaZVKM5C7jTkkiRJDc3CaSKzdaTuZjBHUqazzz47/nvcccdZUpIkqd9qxyycSpmtI3UfgzmSMhnMkSRJ/VWnZOFUymwdqXsYzJGUyWCOJEnqbzo5C6dSZutInc1gjqRMBnMkSVJ/0G1ZOJUyW0fqTAZzJGUymCNJkrpZf8jCqZTZOlLnMJgjKZPBHEmS1G36axZOpczWkdqfwRxJkiRJ/YJZONUzW0dqTwZzJEmSJHUts3Aaw2wdqb0YzJGUad68efHfcePGWVKSJKljmIXTPGbrSK1nMEdSJsfMkSRJncIsnL5lto7UOgZzJGUymCNJktqdWTitZ7aO1LcM5kjKZDBHkiS1I7Nw2pPZOlLfMJgjKZPBHEmS1E7MwukcZutIzWMwR1ImgzmSJKnVzMLpbGbrSI1nMEdSpunTp8d/J0+ebElJkqQ+ZRZO9zFbR2oMgzmSJEmSOiYLZ5/3jg6HTRwbJr6nJwwYMKCl26r6snWunrMwTL/rqTD7icV5740cNiQcutOYMHn82LDJqGEWs1SEwRxJkiRJ7Z2FM2pYmDxhbDjk/WNCz7AhLdtGNcfcZ5eG6bOfCpfduyAseX1F3ns7bzYyBu/22np0GDJooIdA+n8GcyRlmj17dvx3woQJlpQkSWoos3CUZraOVDmDOZIyOQCyJElqNLNwVI7ZOlI2gzmSMhnMkSRJjWAWjmphto5UnMEcSZkM5kiSpHqYhaNGMVtHeofBHEmZDOZIkqRqmYWjZjJbRzKYI6kMgzmSJKlSZuGor5mto/7KzBxJma655pr47z777GNJSZKkXszCUTswW0f9jcEcSZIkSVUzC0ftymwd9QcGcyRJkiRVxCwcdRKzddTNDOZIymQ3K0mSZBaOOp3ZOuo2BnMkZXIAZEmS+iezcNSNzNZRtzCYIymTwRxJkvoXs3DUX5ito05mMEdSJoM5kiR1P7Nw1J+ZraNOZDBHUiaDOZIkdS+zcKR8ZuuoUxjMkZTJYI4kSd3FLBypPLN11O4M5kjKNHv27PjvhAkTLClJkjqYWThSbczWUTsymCNJkiR1KbNwpMYxW0ftxGCOJEmS1GXMwpGay2wdtZrBHEmZpk+fHv+dPHmyJSVJUhszC0fqe2brqFUM5kjK5ADIkiS1N7NwpPZgto76ksEcSZkM5kiS1H7MwpHal9k66gsGcyRlMpgjSVL7MAtH6ixm66hZDOZIymQwR5Kk1jILR+p8Zuuo0QzmSMpkMEeSpNYwC0fqTmbrqBEM5kjKNG/evPjvuHHjLClJkprMLByp/zBbR/UwmCNJkiS1mFk4Uv9mto6qZTBHkiRJagGzcCQVMltHlTKYIymTY+ZIktRYZuFIqoTZOspiMEdSJoM5kiTVzywcSbUyW0fFGMyRlMlgjiRJtTMLR1Ijma2jhMEcSZkM5kiSVB2zcCQ1m9k6MpgjKZPBHEmSKmMWjqRWMFunfzKYIymTwRxJkkozC0dSuzBbp38xmCMp0+LFi+O/PT09lpQkScnz8dXl4bJ7FoTps+eHx154Na9cNh01LEyeMDYc8v4xoWfYEMtMUp8zW6f7GcyRJEmSKmAWjqROY7ZO9zKYIymTmTmSpP7OLBxJ3cBsne5iMEdSJsfMkST1R2bhSOpWZut0B4M5kjIZzJEk9Sdm4UjqT8zW6VwGcyRlMpgjSep2ZuFI6u/M1uk8BnMkZTKYI0nqVmbhSFJvZut0BoM5kjIZzJEkdROzcCSpMmbrtDeDOZIyGcyRJHUDs3AkqXZm67QfgzmSJEnqSmbhSFJjma3TPgzmSJIkqauYhSNJzWe2TmsZzJGUad68efHfcePGWVKSpLZlFo4ktYbZOq1hMEdSJsfMkSS1M7NwJKl9mK3TdwzmSMpkMEeS1IlZOJMnjA0f2LQnDBgwoKXbKkn9OVuH+/TdT7yY997IYUPCoTuNCZPHjw2bjBrWsm3sdAZzJGUymCNJahdm4UhS5zFbpzkM5kjKZDBHktRKZuFIUncwW6exDOZIymQwR5LUCmbhSFL3MlunfgZzJGWaPn16/Hfy5MmWlCSpqczCkaT+xWyd2hnMkSRJUkuZhSNJMlunOgZzJEmS1OfMwpEkFWO2TmUM5kjKNHv27PjvhAkTLClJUt3MwpEkVcpsndIM5kjK5ADIkqR6mYUjSaqH2Tq9GcyRlMlgjiSpVmbhSJIazWydtxnMkZTJYI4kqRpm4UiS+sKyFW+Gq+csDBfdNT/c/cSLee+NHDYkHLrTmDB5/NiwyahhXXlADOZIymQwR5JUCbNwJEmtMvfZpWH67KfCZfcuCEteX5H33s6bjQyHTRwb9tp6dBgyaGDoFgZzJGUymCNJKsUsHElSO1nWj7J1DOZIynTNNdfEf/fZZx9LSpIUmYUjSWp3c7s8W8dgjqSSVqxYEWbNmhUWL14cenp6wi677BIGDx5siUlSP2QWjiSpEy3r0mwdgzmSigZxTj/99HDWWWeFZ599dtXro0ePDscdd1w45ZRTDOpIUj9hFo4kqVvM7aJsHYM5knoFcg488MBw1VVXhQEDBsSW2FU3jP//fdKkSWHGjBkGdCSpS5mFI0nqZsu6IFunrYI5ixYtit05JLUO2TjTpk0ru9yJJ54Ys3QktQZdH0eOHGnxq6F1sCXL3gzXzl0Srnp4SViwJL/FcsyIwWG/LUeEvTYfEUasvpolL0nqCk+++EZ87l039+WwdPlbee/tsP67wv5brRV23nh4GLzagLaqg7VVMIfuHM8//3yrN0Pq11k5e+65Z6zQZ90ayNDhBnbttdeanSO1yDrrrBPWW289y191435/zb2PhovvfirMmv9aWJmqx5JlvsvYNcLe49YM2643NN7/JUnqRm+sfCvcPv+1cPXcV8JDz7+R996IoQPDHpsND3tvPjxsv+kGbVEHa6tgzgU3Phhm3Lew1Zsh9VuL5t0f7jznyxUv/4FjvxdGjtu+qdskqbiDdlg/HP3RbSwe1e3M6x4J026Ym/fahmsOCntvvmbYfdNhZuFIkvqd+S8tD7+f90q44bFXwysF2TrnfWLLsNeO40KrDQptZOHLb4Q5z+VHwCT1ndcWVtfN8dGFi8PCd3vNSq3wgZe99tQY+247OgZzzMKRJOltY9caEj6/U084Yoe18rJ1Rq2xWnjv+sNDO2irYM767x4a3rvu0FZvhtRvLXq5J1TT0XGz9XvCSK9ZqWXPTKkRthr97nDa3u8JW6y5wiwcSZJShg4aGHbbdHj8IVvn+dfeDKsNbI8ux23Vzcoxc6TWcswcqXM4Zo4ayTqYJEmdVQdrq8ycQYMGhaFDbWmUWoXr7/DDDy87mxUxYJYbPrw9Ugyl/ohnptTI88k6mCRJnVMHa6vMHEntkZ1z0EEHhZkzZ8ZZS9K3iOT3SZMmhRkzZjiTlSRJkiS1wMBWfKik9jV48OAYqJkyZUqv9EF+53UDOZIkSZLUOmbmSJIkSZIkdRAzcyRJkiRJkjqIwRxJkiRJkqQOYjBHkiRJkiSpgxjMkSRJkiRJ6iAGcyRJkiRJkjqIwRxJkiRJkqQOYjBHkiRJkiSpgxjMkSRJkiRJ6iAGcyRJkiRJkjqIwRxJkiRJkqTQOf4PFseKyrKSWKoAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + } + } + ], + "source": [ + "import matplotlib.pyplot as plt\n", + "from shapely.geometry import Point\n", + "\n", + "streets = gpd.GeoDataFrame(geometry=[LineString([(0, 0), (10, 0)])], crs=\"EPSG:3857\")\n", + "points = gpd.GeoDataFrame(geometry=[Point(5, 2)], crs=\"EPSG:3857\")\n", + "\n", + "injected = neatnet.inject_points(streets, points, snap_radius=5)\n", + "forced = neatnet.split(points.geometry, streets, streets.crs, eps=5)\n", "\n", + "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 3))\n", + "\n", + "# inject_points: the point is projected onto the line; the linework is kept\n", + "streets.plot(ax=ax1, color=\"lightgray\", linewidth=4)\n", + "injected.plot(ax=ax1, linewidth=1.5)\n", + "points.plot(ax=ax1, color=\"red\", zorder=5)\n", + "ax1.plot([5, 5], [2, 0], color=\"gray\", linestyle=\"--\", linewidth=1)\n", + "ax1.scatter([5], [0], color=\"black\", zorder=6)\n", + "ax1.set_title(\"inject_points: linework kept, just subdivided\")\n", + "ax1.set_axis_off()\n", + "\n", + "# split: the line is snapped onto the point, bending through it\n", + "streets.plot(ax=ax2, color=\"lightgray\", linewidth=4)\n", + "forced.plot(ax=ax2, linewidth=1.5)\n", + "points.plot(ax=ax2, color=\"red\", zorder=5)\n", + "ax2.set_title(\"split: line forced through the point\")\n", + "ax2.set_axis_off()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "On the left, `inject_points` projected the point down onto the line and added a node there (black), keeping the line straight. On the right, `split` pulled the line up to the point, leaving a kink. Use `inject_points` when points lie near the line and you want to preserve its shape; points farther than `snap_radius` are ignored." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ "For more details and further options, see the [API documentation](api.rst)." ] } From 5c4f385731373469a76ba1a91b68bcc81d5d8214 Mon Sep 17 00:00:00 2001 From: jg-codes <53511569+jg-codes@users.noreply.github.com> Date: Sun, 12 Jul 2026 11:28:35 +0200 Subject: [PATCH 8/8] Address review comments on inject_points - split docstring: fix inverted snap direction -- the line is snapped onto the point, not the point onto the line; same fix in the eps description and the pre-existing _snap_n_split helper docstring - inject_points docstring: drop exactness claims; the interpolated vertex lies on the original geometry only up to floating-point rounding - pass radius to sindex.nearest as max_distance, pruning the tree search and replacing the post-hoc distance filter; the bound stays inclusive and radius=0 still means "on-line points only" (geopandas rejects max_distance=0), both pinned by new tests - rename snap_radius -> radius throughout (docstring, tests, docs notebook) --- docs/source/simple_preprocessing.ipynb | 4 +- neatnet/nodes.py | 55 +++++++++++++++----------- neatnet/tests/test_nodes.py | 40 +++++++++++++------ 3 files changed, 60 insertions(+), 39 deletions(-) diff --git a/docs/source/simple_preprocessing.ipynb b/docs/source/simple_preprocessing.ipynb index 22a6a6e..c007ed8 100644 --- a/docs/source/simple_preprocessing.ipynb +++ b/docs/source/simple_preprocessing.ipynb @@ -561,7 +561,7 @@ "streets = gpd.GeoDataFrame(geometry=[LineString([(0, 0), (10, 0)])], crs=\"EPSG:3857\")\n", "points = gpd.GeoDataFrame(geometry=[Point(5, 2)], crs=\"EPSG:3857\")\n", "\n", - "injected = neatnet.inject_points(streets, points, snap_radius=5)\n", + "injected = neatnet.inject_points(streets, points, radius=5)\n", "forced = neatnet.split(points.geometry, streets, streets.crs, eps=5)\n", "\n", "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 3))\n", @@ -587,7 +587,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "On the left, `inject_points` projected the point down onto the line and added a node there (black), keeping the line straight. On the right, `split` pulled the line up to the point, leaving a kink. Use `inject_points` when points lie near the line and you want to preserve its shape; points farther than `snap_radius` are ignored." + "On the left, `inject_points` projected the point down onto the line and added a node there (black), keeping the line straight. On the right, `split` pulled the line up to the point, leaving a kink. Use `inject_points` when points lie near the line and you want to preserve its shape; points farther than `radius` are ignored." ] }, { diff --git a/neatnet/nodes.py b/neatnet/nodes.py index 7890272..d41e16d 100644 --- a/neatnet/nodes.py +++ b/neatnet/nodes.py @@ -55,8 +55,8 @@ def split( ) -> gpd.GeoSeries | gpd.GeoDataFrame: """Split lines at the given points, snapping each line through the point. - Each point is snapped onto the nearest line within ``eps`` (via - ``shapely.snap``) and that line is split there. Because the snap moves the + Each line within ``eps`` of a point is snapped onto that point (via + ``shapely.snap``) and split there. Because the snap moves the line onto the point, a point that does not already lie on the line forces the line to pass through it, introducing a vertex (a kink) at the point. Use :func:`inject_points` instead when points lie off the line and you want @@ -72,8 +72,8 @@ def split( crs : str | pyproj.CRS Anything accepted by ``pyproj.CRS``. eps : float = 1e-4 - Tolerance epsilon for point snapping. Points within ``eps`` of a line - are snapped onto it; the line is moved to pass through them. + Snapping tolerance. A line within ``eps`` of a point is snapped onto + it, i.e. moved to pass through the point. Returns ------- @@ -135,7 +135,7 @@ def split( def _snap_n_split(e: shapely.LineString, s: shapely.Point, tol: float) -> np.ndarray: - """Snap point to edge and return lines to split.""" + """Snap edge to point and return the split parts.""" snapped = shapely.snap(e, s, tolerance=tol) _lines_split = shapely.get_parts(shapely.ops.split(snapped, s)) return _lines_split[~shapely.is_empty(_lines_split)] @@ -364,23 +364,24 @@ def inject_points( streets: gpd.GeoDataFrame, points: gpd.GeoDataFrame, *, - snap_radius: float | None = None, + radius: float | None = None, eps: float = 1e-4, ) -> gpd.GeoDataFrame: """Project external points onto the nearest line and split it at the foot. - Each input point within ``snap_radius`` of the network is projected onto + Each input point within ``radius`` of the network is projected onto its *nearest* LineString in ``streets`` -- at the foot of the perpendicular, via linear referencing -- and that line is split at the - projected location. The new node therefore lands exactly *on* the original - geometry, so the line's shape (and hence its length and any along-line - referencing) is left unchanged. Points farther than ``snap_radius`` from + projected location. The new node lands on the original geometry up to + floating-point precision -- the interpolated vertex may deviate from the + mathematical segment by rounding error -- so the line's shape is preserved + rather than bent toward the point. Points farther than ``radius`` from every line are ignored; no node is injected for them. This makes ``inject_points`` the geometry-preserving way to turn near-line features (e.g. a structure surveyed a few metres off the digitised - centreline) into exact topological nodes: the original linework is kept, - just subdivided. It differs from :func:`split`, which snaps the *line* onto + centreline) into topological nodes: the original linework is kept, just + subdivided. It differs from :func:`split`, which snaps the *line* onto the raw point (via ``shapely.snap``) and so bends the geometry toward an off-line point. The two coincide only when the point already lies on the line. @@ -393,11 +394,12 @@ def inject_points( Point features to project onto ``streets``. Must share the same CRS as ``streets`` (a mismatch raises rather than reprojecting implicitly). Only the geometries are used; point attributes are not carried over. - snap_radius : float | None = None + radius : float | None = None Maximum projection distance, in ``streets`` CRS units. Points beyond this distance from every line are ignored (no node injected). Choose a - value that reflects genuine on-network membership. ``None`` (the - default) injects *every* point, however far off-network it lies. + value that reflects genuine on-network membership. ``0`` keeps only + points that already lie on a line; ``None`` (the default) injects + *every* point, however far off-network it lies. eps : float = 1e-4 Tolerance epsilon passed to :func:`split` for the actual snap. @@ -435,15 +437,20 @@ def inject_points( if len(point_geoms) == 0 or len(line_geoms) == 0: return streets.copy() - # Nearest line per input point, with distance (sindex-accelerated) - (input_idx, line_idx), distances = streets.sindex.nearest( - point_geoms, return_all=False, return_distance=True - ) - - if snap_radius is not None: - keep = distances <= snap_radius - input_idx = input_idx[keep] - line_idx = line_idx[keep] + if radius == 0: + # geopandas rejects ``max_distance=0``: filter to on-line points + (input_idx, line_idx), distances = streets.sindex.nearest( + point_geoms, return_all=False, return_distance=True + ) + on_line = distances == 0 + input_idx = input_idx[on_line] + line_idx = line_idx[on_line] + else: + # ``max_distance`` both prunes the tree search and drops points + # with no line within ``radius`` + input_idx, line_idx = streets.sindex.nearest( + point_geoms, return_all=False, max_distance=radius + ) if len(input_idx) == 0: return streets.copy() diff --git a/neatnet/tests/test_nodes.py b/neatnet/tests/test_nodes.py index 31803a8..94af117 100644 --- a/neatnet/tests/test_nodes.py +++ b/neatnet/tests/test_nodes.py @@ -1334,17 +1334,19 @@ def _inj_network(): @pytest.mark.parametrize( - ("point", "snap_radius", "expected_h_segments"), + ("point", "radius", "expected_h_segments"), [ (shapely.Point(40, 0), 1, 2), (shapely.Point(40, 1), 5, 2), + # distance to H is exactly ``radius`` -- the bound is inclusive + (shapely.Point(40, 5), 5, 2), ], ) -def test_inject_points_projects_to_target_line(point, snap_radius, expected_h_segments): +def test_inject_points_projects_to_target_line(point, 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) + augmented = neatnet.inject_points(streets, points, radius=radius) assert (augmented["road"] == "H").sum() == expected_h_segments assert (augmented["road"] == "V").sum() == 1 h_rows = augmented[augmented["road"] == "H"] @@ -1353,30 +1355,42 @@ def test_inject_points_projects_to_target_line(point, snap_radius, expected_h_se 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.""" +def test_inject_points_radius_rejects_far_points(): + """Points beyond 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) + augmented = neatnet.inject_points(streets, points, 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.""" +def test_inject_points_zero_radius_keeps_only_on_line_points(): + """radius=0 injects nodes only for points already lying on a line.""" + streets = _inj_network() + points = geopandas.GeoDataFrame( + geometry=[shapely.Point(40, 0), shapely.Point(60, 1)], + crs="EPSG:3035", + ) + augmented = neatnet.inject_points(streets, points, radius=0) + assert (augmented["road"] == "H").sum() == 2 + assert (augmented["road"] == "V").sum() == 1 + + +def test_inject_points_none_radius_keeps_all(): + """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. + # radius=None must accept it; radius=10 would reject. points = geopandas.GeoDataFrame( geometry=[shapely.Point(40, -100)], crs="EPSG:3035", ) - augmented = neatnet.inject_points(streets, points, snap_radius=None) + augmented = neatnet.inject_points(streets, points, radius=None) assert len(augmented) > len(streets) @@ -1390,7 +1404,7 @@ def test_inject_points_crs_mismatch_raises(): ) points_4326 = geopandas.GeoDataFrame(geometry=[point_in_4326], crs="EPSG:4326") with pytest.raises(ValueError, match="coordinate reference system"): - neatnet.inject_points(streets_3035, points_4326, snap_radius=10) + neatnet.inject_points(streets_3035, points_4326, radius=10) def test_inject_points_requires_crs(): @@ -1405,7 +1419,7 @@ 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) + augmented = neatnet.inject_points(streets, empty_points, radius=5) assert len(augmented) == len(streets) @@ -1420,5 +1434,5 @@ def test_inject_points_multiple_pois_on_one_line(): geometry=[shapely.Point(30, 1), shapely.Point(70, 1)], crs="EPSG:3035", ) - augmented = neatnet.inject_points(streets, pois, snap_radius=5) + augmented = neatnet.inject_points(streets, pois, radius=5) assert len(augmented) == 3