History/scripts/fill-gaps.mjs

#!/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}`);