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