"""Merge fragmented planting rows back into continuous lines.
Drone-mapped rows are sometimes "picotadas": one physical row is split into
several ``LineString`` fragments separated by holes. This tool reconstructs the
continuous rows by bridging each hole — guided, always, by a continuous neighbor
row (the "rail") that runs alongside it.
The algorithm follows each tip's neighbor to its partner. Every fragment starts
as its own one-piece run; each round, for every dangling tip it first picks the
**rail** — the nearest line running forward past the tip — and then joins the tip
to the next tip lying **in the same lane**: at the tip's own offset from that
rail, on the same side, with the rail still its direct neighbor there too. The
destination follows from riding the rail, so a tip can only ever reach its row's
direct continuation. The nearest such partner along the rail is taken first, and
each freshly merged line becomes a rail for the next round — so the
reconstruction grows outward from solid ground until nothing reconstructable is
left.
Three properties are core to the algorithm (not optional toggles):
* **Same-lane joining.** A tip is joined to another only when its rail is the
*direct* neighbor of *both* tips — same side, same offset, with no other row
between the rail and either tip. A same-row continuation rides alongside the rail
one spacing away at both ends; a link that skips a row sits two spacings off, or
reaches the far side of a row that lies between them, and is rejected. This is
what keeps the merge from joining different rows on curves.
* **Easy-first, outside-in.** A tip is joined to its nearest in-lane partner first,
and only when a rail already runs past it. A hole in the middle of a band that has
broken across several rows has no rail yet, so it waits; once the rows around it
are rebuilt they become its rail and it is bridged in turn. The reconstruction
never starts from an estimate in the middle of a hole, and a single row's several
holes are bridged one at a time as each gains a rail.
* **A fill always needs a rail.** Holes are never bridged with a coarse straight
guess: a hole with no neighbor running alongside it is left open. Every fill in
the output is guided by a real neighbor row.
The public entry point is :func:`merge_lines`.
"""
import geopandas as gpd
import numpy as np
from shapely.geometry import LineString, Point
from shapely.strtree import STRtree
from shapelib.tools.lines.grouping.core import _LineGrouper
from shapelib.tools.lines.merge.models import _Fragment
# Output columns the tool always produces; colliding input attributes are dropped.
_RESERVED_COLUMNS = {
"geometry",
"group_id",
"n_fragments_merged",
"n_gaps_filled",
}
class _Piece:
"""A run of one or more fragments merged so far.
Holds the merged geometry, the original fragment indices that compose it (for
output attributes), and how many holes were bridged inside it. Tip positions
and outward directions are taken from a :class:`_Fragment` wrapping the current
geometry, so a freshly merged piece exposes its new ends for the next round.
"""
__slots__ = ("frag_idxs", "geom", "n_gaps_filled", "_tips")
def __init__(self, frag_idxs: list[int], geom: LineString, n_gaps_filled: int):
self.frag_idxs = frag_idxs
self.geom = geom
self.n_gaps_filled = n_gaps_filled
self._tips = _Fragment(geom)
def tip(self, at_start: bool) -> tuple[np.ndarray, np.ndarray]:
return self._tips.tip(at_start)
class _LineMerger:
"""Bridges fragment holes best-first, each fill guided by a neighbor row."""
def __init__(
self,
lines: gpd.GeoDataFrame,
distance_threshold: float,
gap_tolerance: float,
angle_tolerance: float,
max_reference_distance: float,
) -> None:
self.grouper = _LineGrouper(lines, distance_threshold)
self.lines = self.grouper.lines
self.gap_tolerance = gap_tolerance
self.angle_tolerance = angle_tolerance
self.max_reference_distance = max_reference_distance
# Same-lane cue: a tip and its partner must ride the rail at nearly the same
# offset. Tied to ``distance_threshold`` (≈ the real row spacing) rather than
# ``max_reference_distance`` (often inflated by the caller to reach a rail),
# so the gate stays below one spacing. The "rail is the direct neighbor of
# both tips" check in the pass backstops it where the real spacing is smaller.
self.offset_tolerance = 0.5 * distance_threshold
# A reference is always a *different* row; a millimetre floor on the
# gap-to-reference distance excludes a row's own merged copy (distance ~0)
# without touching any real neighbor.
self.min_reference_distance = 1e-3
# A rail must run *alongside* the tip — extend past its foot on *both* sides,
# not just ahead. A candidate that begins at the tip (nothing behind its foot)
# is the row's own collinear continuation across a small gap (the very partner
# we want to join), not a side neighbor; the back floor rejects it. The forward
# floor keeps the symmetric requirement on the ahead side.
self.min_rail_forwardspan = 1.0
self.min_rail_backspan = 1.0
# The same-lane and direct-neighbor tests are measured a little *inside* each
# tip, not at the geometric end, so a bent or hooked end ("embicada") — or a
# sharper curve right at the break — does not corrupt the offset reading. One
# spacing in is past a typical hook while still local to the tip; clamped to
# half the piece so a short fragment samples its midpoint.
self.lane_inset = distance_threshold
# Two tips may be joined only when the gap runs *outward from both* — each tip
# faces into the hole. This is what stops a tip from linking to the far end of a
# parallel partner (which shares the same rail along its whole length) or to a
# sideways neighbor. The bound is deliberately generous: a real same-row break,
# even on a sharp curve, leaves the gap within a few tens of degrees of each
# tip's heading, while a far-end or sideways link sits past 80°.
self.facing_tolerance = 60.0
# Direct-neighbor gate: the rail must be the partner's *nearest* row too, not a
# row over. Reject if another piece sits more than this factor closer to the
# partner than the rail does — a line lies between them, so the rail is the
# neighbor only at the start tip, not at the landing.
self.neighbor_margin = 1.5
# Filled per run(); the original fragments, for output attribute lookup.
self.fragments: list[_Fragment] = []
# ------------------------------------------------------------------
# Top-level driver
# ------------------------------------------------------------------
def run(self) -> tuple[list[dict], list[dict]]:
"""Return ``(merged, unmerged)`` records.
``merged`` holds one record per piece of >= 2 fragments; ``unmerged`` holds
one record per fragment that was never bridged to a neighbor. The proximity
``group`` is computed only to tag the output.
"""
grouped_lines = self.grouper.classify_lines_by_groups()
fragments: list[_Fragment] = []
for grouped_line in grouped_lines:
fragment = _Fragment(
grouped_line.geom,
original_index=grouped_line.original_index,
group=grouped_line.group,
)
fragments.append(fragment)
self.fragments = fragments
pieces = [_Piece([i], f.geom, 0) for i, f in enumerate(fragments)]
pieces = self._merge_all(pieces)
merged: list[dict] = []
unmerged: list[dict] = []
for piece in pieces:
record = self._record(piece)
(merged if len(piece.frag_idxs) >= 2 else unmerged).append(record)
return merged, unmerged
def _merge_all(self, pieces: list[_Piece]) -> list[_Piece]:
"""Bridge gaps round by round until nothing reconstructable is left."""
while True:
pieces, changed = self._merge_pass(pieces)
if not changed:
return pieces
def _commit(self, pieces: list[_Piece], candidates: list) -> tuple[list[_Piece], bool]:
"""Greedily merge the best non-conflicting candidate gaps of one round.
``candidates`` are ``(cost, a_idx, a_start, b_idx, b_start, fill, bridge)``;
the lowest cost wins, each piece takes part in at most one merge, and a bridge
is also rejected against bridges already committed this round (they are not yet
in the spatial index, so two new fills cannot be allowed to cross).
"""
candidates.sort(key=lambda c: c[0])
used: set = set()
accepted_bridges: list[LineString] = []
new_pieces: list[_Piece] = []
for _, a_idx, a_start, b_idx, b_start, fill, bridge in candidates:
if a_idx in used or b_idx in used:
continue
if any(bridge.crosses(other) for other in accepted_bridges):
continue
new_pieces.append(self._join(pieces[a_idx], a_start, pieces[b_idx], b_start, fill))
used.add(a_idx)
used.add(b_idx)
accepted_bridges.append(bridge)
if not used:
return pieces, False
new_pieces.extend(p for i, p in enumerate(pieces) if i not in used)
return new_pieces, True
def _merge_pass(self, pieces: list[_Piece]) -> tuple[list[_Piece], bool]: # noqa: C901
"""One round: follow each tip's neighbor (its rail) to its partner.
For every dangling tip the rail is chosen *first* — the nearest parallel
neighbor that runs forward past the tip. The partner is then whatever tip
lies ahead **in the same lane**: at the tip's own signed offset from that
rail, on the same side, with the rail still its direct neighbor there. The
destination is a consequence of riding the rail, not a pre-picked pair — so a
tip can only ever join its row's direct continuation, never skip across to a
neighbor. The nearest such partner along the rail is taken first; each freshly
merged piece becomes a rail next round.
The join does *not* require the two tips to be collinear: if both share the
same rail at the same offset and side, they belong to the same row and are
bridged even when one end is hooked or the row curves sharply at the break.
The offset and direct-neighbor tests are read a little inside each tip
(:meth:`_lane_offset`); the only directional rule kept is that the gap must run
outward from both tips, so a tip never links to a partner's far end.
"""
geoms = [p.geom for p in pieces]
tree = STRtree(geoms)
tips: list[tuple[int, bool, np.ndarray, np.ndarray]] = []
for piece_idx, piece in enumerate(pieces):
start_pt, start_dir = piece.tip(True)
end_pt, end_dir = piece.tip(False)
tips.append((piece_idx, True, start_pt, start_dir))
tips.append((piece_idx, False, end_pt, end_dir))
tip_pts = [Point(t[2]) for t in tips]
tip_tree = STRtree(tip_pts)
candidates = []
for ti, (a_idx, a_start, pa, da) in enumerate(tips):
rail = self._find_rail(pa, da, a_idx, geoms, tree)
if rail is None:
continue
rail_geom, proj_a, fwd_sign = rail
pa_in, off_a = self._lane_offset(pieces[a_idx], a_start, rail_geom)
out_a = pa - pa_in # de-hooked outward heading at tip A
for tj in tip_tree.query(tip_pts[ti].buffer(self.gap_tolerance)):
tj = int(tj)
if tj == ti:
continue
b_idx, b_start, pb, _ = tips[tj]
if a_idx == b_idx:
continue
proj_b = rail_geom.project(tip_pts[tj])
if fwd_sign * (proj_b - proj_a) <= 0:
continue # partner must lie ahead along the rail
lo, hi = sorted((proj_a, proj_b))
if lo <= 1e-6 or hi >= rail_geom.length - 1e-6:
continue # rail must span the gap to guide the fill
# Same lane? Compare offsets read a little inside each tip, so a hooked
# or curved end does not throw the reading off.
pb_in, off_b = self._lane_offset(pieces[b_idx], b_start, rail_geom)
if abs(off_a - off_b) > self.offset_tolerance:
continue # partner is in a different lane (different row)
# Tips must face each other: the gap runs outward from both. Using the
# de-hooked headings keeps curves while rejecting a link to a partner's
# far end (gap ~opposite B's heading) or a sideways neighbor.
gap = pb - pa
out_b = pb - pb_in
if (
self._directed_angle(out_a, gap) > self.facing_tolerance
or self._directed_angle(out_b, -gap) > self.facing_tolerance
):
continue
# The rail must be the partner's *direct* neighbor too, not a row
# over: reject if another piece sits notably closer to the partner
# than the rail does — that means a line lies between them, so the
# rail is the neighbor only at the start tip, not at the landing.
# Measured inside the tip too, for the same reason.
pb_in_pt = Point(pb_in)
rail_at_b = rail_geom.distance(pb_in_pt)
nearest_b = self._nearest_other(pb_in_pt, geoms, tree, (a_idx, b_idx))
if rail_at_b > self.neighbor_margin * nearest_b:
continue
ref_dist = rail_geom.distance(Point((pa + pb) / 2))
fill = self._fill_from_reference(pa, pb, (rail_geom, proj_a, proj_b))
bridge = LineString([tuple(pa), *fill, tuple(pb)])
if self._bridge_blocked(bridge, fill, ref_dist, geoms, tree, (a_idx, b_idx)):
continue
gap_rail = abs(proj_b - proj_a)
candidates.append((gap_rail, a_idx, a_start, b_idx, b_start, fill, bridge))
return self._commit(pieces, candidates)
def _find_rail(
self,
pa: np.ndarray,
da: np.ndarray,
exclude_idx: int,
geoms: list[LineString],
tree: STRtree,
) -> tuple[LineString, float, float] | None:
"""The nearest parallel neighbor that runs *alongside* the tip.
Returns ``(rail, proj_a, fwd_sign)``: the rail geometry, the tip's arc-length
foot on it, and the sign of the rail's parameter that advances in the tip's
outgoing direction. ``None`` if no qualifying neighbor runs alongside the tip.
A rail must extend *both* ahead of and behind the tip's foot: a true side
neighbor runs alongside the tip, bracketing it. This rejects the row's own
collinear continuation across a small gap — which sits within
``max_reference_distance`` and looks parallel, but begins *at* the tip
(nothing behind its foot) — so it is never mistaken for the rail and stolen
from the partner role.
"""
pa_pt = Point(pa)
best = None
best_dist = np.inf
for i in tree.query(pa_pt.buffer(self.max_reference_distance)):
i = int(i)
if i == exclude_idx:
continue
rail = geoms[i]
dist = rail.distance(pa_pt)
if dist < self.min_reference_distance or dist > self.max_reference_distance:
continue
proj = rail.project(pa_pt)
here = np.array(rail.interpolate(proj).coords[0])
ahead = np.array(rail.interpolate(min(proj + 5.0, rail.length)).coords[0])
behind = np.array(rail.interpolate(max(proj - 5.0, 0.0)).coords[0])
v_ahead, v_behind = ahead - here, behind - here
if np.dot(v_ahead, da) >= np.dot(v_behind, da):
fdir, fwd_sign, fwd_len, back_len = v_ahead, 1.0, rail.length - proj, proj
else:
fdir, fwd_sign, fwd_len, back_len = v_behind, -1.0, proj, rail.length - proj
norm = float(np.linalg.norm(fdir))
if norm < 1e-9 or fwd_len < self.min_rail_forwardspan:
continue # rail does not extend forward of the tip
if back_len < self.min_rail_backspan:
continue # rail begins at the tip: it is the row's own continuation,
# not a side neighbor running alongside — skip so the genuine rail wins
if self._angle_between(fdir / norm, da) > self.angle_tolerance:
continue
if dist < best_dist:
best_dist = dist
best = (rail, proj, fwd_sign)
return best
def _nearest_other(
self,
pt: Point,
geoms: list[LineString],
tree: STRtree,
exclude: tuple[int, int],
) -> float:
"""Distance from ``pt`` to the nearest piece that is not in ``exclude``."""
best = np.inf
for i in tree.query(pt.buffer(2.0 * self.max_reference_distance)):
i = int(i)
if i in exclude:
continue
dist = geoms[i].distance(pt)
if dist < best:
best = dist
return best
@staticmethod
def _signed_offset(rail: LineString, point: np.ndarray, proj: float) -> float:
"""Signed perpendicular offset of ``point`` from ``rail`` at foot ``proj``.
The sign (left/right of the rail's local forward direction) is what tells two
tips apart by *which side* of the neighbor they ride, not merely how far.
"""
here = np.array(rail.interpolate(proj).coords[0])
ahead = np.array(rail.interpolate(min(proj + 2.0, rail.length)).coords[0])
rdir = ahead - here
norm = float(np.linalg.norm(rdir))
if norm < 1e-9:
behind = np.array(rail.interpolate(max(proj - 2.0, 0.0)).coords[0])
rdir = here - behind
norm = float(np.linalg.norm(rdir))
if norm < 1e-9:
return 0.0
rdir = rdir / norm
return rdir[0] * (point[1] - here[1]) - rdir[1] * (point[0] - here[0])
def _lane_offset(self, piece: _Piece, at_start: bool, rail: LineString) -> tuple[np.ndarray, float]:
"""Tip's signed offset from ``rail``, sampled a little inside the piece.
Reading the offset one ``lane_inset`` in from the geometric end (clamped to
the piece's midpoint) skips a hooked or bent tip, so the value reflects the
row's real lane rather than the stray end. Returns the sampled interior point
and its signed offset from the rail.
"""
geom = piece.geom
inset = min(self.lane_inset, geom.length / 2.0)
along = inset if at_start else geom.length - inset
point = np.array(geom.interpolate(along).coords[0])
proj = rail.project(Point(point))
return point, self._signed_offset(rail, point, proj)
@staticmethod
def _join(
piece_a: _Piece,
a_start: bool,
piece_b: _Piece,
b_start: bool,
fill: list[tuple],
) -> _Piece:
"""Concatenate two pieces tip-to-tip across ``fill`` into one piece."""
a_coords = list(piece_a.geom.coords)
if a_start: # connecting at A's start: orient so that tip comes last
a_coords.reverse()
b_coords = list(piece_b.geom.coords)
if not b_start: # connecting at B's end: orient so that tip comes first
b_coords.reverse()
coords = a_coords + list(fill) + b_coords
return _Piece(
piece_a.frag_idxs + piece_b.frag_idxs,
LineString(coords),
piece_a.n_gaps_filled + piece_b.n_gaps_filled + 1,
)
def _record(self, piece: _Piece) -> dict:
longest = max((self.fragments[i] for i in piece.frag_idxs), key=lambda f: f.length)
return {
"geometry": piece.geom,
"group_id": longest.group,
"n_fragments_merged": len(piece.frag_idxs),
"n_gaps_filled": piece.n_gaps_filled,
"original_index": longest.original_index,
}
# ------------------------------------------------------------------
# Gap reconstruction primitives
# ------------------------------------------------------------------
def _fill_from_reference(
self,
pa: np.ndarray,
pb: np.ndarray,
reference: tuple[LineString, float, float],
) -> list[tuple]:
"""Return interior fill points for the hole from ``pa`` to ``pb``.
Maps the reference row's shape across the gap onto ``pa -> pb`` with a
similarity transform so it meets both tips exactly and follows the row's
curvature.
"""
ref, proj_a, proj_b = reference
gap_len = float(np.linalg.norm(pb - pa))
step = max(0.5, gap_len / 20.0)
n = max(int(abs(proj_b - proj_a) / step), 1)
samples = [ref.interpolate(float(d)) for d in np.linspace(proj_a, proj_b, n + 1)]
q0 = complex(samples[0].x, samples[0].y)
qm = complex(samples[-1].x, samples[-1].y)
if abs(qm - q0) < 1e-9:
return []
za = complex(pa[0], pa[1])
zb = complex(pb[0], pb[1])
scale = (zb - za) / (qm - q0)
offset = za - scale * q0
out = []
for s in samples[1:-1]:
w = scale * complex(s.x, s.y) + offset
out.append((w.real, w.imag))
return out
def _bridge_blocked(
self,
bridge: LineString,
fill: list[tuple],
ref_dist: float,
geoms: list[LineString],
tree: STRtree,
exclude: tuple[int, int],
) -> bool:
"""True if ``bridge`` may not be used: it crosses, or skips, another row.
Two reasons reject a bridge:
* It properly **crosses** another piece (``crosses`` rather than
``intersects`` ignores merely touching its own pieces at the endpoints).
* Its interior runs **too close to another row** — nearer than half the
reference distance. A correct fill rides one row-spacing away from every
other row (it follows its neighbor at that distance); a fill that passes
right through another row has skipped it, connecting tips that aren't the
row's direct continuation. Scaling the threshold to the reference distance
makes this hold at any spacing, and the followed neighbor itself sits a
full reference-distance away, so it never trips the test.
"""
for i in tree.query(bridge):
i = int(i)
if i in exclude:
continue
if bridge.crosses(geoms[i]):
return True
threshold = 0.5 * ref_dist
for point in fill:
pt = Point(point)
for i in tree.query(pt.buffer(threshold)):
i = int(i)
if i in exclude:
continue
if geoms[i].distance(pt) < threshold:
return True
return False
@staticmethod
def _angle_between(d1: np.ndarray, d2: np.ndarray) -> float:
"""Undirected angle in degrees between two directions (in ``[0, 90]``)."""
cos_val = np.clip(np.dot(d1, d2), -1.0, 1.0)
return float(np.degrees(np.arccos(abs(cos_val))))
@staticmethod
def _directed_angle(u: np.ndarray, v: np.ndarray) -> float:
"""Directed angle in degrees between two vectors (in ``[0, 180]``).
Unlike :meth:`_angle_between`, this keeps the sign of the dot product, so a
vector and its opposite are 180° apart — needed to tell a tip facing into the
gap from one facing away (a far-end link).
"""
nu = float(np.linalg.norm(u))
nv = float(np.linalg.norm(v))
if nu < 1e-9 or nv < 1e-9:
return 180.0
cos_val = np.clip(np.dot(u, v) / (nu * nv), -1.0, 1.0)
return float(np.degrees(np.arccos(cos_val)))
[docs]
def merge_lines(
lines: gpd.GeoDataFrame,
distance_threshold: float = 2.5,
gap_tolerance: float = 50.0,
angle_tolerance: float = 10.0,
max_reference_distance: float = 2.0,
) -> gpd.GeoDataFrame:
"""
Merge fragmented rows back into continuous lines, filling the holes.
Takes the input lines and returns the fully fused set: fragments of one
physical row (collinear, end-to-end, with a hole between them) are detected and
concatenated into a single ``LineString``, and every line that was not part of
a merge is passed through unchanged — so the output is one self-contained,
de-fragmented line set, not just the rebuilt rows.
Each tip is joined by following its **rail** — the nearest neighbor row running
alongside it — to the next tip in the same lane (the same offset from the rail,
on the same side, with the rail still its direct neighbor there). The fill
follows the rail's curvature, each freshly merged line becomes a rail for the
rows around it (so curvature propagates outward across holes that span several
rows, and a single row's several holes are bridged one at a time), and a hole
with no neighbor running alongside it is left open rather than guessed straight.
Filtering by talhão is the caller's responsibility — pass only the lines of a
single talhão.
Parameters
----------
lines : geopandas.GeoDataFrame
Line geometries. Must be in a projected CRS (metres), since the tolerances
are in CRS units. MultiLineStrings are flattened.
distance_threshold : float, default 2.5
Maximum spacing (CRS units) for two rows to be grouped as neighbors (used
only to tag the output with a ``group_id``).
gap_tolerance : float, default 50.0
Maximum hole length (CRS units) that may be bridged between two fragments.
angle_tolerance : float, default 10.0
Maximum angular deviation (degrees) for two fragment tips to be considered
collinear, and for a reference row to be considered parallel.
max_reference_distance : float, default 2.0
Maximum perpendicular distance (CRS units) from a tip to a row for that row
to be used as its rail — about one row spacing, so only the immediate
neighbor guides the join. Holes with no row this close are left open.
Returns
-------
geopandas.GeoDataFrame
The full fused line set — one row per merged line plus every line that was
not merged, passed through unchanged. Each row carries the attributes of
its longest contributing fragment plus:
- ``group_id`` : proximity-cluster id the row belongs to.
- ``n_fragments_merged`` : number of fragments combined (>= 2 for merged
rows, ``1`` for passed-through lines).
- ``n_gaps_filled`` : number of holes bridged in the row.
Examples
--------
>>> import shapelib
>>> lines = shapelib.read_file("rows.shp")
>>> lines.to_utm()
>>> fused = merge_lines(lines)
>>> rebuilt = fused[fused["n_fragments_merged"] >= 2] # only the merged rows
"""
base_columns = ["geometry", "group_id", "n_fragments_merged", "n_gaps_filled"]
if lines is None or len(lines) == 0:
return gpd.GeoDataFrame(
{c: [] for c in base_columns},
geometry="geometry",
crs=getattr(lines, "crs", None),
)
merger = _LineMerger(
lines,
distance_threshold,
gap_tolerance,
angle_tolerance,
max_reference_distance,
)
merged, unmerged = merger.run()
records = merged + unmerged
lines_reset = merger.lines
geom_name = lines_reset.geometry.name
attr_cols = [c for c in lines_reset.columns if c != geom_name and c not in _RESERVED_COLUMNS]
out: dict[str, list] = {c: [] for c in [*base_columns, *attr_cols]}
for record in records:
for col in base_columns:
out[col].append(record[col])
source_index = record["original_index"]
source = lines_reset.iloc[source_index] if source_index is not None else None
for col in attr_cols:
out[col].append(source[col] if source is not None else None)
return gpd.GeoDataFrame(out, geometry="geometry", crs=lines.crs)