lgtm

pipeline/transform/postcode_boundaries/README.md

@@ -1,43 +1,43 @@
 # postcode_boundaries
 
-Synthesizes postcode boundary polygons for England and Wales from three datasets. UK postcodes don't have official boundary polygons — Royal Mail defines postcodes as sets of delivery addresses, not geographic areas. This pipeline constructs a plausible polygon for every postcode by combining Output Area boundaries, UPRN point locations, and INSPIRE cadastral parcels.
+Synthesizes postcode boundary polygons for England and Wales from three datasets. UK postcodes don't have official boundary polygons: Royal Mail defines postcodes as sets of delivery addresses, not geographic areas. This pipeline constructs a plausible polygon for every postcode by combining Output Area boundaries, UPRN point locations, and INSPIRE cadastral parcels.
 
 ## The three input datasets
 
-**1. Output Area (OA) boundaries** — ONS Census Output Areas are the smallest geographic unit in the UK census (~125 households each). They tile all of England and Wales with no gaps or overlaps. Stored in a GeoPackage in British National Grid (EPSG:27700, meters). ~190,000 OAs.
+**1. Output Area (OA) boundaries**: ONS Census Output Areas are the smallest geographic unit in the UK census (~125 households each). They tile all of England and Wales with no gaps or overlaps. Stored in a GeoPackage in British National Grid (EPSG:27700, meters). ~190,000 OAs.
 
-**2. UPRN lookup** — Every Unique Property Reference Number in England and Wales, with its grid coordinates (easting/northing in BNG), its postcode (`PCDS`), and its OA code (`OA21CD`). ~37 million rows. This is the critical bridge: it tells you which postcodes exist inside each OA, and where each address physically sits.
+**2. UPRN lookup**: Every Unique Property Reference Number in England and Wales, with its grid coordinates (easting/northing in BNG), its postcode (`PCDS`), and its OA code (`OA21CD`). ~37 million rows. This is the critical bridge: it tells you which postcodes exist inside each OA, and where each address physically sits.
 
-**3. INSPIRE Index Polygons** — Land Registry cadastral parcels covering most of England and Wales. Each ZIP contains a GML file with polygon coordinate lists representing individual land parcels (buildings, plots of land). ~24 million polygons. These give fine-grained building/plot outlines that are much more precise than anything you could derive from point locations alone.
+**3. INSPIRE Index Polygons**: Land Registry cadastral parcels covering most of England and Wales. Each ZIP contains a GML file with polygon coordinate lists representing individual land parcels (buildings, plots of land). ~24 million polygons. These give fine-grained building/plot outlines that are much more precise than anything you could derive from point locations alone.
 
 ## The four phases
 
 ### Phase 1: Loading data
 
-**OA boundaries** (`oa_boundaries.py`): Opens the GeoPackage via SQLite, reads every row from `OA_2021_EW_BGC_V2`. Each row's `SHAPE` column is a GeoPackage binary blob — a standard 8-byte header, then a variable-size envelope (bounding box), then WKB geometry. `parse_gpkg_geometry` reads byte 3 to extract the envelope type (0-4), looks up the envelope size, skips past the header, and hands the remaining WKB bytes to Shapely. Single-polygon MultiPolygons are unwrapped. Result: `dict[oa_code, Polygon]`, all in BNG.
+**OA boundaries** (`oa_boundaries.py`): Opens the GeoPackage via SQLite, reads every row from `OA_2021_EW_BGC_V2`. Each row's `SHAPE` column is a GeoPackage binary blob: a standard 8-byte header, then a variable-size envelope (bounding box), then WKB geometry. `parse_gpkg_geometry` reads byte 3 to extract the envelope type (0-4), looks up the envelope size, skips past the header, and hands the remaining WKB bytes to Shapely. Single-polygon MultiPolygons are unwrapped. Result: `dict[oa_code, Polygon]`, all in BNG.
 
-**UPRNs** (`uprn.py`): The raw parquet has far more columns than needed. The lazy scan selects only four columns, filters out Scotland (OA codes starting with `S`), drops nulls and blank postcodes (stripping whitespace first), then sorts by OA code. The sort uses `sink_parquet` to write to a temp file — this avoids polars doubling memory from an in-memory sort on ~37M rows.
+**UPRNs** (`uprn.py`): The raw parquet has far more columns than needed. The lazy scan selects only four columns, filters out Scotland (OA codes starting with `S`), drops nulls and blank postcodes (stripping whitespace first), then sorts by OA code. The sort uses `sink_parquet` to write to a temp file. This avoids polars doubling memory from an in-memory sort on ~37M rows.
 
-After reading the sorted file back, it builds an offset dictionary. Rather than grouping into Python lists (which would create 37M Python string objects), it detects group boundaries by comparing each row's OA code to the previous row's. The result is `offsets[oa_code] = (start_row, end_row)` — a simple slice into the DataFrame. The OA column is then dropped since it's no longer needed, saving ~400MB.
+After reading the sorted file back, it builds an offset dictionary. Rather than grouping into Python lists (which would create 37M Python string objects), it detects group boundaries by comparing each row's OA code to the previous row's. The result is `offsets[oa_code] = (start_row, end_row)`: a simple slice into the DataFrame. The OA column is then dropped since it's no longer needed, saving ~400MB.
 
 `get_oa_uprns` later retrieves a single OA's data by slicing `df[start:end]` and extracting the coordinates and postcodes.
 
 ### Phase 2: INSPIRE data
 
-INSPIRE comes as ~350 ZIP files, each containing a GML file with thousands of `PREDEFINED` elements. Each element has a `posList` — a flat string of coordinate pairs.
+INSPIRE comes as ~350 ZIP files, each containing a GML file with thousands of `PREDEFINED` elements. Each element has a `posList`: a flat string of coordinate pairs.
 
 **Parsing** (`inspire.py:parse_inspire_zip`): Uses `iterparse` for streaming XML parsing (constant memory per ZIP). For each `PREDEFINED` element, extracts the `posList` text, splits into floats, reshapes to Nx2. Calls `elem.clear()` after each element to free XML nodes immediately.
 
 **Caching** (`inspire.py:cache_inspire`): Parsing 350 ZIPs takes a while, so results are cached as three files:
-- `inspire_coords.bin` — flat binary dump of all float64 coordinate pairs, streamed to disk as each ZIP is parsed
-- `inspire_bboxes.npy` — (N, 4) array of `[min_e, min_n, max_e, max_n]` per polygon
-- `inspire_offsets.npy` — (N, 2) array of `[byte_offset_into_coords_bin, n_points]`
+- `inspire_coords.bin`: flat binary dump of all float64 coordinate pairs, streamed to disk as each ZIP is parsed
+- `inspire_bboxes.npy`: (N, 4) array of `[min_e, min_n, max_e, max_n]` per polygon
+- `inspire_offsets.npy`: (N, 2) array of `[byte_offset_into_coords_bin, n_points]`
 
-Pre-allocates numpy arrays at 25M capacity and grows by 1.5x if needed (using in-place `resize` with `refcheck=False`). This avoids Python list overhead for 24M polygons. The coords file is written sequentially — each polygon's raw bytes are appended, and its byte offset is recorded.
+Pre-allocates numpy arrays at 25M capacity and grows by 1.5x if needed (using in-place `resize` with `refcheck=False`). This avoids Python list overhead for 24M polygons. The coords file is written sequentially: each polygon's raw bytes are appended, and its byte offset is recorded.
 
-**Loading** (`inspire.py:load_inspire`): Bboxes and offsets are loaded into RAM (~1.1GB). Coords are memory-mapped — the OS pages them in on demand from the ~3GB file, never loading the whole thing.
+**Loading** (`inspire.py:load_inspire`): Bboxes and offsets are loaded into RAM (~1.1GB). Coords are memory-mapped: the OS pages them in on demand from the ~3GB file, never loading the whole thing.
 
-**Candidate retrieval** (`inspire.py:InspireIndex`): A uniform 1km grid index is built once over the 24M parcel bboxes (`build_inspire_index`). Each OA lookup (`InspireIndex.candidates`) gathers parcels from the cells its bounding box covers plus a small overflow list of parcels larger than one cell, then applies the exact bbox overlap test — O(cells + candidates) instead of an O(24M) scan per OA (the old linear scan was ~4h of the run on its own). The candidate set and order are identical to the scan. Typically matches 10-500 parcels per OA, materialized as Shapely Polygon objects by reading their coordinate slice from the memory-mapped file; invalid polygons are repaired with `make_valid`.
+**Candidate retrieval** (`inspire.py:InspireIndex`): A uniform 1km grid index is built once over the 24M parcel bboxes (`build_inspire_index`). Each OA lookup (`InspireIndex.candidates`) gathers parcels from the cells its bounding box covers plus a small overflow list of parcels larger than one cell, then applies the exact bbox overlap test: O(cells + candidates) instead of an O(24M) scan per OA (the old linear scan was ~4h of the run on its own). The candidate set and order are identical to the scan. Typically matches 10-500 parcels per OA, materialized as Shapely Polygon objects by reading their coordinate slice from the memory-mapped file; invalid polygons are repaired with `make_valid`.
 
 ### Phase 3: Processing OAs
 
@@ -49,7 +49,7 @@ The main loop in `__main__.py` (`_process_oas`) iterates through every OA that h
 
 #### Stage A: INSPIRE-based claiming
 
-Build an STRtree spatial index over the INSPIRE candidate polygons. Convert all UPRN points to Shapely Point objects and batch-query the tree with `predicate="intersects"`. This returns pairs of (point_index, candidate_index) — which UPRNs fall inside which parcels.
+Build an STRtree spatial index over the INSPIRE candidate polygons. Convert all UPRN points to Shapely Point objects and batch-query the tree with `predicate="intersects"`. This returns pairs of (point_index, candidate_index): which UPRNs fall inside which parcels.
 
 For each INSPIRE parcel that contains at least one UPRN, run a majority vote: whichever postcode has the most UPRNs inside that parcel wins the parcel. Accumulate winning parcels per postcode, union them, and clip to the OA boundary. The result is `claimed[postcode] = polygon_within_oa`.
 
@@ -73,15 +73,15 @@ The effect: every non-parcel patch of OA gets assigned to the nearest postcode b
 
 Each postcode gets its INSPIRE-claimed polygon (if any) plus its Voronoi share (if any). These are unioned together, validated, and stripped of any non-polygonal geometry debris from `make_valid`.
 
-The output of `process_oa` is `list[(postcode, polygon)]` — the per-OA fragments. A single postcode that spans two OAs produces two separate fragments (one from each OA's processing).
+The output of `process_oa` is `list[(postcode, polygon)]`: the per-OA fragments. A single postcode that spans two OAs produces two separate fragments (one from each OA's processing).
 
 ### Phase 4: Merging and writing
 
-**Fragment merging** (`output.py:merge_fragments`): Groups all fragments by postcode, unions them. If the result is a MultiPolygon (meaning the postcode has disconnected pieces — either from spanning OAs with a gap, or algorithm artifacts), applies a 5m buffer-then-unbuffer to close tiny gaps from floating-point mismatches at OA boundary edges. If still a MultiPolygon after that, keeps the largest part **plus any other part ≥ `_MIN_DETACHED_PART_AREA` (100 m²)** (`_keep_polygon_parts`); only sub-100 m² noise slivers are dropped. Keeping substantial detached parts matters because a postcode genuinely split across an OA seam (by a railway, river, or main road wider than the 5m buffer) would otherwise lose a chunk — measured at ~1.8% of merged area left as uncovered gaps (often 3000–5000 m² building blocks) before this change.
+**Fragment merging** (`output.py:merge_fragments`): Groups all fragments by postcode, unions them. If the result is a MultiPolygon (meaning the postcode has disconnected pieces: either from spanning OAs with a gap, or algorithm artifacts), applies a 5m buffer-then-unbuffer to close tiny gaps from floating-point mismatches at OA boundary edges. If still a MultiPolygon after that, keeps the largest part **plus any other part ≥ `_MIN_DETACHED_PART_AREA` (100 m²)** (`_keep_polygon_parts`); only sub-100 m² noise slivers are dropped. Keeping substantial detached parts matters because a postcode genuinely split across an OA seam (by a railway, river, or main road wider than the 5m buffer) would otherwise lose a chunk, measured at ~1.8% of merged area left as uncovered gaps (often 3000–5000 m² building blocks) before this change.
 
-**Greenspace subtraction is connectivity-preserving** (`greenspace.py:subtract_greenspace`): park/water polygons are subtracted from each postcode, but greenspace that *crosses* a postcode (a river, a strip of parkland, a golf course through a village) would otherwise split it into scattered pieces. When the subtraction disconnects a postcode, `_reconnect_split` re-adds the narrowest removed necks — a morphological closing (`_RECONNECT_BRIDGE_M`, 25 m) clipped to the original postcode footprint — so parts ≤ ~50 m apart stay joined by a thin bridge of the postcode's own land (no address moves); genuinely wide barriers stay subtracted and the postcode legitimately splits.
+**Greenspace subtraction is connectivity-preserving** (`greenspace.py:subtract_greenspace`): park/water polygons are subtracted from each postcode, but greenspace that *crosses* a postcode (a river, a strip of parkland, a golf course through a village) would otherwise split it into scattered pieces. When the subtraction disconnects a postcode, `_reconnect_split` re-adds the narrowest removed necks, a morphological closing (`_RECONNECT_BRIDGE_M`, 25 m) clipped to the original postcode footprint, so parts ≤ ~50 m apart stay joined by a thin bridge of the postcode's own land (no address moves); genuinely wide barriers stay subtracted and the postcode legitimately splits.
 
-**GeoJSON output** (`output.py:write_district_geojson`): three passes. Pass 1 converts every postcode from BNG to WGS84 (pyproj), simplifies with 1m tolerance (Douglas-Peucker), and snaps to 6 decimal places (~0.1m precision); multi-part postcodes become `MultiPolygon` (`to_wgs84_geojson_multi`, each part handled independently), single-part stay `Polygon`. The whole set is then made a **partition** (`_resolve_overlaps`): each postcode is trimmed by the union of its higher-priority overlapping neighbours, where **priority = ascending area** (smaller postcodes win contested ground). That single rule handles both seam overlap *and* containment — an enclosed postcode is always smaller than its container, so it keeps its area while the container gets a hole (the query uses both the `overlaps` and `contains` predicates, since `overlaps` alone excludes containment). This runs last, so nothing re-introduces overlap; a postcode that would be emptied keeps its original geometry, so no active postcode is dropped. Pass 2 **de-fragments** the partition (`_eliminate_small_detached_parts`): a detached part that is *both* small in absolute terms (< `_ELIM_ABS_MAX_M2`, 3000 m²) *and* a minor fraction (< `_ELIM_FRAC_MAX`, 15%) of its postcode is absorbed into the neighbouring postcode it shares the most boundary with — the classic GIS *eliminate*. This removes the Voronoi/INSPIRE/seam *scatter* that left ~1/3 of postcodes non-contiguous, while a genuine bisection (two substantial parts split by a river/railway) keeps both parts. The land is **reassigned**, never dropped, so the output stays a gapless partition and coverage is conserved; the largest part of every postcode is always retained, so no active postcode is dropped (a tiny neighbour-less sliver in removed greenspace is dropped, a larger isolated patch is kept). Pass 3 groups postcodes by district (the outward code, e.g. `SW1A` from `SW1A 1AA`), rounds coordinates to 6dp, and writes a `{district}.geojson` FeatureCollection. Each Feature has `postcodes` (formatted like `"SW1A 1AA"`) and `mapit_code` (no space: `"SW1A1AA"`) in its properties.
+**GeoJSON output** (`output.py:write_district_geojson`): three passes. Pass 1 converts every postcode from BNG to WGS84 (pyproj), simplifies with 1m tolerance (Douglas-Peucker), and snaps to 6 decimal places (~0.1m precision); multi-part postcodes become `MultiPolygon` (`to_wgs84_geojson_multi`, each part handled independently), single-part stay `Polygon`. The whole set is then made a **partition** (`_resolve_overlaps`): each postcode is trimmed by the union of its higher-priority overlapping neighbours, where **priority = ascending area** (smaller postcodes win contested ground). That single rule handles both seam overlap *and* containment: an enclosed postcode is always smaller than its container, so it keeps its area while the container gets a hole (the query uses both the `overlaps` and `contains` predicates, since `overlaps` alone excludes containment). This runs last, so nothing re-introduces overlap; a postcode that would be emptied keeps its original geometry, so no active postcode is dropped. Pass 2 **de-fragments** the partition (`_eliminate_small_detached_parts`): a detached part that is *both* small in absolute terms (< `_ELIM_ABS_MAX_M2`, 3000 m²) *and* a minor fraction (< `_ELIM_FRAC_MAX`, 15%) of its postcode is absorbed into the neighbouring postcode it shares the most boundary with: the classic GIS *eliminate*. This removes the Voronoi/INSPIRE/seam *scatter* that left ~1/3 of postcodes non-contiguous, while a genuine bisection (two substantial parts split by a river/railway) keeps both parts. The land is **reassigned**, never dropped, so the output stays a gapless partition and coverage is conserved; the largest part of every postcode is always retained, so no active postcode is dropped (a tiny neighbour-less sliver in removed greenspace is dropped, a larger isolated patch is kept). Pass 3 groups postcodes by district (the outward code, e.g. `SW1A` from `SW1A 1AA`), rounds coordinates to 6dp, and writes a `{district}.geojson` FeatureCollection. Each Feature has `postcodes` (formatted like `"SW1A 1AA"`) and `mapit_code` (no space: `"SW1A1AA"`) in its properties.
 
 ## Memory architecture
 
@@ -97,7 +97,7 @@ The pipeline is designed to run in <12GB:
 | Fragments | Python list of (str, Shapely) | grows during processing |
 
 Key design choices:
-- INSPIRE coords are memory-mapped, not loaded — the OS pages in only the ~100-500 polygons needed per OA
+- INSPIRE coords are memory-mapped, not loaded: the OS pages in only the ~100-500 polygons needed per OA
 - UPRNs sorted + offset dict avoids per-OA groupby allocation
 - `sink_parquet` for the sort avoids doubling memory
 - `release_memory()` calls `gc.collect()` + glibc `malloc_trim(0)` to return freed pages to the OS between phases
@@ -105,25 +105,25 @@ Key design choices:
 
 ## Key invariants
 
-1. **No two postcodes cover the same ground in the output** — within an OA the INSPIRE claiming + Voronoi tile it with no overlap, and a final `_resolve_overlaps` partition pass removes the thin overlap strips that the merge buffer + per-postcode simplification introduce across OA seams (measured residual overlap ~0.01% of area)
-2. **Every postcode that exists in the UPRN data gets a polygon** — unless all its UPRNs share coordinates with another postcode's UPRNs (handled by jitter) or it has zero UPRNs
-3. **Postcode polygons never extend outside their OA(s)** — all geometry is clipped to OA boundaries
-4. **A postcode split across an OA seam keeps all its substantial parts** — `merge_fragments` keeps every part ≥ 100 m² and the output is emitted as a `MultiPolygon` (the Rust server `postcodes.rs` and `loader.py` both parse MultiPolygon); only sub-100 m² noise slivers are dropped
-5. **Postcodes are contiguous unless genuinely split** — most non-contiguity is *scatter* (a unit drawn as many disconnected specks) from point-Voronoi over sparse/interleaved UPRNs, greenspace cutting across a unit, and overlap/seam slivers. Connectivity-preserving greenspace subtraction + the `_eliminate_small_detached_parts` de-fragmentation pass absorb that scatter into neighbours (coverage-conserving), cutting the share of multi-part postcodes roughly in half (~30% → ~14% measured on the worst rural/coastal districts) without dropping any postcode or leaving coverage gaps. Genuine bisections (river/railway/major road, or a detached part above the absolute+fraction thresholds) are preserved.
+1. **No two postcodes cover the same ground in the output**: within an OA the INSPIRE claiming + Voronoi tile it with no overlap, and a final `_resolve_overlaps` partition pass removes the thin overlap strips that the merge buffer + per-postcode simplification introduce across OA seams (measured residual overlap ~0.01% of area)
+2. **Every postcode that exists in the UPRN data gets a polygon**: unless all its UPRNs share coordinates with another postcode's UPRNs (handled by jitter) or it has zero UPRNs
+3. **Postcode polygons never extend outside their OA(s)**: all geometry is clipped to OA boundaries
+4. **A postcode split across an OA seam keeps all its substantial parts**: `merge_fragments` keeps every part ≥ 100 m² and the output is emitted as a `MultiPolygon` (the Rust server `postcodes.rs` and `loader.py` both parse MultiPolygon); only sub-100 m² noise slivers are dropped
+5. **Postcodes are contiguous unless genuinely split**: most non-contiguity is *scatter* (a unit drawn as many disconnected specks) from point-Voronoi over sparse/interleaved UPRNs, greenspace cutting across a unit, and overlap/seam slivers. Connectivity-preserving greenspace subtraction + the `_eliminate_small_detached_parts` de-fragmentation pass absorb that scatter into neighbours (coverage-conserving), cutting the share of multi-part postcodes roughly in half (~30% → ~14% measured on the worst rural/coastal districts) without dropping any postcode or leaving coverage gaps. Genuine bisections (river/railway/major road, or a detached part above the absolute+fraction thresholds) are preserved.
 
 ## Module structure
 
 ```
 postcode_boundaries/
-  __init__.py         — Package docstring
-  __main__.py         — CLI entry point, four-phase orchestration
-  memory.py           — release_memory() glibc malloc_trim helper
-  oa_boundaries.py    — GeoPackage parsing, OA boundary loading
-  uprn.py             — UPRN loading (sorted DataFrame + offset dict), per-OA access
-  inspire.py          — INSPIRE GML parsing, caching, loading, bbox candidate retrieval
-  voronoi.py          — Voronoi region computation clipped to boundary
-  process_oa.py       — Per-OA processing (INSPIRE assignment + Voronoi fallback)
-  output.py           — BNG to WGS84 transform, fragment merging, GeoJSON writing
+  __init__.py         : Package docstring
+  __main__.py         : CLI entry point, four-phase orchestration
+  memory.py           : release_memory() glibc malloc_trim helper
+  oa_boundaries.py    : GeoPackage parsing, OA boundary loading
+  uprn.py             : UPRN loading (sorted DataFrame + offset dict), per-OA access
+  inspire.py          : INSPIRE GML parsing, caching, loading, bbox candidate retrieval
+  voronoi.py          : Voronoi region computation clipped to boundary
+  process_oa.py       : Per-OA processing (INSPIRE assignment + Voronoi fallback)
+  output.py           : BNG to WGS84 transform, fragment merging, GeoJSON writing
 ```
 
 Invoked as: