Fill no-one's-land gaps between merged sub-regions (v0.17.0)

Part of #16

scripts/fill-gaps.mjs

@@ -0,0 +1,364 @@
+#!/usr/bin/env node
+/**
+ * Fill "no-one's-land" gaps left in `public/regions.geo.json` by the v0.14 coarsening merge.
+ *
+ * For each country, the gap = (country landmass) − (union of the country's sub-regions):
+ * area inside the country's authoritative landmass that no sub-region covers. Each gap
+ * fragment is assigned to the **nearest** sub-region (centroid haversine distance) and
+ * appended to it as extra polygon part(s). A fragment is only added if:
+ *
+ *   • it is not a sub-tiny sliver (≥ GAP_MIN_KM2),
+ *   • it does not fall inside an existing region's *hole* (water) — this keeps the
+ *     Aral-basin hole in KZ-KZY open,
+ *   • it is not byte-identical to a part the file already contains (polyclip sometimes
+ *     re-emits a region's own outer ring with a phantom sliver hole; once the hole is
+ *     cleaned the piece equals the existing part, and appending it would only grow the
+ *     file on every re-run),
+ *   • it has been clipped against every region of *other* countries (and against pieces
+ *     already appended this run) that touch its bbox, so it never overlaps an existing
+ *     region — no new cross-border overlaps are introduced.
+ *
+ * Holes found in gap fragments are kept only when they are real: pre-existing region
+ * holes (matched by centroid) or enclaves of another country (e.g. Lesotho inside ZAF,
+ * so the ZAF fragment stays a ring and never covers the LSO region); every other hole
+ * is a polyclip phantom and is dropped.
+ *
+ * Standard GeoJSON (`{type, coordinates}`) throughout — the same shape the file and
+ * MapLibre already use, and what @turf v7 expects.
+ *
+ * Usage:   node scripts/fill-gaps.mjs [inFile [outFile]]
+ *   inFile  default: public/regions.geo.json
+ *   outFile default: inFile (edit in place)
+ * Env:      GAP_MIN_KM2 (default 25) — gap fragments below this area are ignored slivers.
+ *
+ * Deterministic: same input → same output. Re-running on an already-filled file is a
+ * no-op (every re-discovered fragment is either a sub-threshold sliver, a phantom that
+ * byte-matches an existing part, or already covered by the enlarged union), so the
+ * result is idempotent.
+ */
+import fs from "node:fs";
+
+import difference from "@turf/difference";
+import union from "@turf/union";
+import centroid from "@turf/centroid";
+import bbox from "@turf/bbox";
+import area from "@turf/area";
+import booleanPointInPolygon from "@turf/boolean-point-in-polygon";
+import distance from "@turf/distance";
+import { feature, featureCollection, multiPolygon, polygon, point } from "@turf/helpers";
+
+const GAP_MIN_KM2 = Number(process.env.GAP_MIN_KM2 ?? 25);
+const RING_MIN_KM2 = 1e-4; // drop degenerate (numerical-noise) rings
+
+// --------------------------------------------------------------------------- helpers
+
+/**
+ * Round coordinates in place to 6 decimals (file convention). Accepts either a single
+ * ring (array of positions) or a polygon part (array of rings) and detects the depth.
+ */
+function roundCoords(coords) {
+  const isRing = Array.isArray(coords[0]) && Array.isArray(coords[0][0]) ? false : true;
+  const rings = isRing ? [coords] : coords;
+  for (const ring of rings) for (const c of ring) {
+    c[0] = Number(c[0].toFixed(6));
+    c[1] = Number(c[1].toFixed(6));
+  }
+  return coords;
+}
+
+/** Area (km²) of a closed ring, 0 when degenerate/NaN. */
+function ringKm2(ring) {
+  if (!isValidRing(ring)) return 0;
+  try { return area(polygon([ring])) * 1e-6; } catch { return 0; }
+}
+
+/** A valid GeoJSON ring: ≥4 closed positions, each [finiteNumber, finiteNumber]. */
+function isValidRing(ring) {
+  if (!Array.isArray(ring) || ring.length < 4) return false;
+  for (const c of ring) {
+    if (!Array.isArray(c) || c.length < 2) return false;
+    if (typeof c[0] !== "number" || !Number.isFinite(c[0])) return false;
+    if (typeof c[1] !== "number" || !Number.isFinite(c[1])) return false;
+  }
+  return true;
+}
+
+// Module-level refs populated in main before use (see keepHole).
+let preHoleCentroids = []; // centroids of holes present in the input file (kept verbatim)
+let otherLandmasses = [];  // { iso, feature } for every country landmass
+
+/**
+ * Is this hole a "real" one that must be preserved? Kept when it is a pre-existing
+ * region hole (centroid match, e.g. the Aral) or lies inside another country's
+ * landmass (enclave, e.g. Lesotho). Otherwise it is a polyclip phantom → dropped.
+ */
+function keepHole(hole, currentIso) {
+  const hKm2 = ringKm2(hole);
+  if (!(hKm2 >= RING_MIN_KM2)) return false; // degenerate / NaN
+  const cm = centroid(polygon([hole])).geometry.coordinates;
+  if (preHoleCentroids.some((p) => distance(cm, p) < 10)) return true;
+  if (otherLandmasses.some((l) => l.iso !== currentIso && booleanPointInPolygon(point(cm), l.feature))) return true;
+  return false;
+}
+
+/**
+ * Clean a GeoJSON polygon part (`[outerRing, ...holeRings]`): round coordinates, keep
+ * only real holes (see keepHole). Returns the cleaned part, or null if the outer ring
+ * is degenerate.
+ */
+function cleanPolygonPart(part, currentIso) {
+  if (!Array.isArray(part) || part.length === 0) return null;
+  const outer = part[0];
+  if (!isValidRing(outer)) return null;
+  const outerKm2 = ringKm2(outer);
+  if (!(outerKm2 >= RING_MIN_KM2)) return null; // covers NaN too
+  const holes = part.slice(1).filter((h) => isValidRing(h) && keepHole(h, currentIso));
+  roundCoords([outer, ...holes]);
+  return [outer, ...holes];
+}
+
+/** Split a geometry (Polygon|MultiPolygon) into cleaned polygon parts. */
+function partsOf(geom, iso) {
+  if (!geom) return [];
+  const raw = geom.type === "MultiPolygon" ? geom.coordinates : [geom.coordinates];
+  const out = [];
+  for (const p of raw) { const c = cleanPolygonPart(p, iso); if (c) out.push(c); }
+  return out;
+}
+
+/** Rebuild a geometry (NOT a feature) from cleaned parts. */
+function toGeometry(parts) {
+  // turf helpers return Features; take their .geometry ({type, coordinates}).
+  return parts.length > 1 ? multiPolygon(parts).geometry : polygon(parts[0]).geometry;
+}
+
+/** Bbox of a polygon part, as [west, south, east, north]. */
+function partBBox(part) { return bbox(polygon(part)); }
+
+function bBoxesOverlap(a, b) {
+  return !(a[2] < b[0] || a[0] > b[2] || a[3] < b[1] || a[1] > b[3]);
+}
+
+/** Sample points along a part's outer ring (start + every ~max(1, n/16)th vertex). */
+function samplePoints(part, maxN = 16) {
+  const ring = part[0];
+  const n = ring.length;
+  if (n <= maxN) return ring;
+  const step = Math.ceil(n / maxN);
+  const out = [ring[0]];
+  for (let i = step; i < n; i += step) out.push(ring[i]);
+  return out;
+}
+
+// Module-level list of pre-existing region holes (water guard), built in main.
+let holeRings = []; // { bbox, poly }
+
+/** Does a gap piece overlap any pre-existing region hole (water)? */
+function overlapsAnyHole(piece, pieceCm) {
+  const pBBox = partBBox(piece);
+  const candidates = holeRings.filter((h) => bBoxesOverlap(pBBox, h.bbox));
+  if (candidates.length === 0) return false;
+  const pts = [pieceCm, ...samplePoints(piece)];
+  return candidates.some((h) => pts.some((p) => booleanPointInPolygon(point(p), h.poly)));
+}
+
+// --------------------------------------------------------------------------- main
+
+const inFile = process.argv[2] ?? "public/regions.geo.json";
+const outFile = process.argv[3] ?? inFile;
+
+const geo = JSON.parse(fs.readFileSync(inFile, "utf8"));
+const countries = JSON.parse(fs.readFileSync(new URL("../public/countries.geo.json", import.meta.url), "utf8"));
+
+otherLandmasses = countries.features.map((f) => ({ iso: f.id, feature: f }));
+
+const features = geo.features;
+const regionsByParent = new Map();
+for (const f of features) {
+  const p = f.properties.parent;
+  if (!regionsByParent.has(p)) regionsByParent.set(p, []);
+  regionsByParent.get(p).push(f);
+}
+
+// Region inner-ring (hole) list for the water-hole guard; also the centroids of every
+// pre-existing hole (so the hole policy never destroys one, e.g. KZ-KZY's Aral hole).
+for (const f of features) {
+  const parts = f.geometry.type === "MultiPolygon" ? f.geometry.coordinates : [f.geometry.coordinates];
+  for (const p of parts) for (const h of p.slice(1)) {
+    if (Array.isArray(h) && h.length >= 4 && ringKm2(h) >= RING_MIN_KM2) {
+      const poly = polygon([h]);
+      holeRings.push({ bbox: bbox(poly), poly });
+      preHoleCentroids.push(centroid(poly).geometry.coordinates);
+    }
+  }
+}
+
+// Every part already in the input file (string key) — used to skip phantom pieces that
+// re-emiss an existing part byte-for-byte.
+const existingPartSet = new Set();
+for (const f of features) {
+  const parts = f.geometry.type === "MultiPolygon" ? f.geometry.coordinates : [f.geometry.coordinates];
+  for (const p of parts) existingPartSet.add(JSON.stringify(p));
+}
+
+// All regions with bboxes, for cross-border clipping.
+const allRegions = features.map((f) => ({
+  feature: f,
+  parent: f.properties.parent,
+  bbox: bbox(feature(f.geometry)),
+}));
+
+console.log(`input: ${inFile}`);
+console.log(`regions: ${features.length}  countries(with regions): ${regionsByParent.size}  region holes: ${holeRings.length}`);
+
+// Per-region running totals for the report.
+const added = new Map(); // feature -> { km2, parts }
+let totalSliverKm2 = 0, totalWaterKm2 = 0, totalFilledKm2 = 0, totalDupN = 0, totalClippedKm2 = 0, totalClipFailN = 0;
+
+// Pieces already appended this run (any country) — later pieces clip against these so
+// appended pieces never overlap each other either.
+const appendedSoFar = []; // { parent, bbox, part }
+
+for (const [iso, countryRegions] of regionsByParent) {
+  const landmass = countries.features.find((c) => c.id === iso);
+  if (!landmass) continue; // no authoritative landmass → nothing to fill against
+  const feats = countryRegions.map((r) => feature(r.geometry));
+
+  let unionGeom, diffGeom;
+  try {
+    // A country with a single sub-region has no union to compute: it IS the union.
+    unionGeom = countryRegions.length === 1
+      ? countryRegions[0].geometry
+      : union(featureCollection(feats)).geometry;
+    diffGeom = difference(featureCollection([landmass, { ...landmass, geometry: unionGeom }])).geometry;
+  } catch (e) {
+    console.log(`${iso}\tSKIPPED (${unionGeom ? "difference" : "union"} failed: ${e.message})`);
+    continue;
+  }
+  if (!diffGeom) continue;
+
+  const rawPieces = partsOf(diffGeom, iso); // cleaned gap fragments inside the landmass
+
+  // Pass 1: filter (sliver / water / duplicate) and clip against foreign territory.
+  const fillable = [];
+  let sliverKm2 = 0, waterKm2 = 0, dupN = 0, clippedKm2 = 0, sliverN = 0, waterN = 0, clipFailN = 0;
+  for (const piece of rawPieces) {
+    const piecePoly = polygon(piece);
+    const aKm2 = area(piecePoly) * 1e-6;
+    if (aKm2 < GAP_MIN_KM2) { sliverKm2 += aKm2; sliverN++; continue; }
+
+    const pieceCm = centroid(piecePoly).geometry.coordinates;
+    if (overlapsAnyHole(piece, pieceCm)) { waterKm2 += aKm2; waterN++; continue; }
+    if (existingPartSet.has(JSON.stringify(piece))) { dupN++; continue; }
+
+    // Clip: remove any part of this piece that ANY existing region (own or foreign
+    // country) or an already-appended piece claims. Clipped *sequentially* against
+    // each touching target (more robust than one big union-difference): a target that
+    // fails to clip is skipped and counted, not allowed to abort the whole piece.
+    // Clipping against the country's OWN regions as well as foreign ones is what kills
+    // the same-country overlaps: the gap is `landmass − union(own regions)`, and
+    // polychip's difference can leave thin slivers of an own region inside that "gap";
+    // re-differencing each fragment against the individual regions (simpler operands
+    // than the big union) strips them out. Both operands are built as proper Features
+    // carrying a real geometry — never `feature()` on an already-Feature object (that
+    // double-wraps the geometry into the Feature itself and breaks polychip).
+    const pBBox = partBBox(piece);
+    const targets = [
+      ...allRegions.filter((r) => bBoxesOverlap(r.bbox, pBBox)).map((r) => ({ feature: r.feature, bbox: r.bbox })),
+      ...appendedSoFar.filter((a) => bBoxesOverlap(a.bbox, pBBox)).map((a) => ({ feature: polygon(a.part), bbox: a.bbox })),
+    ];
+    if (targets.length === 0) {
+      fillable.push({ piece, aKm2, cm: pieceCm });
+      continue;
+    }
+    let cur = { type: "Feature", properties: {}, geometry: piecePoly.geometry };
+    let curBBox = pBBox;
+    let covered = false;
+    for (const t of targets) {
+      if (!bBoxesOverlap(curBBox, t.bbox)) continue;
+      let d;
+      try {
+        d = difference(featureCollection([cur, t.feature]));
+      } catch {
+        clipFailN++; // clip against this target failed: keep current (safe fallback)
+        continue;
+      }
+      if (!d || !d.geometry) { covered = true; break; } // piece fully covered by this target
+      cur = { type: "Feature", properties: {}, geometry: d.geometry };
+      curBBox = bbox(cur);
+    }
+    if (covered || !cur.geometry) { clippedKm2 += aKm2; continue; } // fully covered by foreign territory
+    const kept = [];
+    for (const cp of partsOf(cur.geometry, iso)) {
+      const cpArea = area(polygon(cp)) * 1e-6;
+      if (cpArea < GAP_MIN_KM2) { sliverKm2 += cpArea; sliverN++; continue; }
+      if (existingPartSet.has(JSON.stringify(cp))) { dupN++; continue; }
+      kept.push({ piece: cp, aKm2: cpArea, cm: centroid(polygon(cp)).geometry.coordinates });
+    }
+    clippedKm2 += Math.max(0, aKm2 - kept.reduce((s, k) => s + k.aKm2, 0));
+    fillable.push(...kept);
+  }
+
+  // Pass 2: nearest-region assignment and append.
+  const regionInfo = countryRegions.map((r) => {
+    const f = feature(r.geometry);
+    return { r, cm: centroid(f).geometry.coordinates, areaKm2: area(f) * 1e-6 };
+  });
+  let realKm2 = 0;
+  for (const { piece, aKm2, cm } of fillable) {
+    // Nearest region by centroid haversine distance; tie → larger area, then lexicographic id.
+    let best = null;
+    for (const info of regionInfo) {
+      const d = distance(cm, info.cm);
+      if (
+        best === null ||
+        d < best.d ||
+        (d === best.d && (info.areaKm2 > best.areaKm2 || (info.areaKm2 === best.areaKm2 && info.r.id < best.id)))
+      ) best = { info, d };
+    }
+    if (!best) continue; // no regions at all (should not happen)
+    realKm2 += aKm2;
+    const target = best.info.r;
+    const parts = (target.geometry.type === "MultiPolygon" ? target.geometry.coordinates : [target.geometry.coordinates]).concat([piece]);
+    target.geometry = toGeometry(parts);
+    const a = added.get(target) ?? { km2: 0, parts: 0 };
+    a.km2 += aKm2; a.parts += 1;
+    added.set(target, a);
+    appendedSoFar.push({ parent: iso, bbox: partBBox(piece), part: piece });
+  }
+
+  totalSliverKm2 += sliverKm2;
+  totalWaterKm2 += waterKm2;
+  totalFilledKm2 += realKm2;
+  totalDupN += dupN;
+  totalClippedKm2 += clippedKm2;
+  totalClipFailN += clipFailN;
+
+  // Verify the fill for countries we actually changed: recompute the residual gap.
+  let afterKm2 = null;
+  if (realKm2 > 0) {
+    try {
+      const u2 = countryRegions.length === 1
+        ? countryRegions[0].geometry
+        : union(featureCollection(countryRegions.map((r) => feature(r.geometry)))).geometry;
+      afterKm2 = area(difference(featureCollection([landmass, { ...landmass, geometry: u2 }]))) * 1e-6;
+    } catch { afterKm2 = null; }
+  }
+  console.log(
+    `${iso}\tregions=${countryRegions.length}\tgap-before=${(area(diffGeom) * 1e-6).toFixed(0)}` +
+    (afterKm2 === null ? "" : `\tgap-after=${afterKm2.toFixed(0)}`) +
+    `\tfilled=${realKm2.toFixed(0)}\tslivers=${sliverN}(${sliverKm2.toFixed(0)})\twater-skipped=${waterN}(${waterKm2.toFixed(0)})` +
+    `\tdup-skipped=${dupN}${clippedKm2 > 0 ? `\tforeign-overlap-removed=${clippedKm2.toFixed(0)}` : ""}${clipFailN > 0 ? `\tclip-failed=${clipFailN}` : ""}`,
+  );
+}
+
+const totalAddedKm2 = [...added.values()].reduce((s, a) => s + a.km2, 0);
+const totalAddedParts = [...added.values()].reduce((s, a) => s + a.parts, 0);
+console.log(`\nfilled ${totalAddedKm2.toFixed(0)} km² into ${added.size} regions (+${totalAddedParts} parts)`);
+console.log(
+  `totals: filled=${totalFilledKm2.toFixed(0)}  slivers-ignored=${totalSliverKm2.toFixed(0)}  water-kept-open=${totalWaterKm2.toFixed(0)}` +
+  `  duplicate-phantoms-skipped=${totalDupN}  foreign-overlap-removed=${totalClippedKm2.toFixed(0)}  clip-failed=${totalClipFailN}`,
+);
+
+fs.writeFileSync(outFile, JSON.stringify(geo));
+console.log(`wrote ${outFile}`);