"""Cluster nearby, roughly parallel lines into groups.
This is a port of ``optimal-interlines``' ``LineGroupSeparator``. The clustering
engine (:class:`_LineGrouper`) is shared with
:mod:`shapelib.tools.lines.interlines`; the public entry point is
:func:`group_lines`.
"""
import itertools
from collections import deque
from typing import Dict, List, Optional
import geopandas as gpd
import numpy as np
import shapely
from shapelib.tools.lines.grouping.models import _GroupedLine
# A "matacão" is a stray/blob line that sits close to many others at once; lines
# that are close to at least this many neighbors are excluded from clustering.
_MATACAO_NEIGHBOR_LIMIT = 10
class _LineGrouper:
"""Group lines by proximity.
Clusters lines via a BFS over a proximity graph, then refines the clusters by
merging groups whose lines are mutually close. Callers that need per-talhão
behavior should pass only the lines that belong to a single talhão.
Parameters
----------
lines : geopandas.GeoDataFrame
Lines to group. Expected to be in a projected CRS (metres).
distance_threshold : float, default 2.5
Maximum spacing (in the CRS unit) for two lines to be considered
neighbors.
"""
def __init__(
self,
lines: gpd.GeoDataFrame,
distance_threshold: float = 2.5,
) -> None:
self.lines = self._drop_duplicate_geometries(lines)
self.distance_threshold = distance_threshold
self._id_gen = itertools.count()
@staticmethod
def _drop_duplicate_geometries(lines: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
"""Drop exact-duplicate geometries before clustering.
Identical geometries would otherwise inflate clusters and skew the
neighbor statistics. De-duplication is cheap (a single hash pass) so it is
always applied; the last occurrence of each geometry is kept, in order of
first appearance, matching the original implementation.
"""
seen: Dict[object, int] = {}
for idx in range(len(lines)):
seen[lines.geometry.iloc[idx]] = idx
return lines.iloc[list(seen.values())].reset_index(drop=True)
def next_id(self) -> int:
return next(self._id_gen)
def _grouped_line(self, idx: int) -> _GroupedLine:
return _GroupedLine(
self.lines.geometry.iloc[idx], self.next_id(), original_index=idx
)
# ------------------------------------------------------------------
# Clustering pipeline
# ------------------------------------------------------------------
def classify_lines_by_groups(self) -> List[_GroupedLine]:
"""Build grouped lines for all input lines, cluster, then refine."""
lines = [self._grouped_line(idx) for idx in range(len(self.lines))]
clustered = self.cluster_lines(lines)
return self.refine_clustering(clustered)
def cluster_lines(self, lines: List[_GroupedLine]) -> List[_GroupedLine]:
"""Assign a ``group`` to each line via BFS over the proximity graph."""
lines_dict = {line.id: line for line in lines}
visiteds: set = set()
for line in lines:
if line.id in visiteds:
continue
visiteds.add(line.id)
to_classify = deque([line])
lines_dict[line.id].group = line.id
while to_classify:
current_line = to_classify.popleft()
neighbors = self.find_line_neighbors(current_line, lines, visiteds)
for neighbor in neighbors:
if self.ensure_neighbors(
current_line, neighbor
) and not self.check_if_is_matacao(neighbor, lines):
visiteds.add(neighbor.id)
to_classify.append(neighbor)
lines_dict[neighbor.id].group = line.id
return list(lines_dict.values())
def refine_clustering(self, lines: List[_GroupedLine]) -> List[_GroupedLine]:
"""Merge distinct groups whose lines are mutually close neighbors."""
groups_to_merge = []
for line in lines:
for other_line in self.find_line_neighbors(line, lines):
if line is other_line or line.group == other_line.group:
continue
if self.ensure_neighbors(line, other_line, threshold=0.5):
groups_to_merge.append((line.group, other_line.group))
all_groups_to_merge = {g for pair in groups_to_merge for g in pair}
for line in lines:
if line.group in all_groups_to_merge:
for group1, group2 in groups_to_merge:
if line.group == group2:
line.group = group1
break
return lines
# ------------------------------------------------------------------
# Neighbor primitives (shared with the interline generator)
# ------------------------------------------------------------------
def find_line_neighbors(
self,
current_line: _GroupedLine,
lines: List[_GroupedLine],
visiteds: Optional[set] = None,
top_closests: int = 2,
max_distance: float = 9999,
) -> List[_GroupedLine]:
"""Return up to ``top_closests`` nearest unvisited lines to ``current_line``.
Distance is the smaller of the two centroid-to-line distances, matching
the original asymmetric-robust measure. The per-candidate distances are
computed in a single vectorized GEOS call, then sorted ascending (stable)
and truncated to ``top_closests`` — equivalent to the original
incremental selection, but far cheaper at scale.
"""
if visiteds is None:
visiteds = set()
candidates = [
line
for line in lines
if line.id not in visiteds and line is not current_line
]
if not candidates:
return []
geoms = np.array([c.geom for c in candidates], dtype=object)
mids = np.array([c.middle_point for c in candidates], dtype=object)
neighbor_center = shapely.distance(current_line.geom, mids)
current_center = shapely.distance(geoms, current_line.middle_point)
min_distance = np.minimum(neighbor_center, current_center)
valid = np.nonzero(min_distance <= max_distance)[0]
if valid.size == 0:
return []
order = valid[np.argsort(min_distance[valid], kind="stable")][:top_closests]
return [candidates[i] for i in order]
def ensure_neighbors(
self,
line1: _GroupedLine,
line2: _GroupedLine,
points_to_check: int = 10,
threshold: float = 0.7,
) -> bool:
"""Confirm two lines run alongside each other.
At least ``round(points_to_check * threshold)`` sampled points of one
line must lie within ``distance_threshold`` of the other line.
"""
min_points_to_beat = round(points_to_check * threshold)
fractions = [i / points_to_check for i in range(1, points_to_check)]
line1_points = np.array(
[line1.geom.interpolate(f, normalized=True) for f in fractions],
dtype=object,
)
line2_points = np.array(
[line2.geom.interpolate(f, normalized=True) for f in fractions],
dtype=object,
)
line1_score = int(
(shapely.distance(line2.geom, line1_points) < self.distance_threshold).sum()
)
line2_score = int(
(shapely.distance(line1.geom, line2_points) < self.distance_threshold).sum()
)
return line1_score >= min_points_to_beat or line2_score >= min_points_to_beat
def check_if_is_matacao(
self,
line: _GroupedLine,
lines: List[_GroupedLine],
matacao_thr: int = _MATACAO_NEIGHBOR_LIMIT,
) -> bool:
"""True if ``line`` sits close to at least ``matacao_thr`` other lines.
Counts other lines whose nearest endpoint is within
``distance_threshold`` of ``line``; the per-line distances are computed
in a single vectorized GEOS call.
"""
others = [other for other in lines if other is not line]
if len(others) < matacao_thr:
return False
firsts = np.array([o.first_point for o in others], dtype=object)
lasts = np.array([o.last_point for o in others], dtype=object)
nearest_endpoint = np.minimum(
shapely.distance(line.geom, firsts),
shapely.distance(line.geom, lasts),
)
close_count = int((nearest_endpoint < self.distance_threshold).sum())
return close_count >= matacao_thr
[docs]
def group_lines(
lines: gpd.GeoDataFrame,
distance_threshold: float = 2.5,
) -> gpd.GeoDataFrame:
"""
Cluster nearby, roughly parallel lines into groups.
Lines are clustered by proximity: stray "matacão" lines (close to many others
at once) are excluded from a cluster, and clusters whose lines are mutually
close are merged. Exact-duplicate geometries are dropped before clustering.
Filtering lines by talhão is the caller's responsibility — pass only the lines
belonging to a single talhão to group them in isolation.
Parameters
----------
lines : geopandas.GeoDataFrame
Line geometries to group. Must be in a projected CRS (metres), since
``distance_threshold`` is expressed in CRS units. MultiLineStrings are
flattened to single LineStrings in the output.
distance_threshold : float, default 2.5
Maximum spacing (in CRS units) for two lines to be treated as neighbors.
Returns
-------
geopandas.GeoDataFrame
One row per de-duplicated input line, carrying the original attribute
columns plus:
- ``group_id`` : integer cluster id (the seed line's id within the call).
- ``line_id`` : per-call unique integer id of the line.
Examples
--------
>>> import shapelib
>>> lines = shapelib.read_file("lines.shp")
>>> lines.to_utm()
>>> grouped = group_lines(lines, distance_threshold=3.5)
>>> grouped["group_id"].nunique()
2
"""
columns = ["geometry", "group_id", "line_id"]
if lines is None or len(lines) == 0:
return gpd.GeoDataFrame(
{c: [] for c in columns},
geometry="geometry",
crs=getattr(lines, "crs", None),
)
grouper = _LineGrouper(lines, distance_threshold)
lines_reset = grouper.lines
geom_name = lines_reset.geometry.name
reserved = set(columns)
attr_cols = [
c for c in lines_reset.columns if c != geom_name and c not in reserved
]
records: Dict[str, list] = {c: [] for c in columns}
for col in attr_cols:
records[col] = []
for grouped_line in grouper.classify_lines_by_groups():
records["geometry"].append(grouped_line.geom)
records["group_id"].append(grouped_line.group)
records["line_id"].append(grouped_line.id)
source = lines_reset.iloc[grouped_line.original_index]
for col in attr_cols:
records[col].append(source[col])
return gpd.GeoDataFrame(records, geometry="geometry", crs=lines.crs)