lgtm

pipeline/transform/postcode_boundaries/output.py

@@ -1,4 +1,5 @@
 import json
+import math
 import shutil
 from collections import defaultdict
 from pathlib import Path
@@ -383,6 +384,164 @@ def _resolve_overlaps(
     return [(pc, out[i]) for i, (pc, _) in enumerate(items)]
 
 
+# A detached (non-largest) part of a postcode is absorbed into the neighbouring
+# postcode it shares the most boundary with, when the part is BOTH small in
+# absolute terms AND a minor fraction of its postcode. This removes the
+# Voronoi/INSPIRE/overlap-seam SCATTER that otherwise left ~1/3 of postcodes
+# non-contiguous (a single unit drawn as many disconnected specks), while keeping
+# genuine splits: a postcode bisected by a river/railway into two SUBSTANTIAL
+# parts has both parts above these thresholds and is left alone. The land is
+# REASSIGNED to an adjacent postcode wherever a neighbour exists, so coverage is
+# conserved and the output stays a partition; only a tiny neighbour-less sliver
+# (< _ELIM_DROP_BELOW_M2, snap/overlap debris floating in removed greenspace) is
+# dropped, and a larger isolated part (a genuine detached hamlet) is kept. The
+# largest part of every postcode is always retained, so no postcode is dropped.
+_ELIM_ABS_MAX_M2 = 3000.0      # absolute size below which a minor part reads as scatter
+_ELIM_FRAC_MAX = 0.15          # ...and only when it is < this fraction of the postcode
+_ELIM_DROP_BELOW_M2 = 200.0    # a neighbour-less sliver this small is snap debris -> drop
+# WGS84-degree distances at UK latitudes (1e-6 deg ~ 0.11 m at ~53N)
+_ELIM_SHARED_EPS_DEG = 5e-7    # ~0.05 m: thin probe whose overlap area ~ shared-edge length
+# Candidate-gather + nearest-neighbour fallback radius, in degrees so it needs no
+# per-part reprojection. The metric reach is mildly anisotropic (~2.2 m N-S,
+# ~1.3-1.4 m E-W at England's latitudes); this only bounds the rare gapped-seam
+# fallback (the dominant border-sharing path is unaffected), so the approximation
+# is harmless. Gather and accept use the SAME radius, so there is no dead band.
+_ELIM_NEAREST_MAX_DEG = 2e-5
+_M_PER_DEG_LAT = 111_320.0
+_ELIM_ITERATIONS = 2
+
+
+def _approx_area_m2(geom_deg) -> float:
+    """Metric area of a WGS84 polygon via a latitude-scaled planar approximation.
+
+    Accurate to a few percent at the few-hundred-to-few-thousand m^2 scale these
+    thresholds work at, and far cheaper than a per-part CRS transform over the
+    full ~1.8M postcodes.
+    """
+    area_deg2 = geom_deg.area
+    if area_deg2 == 0.0:
+        return 0.0
+    centroid = geom_deg.centroid
+    if centroid.is_empty:
+        lat = (geom_deg.bounds[1] + geom_deg.bounds[3]) / 2
+    else:
+        lat = centroid.y
+    return area_deg2 * _M_PER_DEG_LAT * _M_PER_DEG_LAT * math.cos(math.radians(lat))
+
+
+def _geom_parts(geom):
+    return list(geom.geoms) if geom.geom_type == "MultiPolygon" else [geom]
+
+
+def _eliminate_small_detached_parts(
+    items: list[tuple[str, Polygon | MultiPolygon]],
+) -> list[tuple[str, Polygon | MultiPolygon]]:
+    """De-fragment the partition: absorb small detached parts into the neighbour
+    they border most (see the module note above).
+
+    Runs on the de-overlapped WGS84 geometries as the last shaping step. The
+    largest part of every postcode is always retained, so no active postcode is
+    dropped; parts are moved between postcodes (or, for tiny neighbour-less
+    slivers, dropped), so total covered area is conserved to within rounding.
+    """
+    geoms: dict[str, Polygon | MultiPolygon] = {pc: g for pc, g in items}
+    for _ in range(_ELIM_ITERATIONS):
+        recs: list[tuple[str, Polygon, float]] = []
+        total: dict[str, float] = defaultdict(float)
+        largest: dict[str, float] = defaultdict(float)
+        for pc, geom in geoms.items():
+            for part in _geom_parts(geom):
+                if part.is_empty:
+                    continue
+                area = _approx_area_m2(part)
+                recs.append((pc, part, area))
+                total[pc] += area
+                if area > largest[pc]:
+                    largest[pc] = area
+
+        tree = STRtree([part for _, part, _ in recs])
+        assign: dict[int, str | None] = {}
+        for i, (pc, part, area) in enumerate(recs):
+            if area >= largest[pc]:
+                continue  # the largest part always stays -> a postcode never vanishes
+            if not (area < _ELIM_ABS_MAX_M2 and area < _ELIM_FRAC_MAX * total[pc]):
+                continue  # substantial or major-fraction part = genuine split, keep
+
+            best_pc = None
+            best_score = 0.0
+            nearest_pc = None
+            nearest_d = float("inf")
+            # Gather candidates out to the nearest-fallback radius (not just the
+            # smaller border-probe radius), so a gapped seam in the
+            # (border, nearest] band can actually be reassigned rather than dropped.
+            for j in tree.query(part.buffer(_ELIM_NEAREST_MAX_DEG)):
+                if j == i:
+                    continue
+                other_pc, other_geom, _ = recs[j]
+                if other_pc == pc:
+                    continue
+                try:
+                    score = part.buffer(_ELIM_SHARED_EPS_DEG).intersection(other_geom).area
+                except GEOSException:
+                    score = 0.0
+                # Ties broken by the lexicographically smaller postcode so the
+                # result is independent of STRtree traversal order (matches the
+                # stable tie-break in _resolve_overlaps; keeps output byte-stable
+                # across shapely/GEOS upgrades for content-hash caching).
+                if score > best_score or (
+                    score == best_score and best_pc is not None and other_pc < best_pc
+                ):
+                    best_score = score
+                    best_pc = other_pc
+                dist = part.distance(other_geom)
+                if dist < nearest_d or (
+                    dist == nearest_d
+                    and nearest_pc is not None
+                    and other_pc < nearest_pc
+                ):
+                    nearest_d = dist
+                    nearest_pc = other_pc
+
+            if best_pc is not None and best_score > 0:
+                assign[i] = best_pc  # shares a border -> absorb into that neighbour
+            elif nearest_pc is not None and nearest_d <= _ELIM_NEAREST_MAX_DEG:
+                assign[i] = nearest_pc  # gapped seam -> nearest neighbour
+            elif area < _ELIM_DROP_BELOW_M2:
+                assign[i] = None  # neighbour-less snap/overlap sliver -> drop
+            # else: keep with its own postcode (a genuine isolated patch)
+
+        if not assign:
+            break
+
+        new_parts: dict[str, list] = defaultdict(list)
+        for i, (pc, part, _) in enumerate(recs):
+            if i in assign:
+                target = assign[i]
+                if target is None:
+                    continue
+                new_parts[target].append(part)
+            else:
+                new_parts[pc].append(part)
+
+        rebuilt: dict[str, Polygon | MultiPolygon] = {}
+        for pc, parts in new_parts.items():
+            if len(parts) == 1:
+                rebuilt[pc] = parts[0]
+            else:
+                merged = safe_union(parts, grid=_OUTPUT_PRECISION_DEG)
+                if not merged.is_empty:
+                    rebuilt[pc] = merged
+        # Carry forward any postcode that contributed no parts to recs (e.g. an
+        # all-empty input geometry, which the part loop skips): never drop a
+        # postcode just because reassignment fired elsewhere in this pass.
+        for pc, geom in geoms.items():
+            if pc not in rebuilt:
+                rebuilt[pc] = geom
+        geoms = rebuilt
+
+    return list(geoms.items())
+
+
 def _round_coords(coords, ndigits=6):
     if coords and isinstance(coords[0], (int, float)):
         return [round(coords[0], ndigits), round(coords[1], ndigits)]
@@ -501,6 +660,11 @@ def write_district_geojson(
     # Remove overlap strips so the output is a clean partition.
     converted = _resolve_overlaps(converted)
 
+    # De-fragment: absorb small detached parts (Voronoi/INSPIRE/seam scatter) into
+    # the neighbour they border most, so postcodes stop shipping as many
+    # disconnected specks. Coverage-preserving; runs on the final partition.
+    converted = _eliminate_small_detached_parts(converted)
+
     by_district: dict[str, list[tuple[str, Polygon | MultiPolygon]]] = defaultdict(list)
     for pc, geom in converted:
         parts = pc.split()