Stroke classification & typology example demo

neatnet/continuity.py

@@ -47,7 +47,8 @@
 def get_stroke_info(
     artifacts: geopandas.GeoSeries | geopandas.GeoDataFrame,
     roads: geopandas.GeoSeries | geopandas.GeoDataFrame,
-) -> tuple[list[int], list[int], list[int], list[int]]:
+    typify: bool = False,
+) -> tuple[list[int], list[int], list[int], list[int]] | geopandas.GeoDataFrame:
     """Generate information about strokes within ``artifacts`` and the
     resulting lists can be assigned as columns to ``artifacts``. Classifies
     the strokes within the CES typology.
@@ -62,6 +63,8 @@
         Polygons representing the artifacts.
     roads : geopandas.GeoSeries | geopandas.GeoDataFrame
         LineStrings representing the road network.
+    typify : bool = False
+        ....
 
     Returns
     -------
@@ -73,35 +76,67 @@
         Counts for 'ending' strokes - continues only at one end.
     s_ : list[int]
         Counts for 'single' strokes - does not continue.
+    roads : geopandas.GeoDataFrame
+        ...
     """
     strokes = []
     c_ = []
     e_ = []
     s_ = []
+
+    if typify:
+        roads["typology"] = "x"
+
     for geom in artifacts.geometry:
         singles = 0
         ends = 0
         edges = roads.iloc[roads.sindex.query(geom, predicate="covers")]
         ecg = edges.coins_group
+        mains_index = []
+        if typify:
+            ends_index = []
+            singles_index = []
         if ecg.nunique() == 1 and edges.shape[0] == edges.coins_count.iloc[0]:
             # roundabout special case
             singles = 1
             mains = 0
         else:
             all_ends = edges[edges.coins_end]
             ae_cg = all_ends.coins_group
-            mains = edges[~ecg.isin(ae_cg)].coins_group.nunique()
+
+            mains_slice = edges[~ecg.isin(ae_cg)]
+            mains_index += mains_slice.index.to_list()
+            mains = mains_slice.coins_group.nunique()
+
             visited = []
-            for coins_count, group in zip(all_ends.coins_count, ae_cg, strict=True):
+            for ix, coins_count, group in zip(
+                all_ends.index, all_ends.coins_count, ae_cg, strict=True
+            ):
                 if group not in visited:
                     if coins_count == (ecg == group).sum():
                         singles += 1
                         visited.append(group)
+                        if typify:
+                            singles_index += [ix]
                     else:
                         # do not add to visited -- may be disjoint within the artifact
                         ends += 1
+                        if typify:
+                            ends_index += [ix]
         strokes.append(ecg.nunique())
         c_.append(mains)
         e_.append(ends)
         s_.append(singles)
-    return strokes, c_, e_, s_
+        if typify:
+            roads.loc[mains_index, "typology"] = "C"
+            roads.loc[ends_index, "typology"] = "E"
+            roads.loc[singles_index, "typology"] = "S"
+
+    if typify:
+        group_mapper = dict(
+            roads[roads["typology"] != "x"][["coins_group", "typology"]].values
+        )
+        roads["typology"] = roads["coins_group"].map(group_mapper)
+        return roads
+    else:
+        return strokes, c_, e_, s_

neatnet/simplify.py

@@ -96,6 +96,36 @@ def _identify_non_planar(
     return artifacts
 
 
+def classification_sequence(
+    roads, artifacts, compute_coins, eps: float = 1e-4, step="singletons"
+):
+    """Classify singletons into ``{CES}`` typology."""
+
+    # Extract network nodes and relate to artifacts
+    nodes, artifacts = _link_nodes_artifacts(step, roads, artifacts, eps)
+
+    # Compute number of stroke groups per artifact
+    if compute_coins:
+        roads, _ = continuity(roads)
+
+    artifacts = _classify_strokes(artifacts, roads)
+
+    # Filter artifacts caused by non-planar intersections
+    artifacts = _identify_non_planar(artifacts, roads)
+
+    # Count intersititial nodes (primes)
+    _prime_count = artifacts["node_count"] - artifacts[["C", "E", "S"]].sum(axis=1)
+    artifacts["interstitial_nodes"] = _prime_count
+
+    # Define the type label
+    ces_type = []
+    for x in artifacts[["node_count", "C", "E", "S"]].itertuples():
+        ces_type.append(f"{x.node_count}{'C' * x.C}{'E' * x.E}{'S' * x.S}")
+    artifacts["ces_type"] = ces_type
+
+    return nodes, artifacts, roads
+
+
 def simplify_singletons(
     artifacts: gpd.GeoDataFrame,
     roads: gpd.GeoDataFrame,
@@ -159,26 +189,10 @@ def simplify_singletons(
         The road network line data following the singleton procedure.
     """
 
-    # Extract network nodes and relate to artifacts
-    nodes, artifacts = _link_nodes_artifacts("singletons", roads, artifacts, eps)
-
-    # Compute number of stroke groups per artifact
-    if compute_coins:
-        roads, _ = continuity(roads)
-    artifacts = _classify_strokes(artifacts, roads)
-
-    # Filter artifacts caused by non-planar intersections
-    artifacts = _identify_non_planar(artifacts, roads)
-
-    # Count intersititial nodes (primes)
-    _prime_count = artifacts["node_count"] - artifacts[["C", "E", "S"]].sum(axis=1)
-    artifacts["interstitial_nodes"] = _prime_count
-
-    # Define the type label
-    ces_type = []
-    for x in artifacts[["node_count", "C", "E", "S"]].itertuples():
-        ces_type.append(f"{x.node_count}{'C' * x.C}{'E' * x.E}{'S' * x.S}")
-    artifacts["ces_type"] = ces_type
+    # Classify artifact typology
+    nodes, artifacts, roads = classification_sequence(
+        roads, artifacts, compute_coins, eps=eps
+    )
 
     # Collect changes
     to_drop: list[int] = []