first round geometry.voronoi_skeleton() tests

README.md

@@ -1,6 +1,6 @@
 # `sgeop`: Street Geometry Processing Toolkit
 
-[![Continuous Integration](https://github.com/uscuni/sgeop/actions/workflows/testing.yml/badge.svg)](https://github.com/uscuni/sgeop/actions/workflows/testing.yml)
+[![Continuous Integration](https://github.com/uscuni/sgeop/actions/workflows/testing.yml/badge.svg)](https://github.com/uscuni/sgeop/actions/workflows/testing.yml) [![codecov](https://codecov.io/gh/uscuni/sgeop/branch/main/graph/badge.svg?token=VNn0WR5JWT)](https://codecov.io/gh/uscuni/sgeop)
 
 Street geometry processing toolkit
 
@@ -15,7 +15,6 @@ This package developed & and maintained by:
 * [Anastassia Vybornova](https://github.com/anastassiavybornova)
 * [James D. Gaboardi](https://github.com/jGaboardi)
 
-
 ## Introduction
 
 ## Documentation

codecov.yml

@@ -0,0 +1,22 @@
+codecov:
+  notify:
+    after_n_builds: 6
+coverage:
+  range: 50..95
+  round: nearest
+  precision: 1
+  status:
+    project:
+      default:
+        threshold: 2%
+    patch:
+      default:
+        threshold: 2%
+        target: 80%
+  ignore:
+    - "tests/*"
+comment:
+  layout: "reach, diff, files"
+  behavior: once
+  after_n_builds: 6
+  require_changes: true

sgeop/geometry.py

@@ -13,9 +13,7 @@
 from .nodes import consolidate_nodes
 
 
-def _is_within(
-    line: np.ndarray[shapely.Geometry], poly: shapely.Polygon, rtol: float = 1e-4
-) -> bool:
+def _is_within(line: np.ndarray, poly: shapely.Polygon, rtol: float = 1e-4) -> bool:
     """Check if the line is within a polygon with a set relative tolerance.
 
     Parameters
@@ -97,33 +95,39 @@ def angle_between_two_lines(
 
 
 def voronoi_skeleton(
-    lines,
-    poly=None,
-    snap_to=None,
-    max_segment_length=1,
-    buffer=None,
-    secondary_snap_to=None,
-    limit_distance=2,
-    consolidation_tolerance=None,
+    lines: list | np.ndarray | gpd.GeoSeries,
+    poly: None | shapely.Polygon = None,
+    snap_to: None | gpd.GeoSeries = None,
+    max_segment_length: int = 1,
+    buffer: None | float | int = None,
+    secondary_snap_to: None | gpd.GeoSeries = None,
+    limit_distance: None | int = 2,
+    consolidation_tolerance: None | float = None,
 ):
     """
     Returns average geometry.
 
     Parameters
     ----------
     lines : array_like
-        LineStrings connected at endpoints
-    poly : shapely.geometry.Polygon
-        polygon enclosed by `lines`
-    snap_to : gpd.GeoSeries
-        series of geometries that shall be connected to the skeleton
-    distance : float
-        distance for interpolation
-    buffer : float
-        optional custom buffer distance for dealing with Voronoi infinity issues
-    consolidation_tolerance : float
-        tolerance passed to node consolidation within the resulting skeleton. If None,
-        no consolidation happens
+        LineStrings connected at endpoints. If ``poly`` is passed in, ``lines``
+        must be a ``geopandas.GeoSeries``.
+    poly : None | shapely.Polygon = None
+        Polygon enclosed by ``lines``.
+    snap_to : None | gpd.GeoSeries = None
+        Series of geometries that shall be connected to the skeleton.
+    max_segment_length: int = 1
+        Additional vertices will be added so that all line segments
+        are no longer than this value. Must be greater than 0.
+    buffer : None | float | int = None
+        Optional custom buffer distance for dealing with Voronoi infinity issues.
+    secondary_snap_to : None | gpd.GeoSeries = None
+        ...
+    limit_distance : None | int = 2
+        ...
+    consolidation_tolerance : None | float = None
+        Tolerance passed to node consolidation within the resulting skeleton.
+        If ``None``, no consolidation happens.
 
     Returns
     -------
@@ -133,6 +137,8 @@ def voronoi_skeleton(
     if buffer is None:
         buffer = max_segment_length * 20
     if not poly:
+        if not isinstance(lines, gpd.GeoSeries):
+            lines = gpd.GeoSeries(lines)
         poly = shapely.box(*lines.total_bounds)
     # get an additional line around the lines to avoid infinity issues with Voronoi
     extended_lines = list(lines) + [poly.buffer(buffer).boundary]
@@ -196,9 +202,7 @@ def voronoi_skeleton(
             edgeline = shapely.intersection(edgeline, limit)
 
             # in edge cases, this can result in a MultiLineString with one sliver part
-            if edgeline.geom_type == "MultiLineString":
-                parts = shapely.get_parts(edgeline)
-                edgeline = parts[np.argmax(shapely.length(parts))]
+            edgeline = _remove_sliver(edgeline)
 
         # check if a, b lines share a node
         intersection = shapely_lines[b].intersection(shapely_lines[a])
@@ -258,15 +262,38 @@ def voronoi_skeleton(
     edgelines = edgelines[edgelines != None]  # noqa: E711
 
     edgelines = shapely.line_merge(edgelines[shapely.length(edgelines) > 0])
+    edgelines = _as_parts(edgelines)
+    edgelines = _consolidate(edgelines, consolidation_tolerance)
+
+    return edgelines, splitters
+
+
+def _remove_sliver(
+    edgeline: shapely.LineString | shapely.MultiLineString,
+) -> shapely.LineString:
+    """Remove sliver(s) if present."""
+    if edgeline.geom_type == "MultiLineString":
+        parts = shapely.get_parts(edgeline)
+        edgeline = parts[np.argmax(shapely.length(parts))]
+    return edgeline
+
+
+def _as_parts(edgelines: np.ndarray) -> np.ndarray:
+    """Return constituent LineStrings if MultiLineString present."""
     if np.unique(shapely.get_type_id(edgelines)).shape[0] > 1:
         edgelines = shapely.get_parts(edgelines)
+    return edgelines
+
 
+def _consolidate(
+    edgelines: np.ndarray, consolidation_tolerance: float | int
+) -> np.ndarray:
+    """Return ``edgelines`` from consolidated nodes, if criteria met."""
     if consolidation_tolerance and edgelines.shape[0] > 0:
         edgelines = consolidate_nodes(
             edgelines, tolerance=consolidation_tolerance, preserve_ends=True
         ).geometry.to_numpy()
-
-    return edgelines, splitters
+    return edgelines
 
 
 def snap_to_targets(edgelines, poly, snap_to, secondary_snap_to=None):

sgeop/tests/conftest.py

@@ -0,0 +1,66 @@
+import geopandas.testing
+import numpy
+import pandas
+import pytest
+import shapely
+
+line_collection = (
+    list[shapely.LineString]
+    | tuple[shapely.LineString]
+    | numpy.ndarray
+    | pandas.Series
+    | geopandas.GeoSeries
+)
+
+geometry_collection = (
+    list[shapely.GeometryCollection]
+    | tuple[shapely.GeometryCollection]
+    | numpy.ndarray
+    | pandas.Series
+    | geopandas.GeoSeries
+)
+
+
+def polygonize(
+    collection: line_collection, as_geom: bool = True
+) -> shapely.Polygon | geopandas.GeoSeries:
+    """Testing helper -- Create polygon from collection of lines."""
+    if isinstance(collection, pandas.Series | geopandas.GeoSeries):
+        _poly = geopandas.GeoSeries(collection).polygonize()
+        if as_geom:
+            return _poly.squeeze()
+        else:
+            return _poly
+    else:
+        return shapely.polygonize(collection).buffer(0)
+
+
+def is_geopandas(collection: geometry_collection) -> bool:
+    return isinstance(collection, geopandas.GeoSeries | geopandas.GeoDataFrame)
+
+
+def geom_test(
+    collection1: geometry_collection,
+    collection2: geometry_collection,
+    tolerance: float = 1e-1,
+) -> bool:
+    """Testing helper -- geometry verification."""
+
+    if not is_geopandas(collection1):
+        collection1 = geopandas.GeoSeries(collection1)
+
+    if not is_geopandas(collection2):
+        collection2 = geopandas.GeoSeries(collection2)
+
+    assert shapely.equals_exact(
+        collection1.geometry.normalize(),
+        collection2.geometry.normalize(),
+        tolerance=tolerance,
+    ).all()
+
+
+def pytest_configure(config):  # noqa: ARG001
+    """PyTest session attributes, methods, etc."""
+
+    pytest.polygonize = polygonize
+    pytest.geom_test = geom_test

sgeop/tests/test_geometry.py

@@ -1,3 +1,6 @@
+import geopandas.testing
+import numpy
+import pandas
 import pytest
 import shapely
 
@@ -73,7 +76,7 @@ def test_not_within(self):
         assert known == observed
 
 
-class TestAngelBetween2Lines:
+class TestAngleBetween2Lines:
     def setup_method(self):
         self.line1 = shapely.LineString(((0, 0), (1, 0)))
         self.line2 = shapely.LineString(((1, 0), (1, 1)))
@@ -116,3 +119,109 @@ def test_not_adjacent(self):
         ):
             observed = sgeop.geometry.angle_between_two_lines(self.line1, self.line4)
         assert observed == known
+
+
+voronoi_skeleton_params = pytest.mark.parametrize(
+    "lines_type,as_poly,buffer",
+    [
+        (list, False, None),
+        (list, True, 0.001),
+        (numpy.array, False, 0.01),
+        (numpy.array, True, 0.1),
+        (pandas.Series, False, 1),
+        (pandas.Series, True, 2.0),
+        (geopandas.GeoSeries, False, 5),
+        (geopandas.GeoSeries, True, 10.314),
+    ],
+)
+
+
+class TestVoronoiSkeleton:
+    def setup_method(self):
+        self.square = [
+            shapely.LineString(((0, 0), (1000, 0))),
+            shapely.LineString(((1000, 0), (1000, 1000))),
+            shapely.LineString(((0, 0), (0, 1000))),
+            shapely.LineString(((0, 1000), (1000, 1000))),
+        ]
+        self.known_square_skeleton_edges = numpy.array(
+            [
+                shapely.LineString(((1000, 0), (998, 2), (500, 500))),
+                shapely.LineString(((0, 0), (2, 2), (500, 500))),
+                shapely.LineString(((1000, 1000), (998, 998), (500, 500))),
+                shapely.LineString(((0, 1000), (2, 998), (500, 500))),
+            ]
+        )
+        self.known_square_skeleton_splits = [shapely.Point(0, 0)]
+        self.known_square_skeleton_splits_snap_to = [
+            shapely.Point(1000, 0),
+            shapely.Point(0, 0),
+            shapely.Point(0, 1000),
+            shapely.Point(1000, 1000),
+        ]
+
+    @voronoi_skeleton_params
+    def test_square(self, lines_type, as_poly, buffer):
+        known_edges = self.known_square_skeleton_edges
+        known_splits = self.known_square_skeleton_splits
+
+        lines = lines_type(self.square)
+        poly = pytest.polygonize(lines) if as_poly else None
+        observed_edges, observed_splits = sgeop.geometry.voronoi_skeleton(
+            lines,
+            poly=poly,
+            buffer=buffer,
+        )
+
+        pytest.geom_test(observed_edges, known_edges)
+        pytest.geom_test(observed_splits, known_splits)
+
+    @voronoi_skeleton_params
+    def test_square_snap_to(self, lines_type, as_poly, buffer):
+        known_edges = self.known_square_skeleton_edges
+        known_splits = self.known_square_skeleton_splits_snap_to
+
+        lines = lines_type(self.square)
+        poly = pytest.polygonize(lines) if as_poly else None
+        observed_edges, observed_splits = sgeop.geometry.voronoi_skeleton(
+            lines,
+            poly=poly,
+            buffer=buffer,
+            snap_to=(
+                pytest.polygonize(geopandas.GeoSeries(lines), as_geom=False)
+                .extract_unique_points()
+                .explode()
+            ),
+        )
+
+        pytest.geom_test(observed_edges, known_edges)
+        pytest.geom_test(observed_splits, known_splits)
+
+
+line_100_900 = shapely.LineString(((1000, 1000), (1000, 9000)))
+line_100_120 = shapely.LineString(((1000, 1020), (1020, 1020)))
+lines_100_900_100_120 = shapely.MultiLineString((line_100_900, line_100_120))
+line_110_900 = shapely.LineString(((1000, 9000), (1100, 9000)))
+
+
+def test_remove_sliver():
+    known = line_100_900
+    observed = sgeop.geometry._remove_sliver(lines_100_900_100_120)
+    assert observed == known
+
+
+def test_as_parts():
+    known = numpy.array([line_100_900, line_100_120, line_110_900])
+    observed = sgeop.geometry._as_parts(
+        numpy.array([lines_100_900_100_120, line_110_900])
+    )
+    numpy.testing.assert_array_equal(observed, known)
+
+
+@pytest.mark.parametrize("tolerance", [0.1, 1, 10, 100, 1_000, 10_000, 100_000])
+def test_consolidate(tolerance):
+    known = numpy.array([line_100_900, line_100_120, line_110_900])
+    observed = sgeop.geometry._consolidate(
+        numpy.array([line_100_900, line_100_120, line_110_900]), tolerance
+    )
+    numpy.testing.assert_array_equal(observed, known)

sgeop/tests/test_simplify.py

@@ -1,7 +1,7 @@
 import pathlib
 
 import geopandas
-import shapely
+import pytest
 from pandas.testing import assert_series_equal
 
 import sgeop
@@ -17,7 +17,5 @@ def test_simplify():
         geopandas.read_parquet(test_data / f"{ac}_original.parquet")
     )
 
-    assert shapely.equals_exact(
-        known.geometry.normalize(), observed.geometry.normalize(), tolerance=1e-1
-    ).all()
+    pytest.geom_test(known, observed)
     assert_series_equal(known._status, observed._status)