Snapping Near-Duplicate OSM Nodes with a Tolerance Jump to heading#

#!/usr/bin/env python3
"""Snap near-duplicate nodes within a tolerance, deterministically."""
from __future__ import annotations

import logging
from collections import defaultdict
from dataclasses import dataclass, field

import numpy as np
from scipy.sparse import coo_matrix
from scipy.sparse.csgraph import connected_components
from scipy.spatial import cKDTree

logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s")
logger = logging.getLogger(__name__)


@dataclass
class SnapResult:
    """Old node id → surviving node id, plus what happened."""
    mapping: dict[int, int]
    clusters: int = 0
    merged: int = 0
    max_move_m: float = 0.0
    oversized: list[tuple[int, float]] = field(default_factory=list)


def snap_nodes(node_ids: np.ndarray,
               coords: np.ndarray,
               tolerance_m: float,
               max_cluster_diameter_m: float | None = None) -> SnapResult:
    """Cluster nodes within `tolerance_m` and collapse each cluster to one node.

    `coords` must be projected metres. Running this on degrees makes the tolerance
    mean a different distance at every latitude.
    """
    if coords.ndim != 2 or coords.shape[1] != 2:
        raise ValueError("coords must be an (n, 2) array of projected metres")
    if max_cluster_diameter_m is None:
        # A cluster wider than a few tolerances is a chain that ran away.
        max_cluster_diameter_m = tolerance_m * 3

    tree = cKDTree(coords)
    pairs = tree.query_pairs(r=tolerance_m, output_type="ndarray")
    if pairs.size == 0:
        logger.info("no node pairs within %.3f m", tolerance_m)
        return SnapResult(mapping={int(i): int(i) for i in node_ids})

    n = len(node_ids)
    adjacency = coo_matrix(
        (np.ones(len(pairs)), (pairs[:, 0], pairs[:, 1])), shape=(n, n))
    count, labels = connected_components(adjacency, directed=False)

    mapping: dict[int, int] = {}
    merged = 0
    max_move = 0.0
    oversized: list[tuple[int, float]] = []

    for label in range(count):
        members = np.flatnonzero(labels == label)
        if members.size == 1:
            idx = int(members[0])
            mapping[int(node_ids[idx])] = int(node_ids[idx])
            continue

        member_coords = coords[members]
        # Diameter check: transitive chaining can produce a cluster far wider than
        # the tolerance, and that is a merge nobody asked for.
        spread = float(np.max(np.ptp(member_coords, axis=0)))
        if spread > max_cluster_diameter_m:
            oversized.append((label, spread))
            for idx in members:                       # leave the whole cluster alone
                mapping[int(node_ids[idx])] = int(node_ids[idx])
            continue

        # Deterministic representative: the lowest node id in the cluster.
        # A centroid would move every node; the lowest id keeps one of them exact
        # and does not depend on iteration order.
        survivor_idx = int(members[np.argmin(node_ids[members])])
        survivor_id = int(node_ids[survivor_idx])
        survivor_xy = coords[survivor_idx]
        for idx in members:
            mapping[int(node_ids[idx])] = survivor_id
            if idx != survivor_idx:
                merged += 1
                max_move = max(max_move, float(np.hypot(*(coords[idx] - survivor_xy))))

    logger.info("%d cluster(s), %d node(s) merged, max move %.3f m, %d oversized skipped",
                count, merged, max_move, len(oversized))
    return SnapResult(mapping=mapping, clusters=count, merged=merged,
                      max_move_m=max_move, oversized=oversized)


def rewrite_way_refs(ways: dict[int, list[int]],
                     mapping: dict[int, int]) -> tuple[dict[int, list[int]], int]:
    """Point every way at the surviving nodes, collapsing repeats the snap created.

    Snapping two adjacent nodes of the same way to one leaves the way referencing
    the same node twice in a row — a zero-length segment that fails validity.
    """
    rewritten: dict[int, list[int]] = {}
    degenerate = 0
    for way_id, refs in ways.items():
        new_refs = [mapping.get(r, r) for r in refs]
        collapsed = [new_refs[0]]
        for ref in new_refs[1:]:
            if ref != collapsed[-1]:
                collapsed.append(ref)
        # A closed way that collapses below four nodes is no longer a ring.
        if len(collapsed) < 2 or (refs[0] == refs[-1] and len(collapsed) < 4):
            degenerate += 1
            continue
        rewritten[way_id] = collapsed
    if degenerate:
        logger.warning("%d way(s) became degenerate and were dropped", degenerate)
    return rewritten, degenerate

def test_no_node_moves_further_than_the_tolerance():
    result = snap_nodes(ids, coords, tolerance_m=0.10)
    assert result.max_move_m <= 0.10 * 3          # bounded by the diameter guard

def test_snapping_is_order_independent():
    order = np.random.permutation(len(ids))
    a = snap_nodes(ids, coords, 0.10).mapping
    b = snap_nodes(ids[order], coords[order], 0.10).mapping
    assert a == b

def test_ways_have_no_dangling_refs():
    rewritten, _ = rewrite_way_refs(ways, result.mapping)
    survivors = set(result.mapping.values())
    assert all(ref in survivors for refs in rewritten.values() for ref in refs)

def test_rings_stay_closed():
    rewritten, _ = rewrite_way_refs(ways, result.mapping)
    for way_id, refs in rewritten.items():
        if ways[way_id][0] == ways[way_id][-1]:
            assert refs[0] == refs[-1], f"way {way_id} was closed and is not now"

before = gdf.geometry.area.sum()
after = repaired.geometry.area.sum()
logger.info("total area changed by %.4f%%", 100 * (after - before) / before)