324 lines
15 KiB
JavaScript
324 lines
15 KiB
JavaScript
#!/usr/bin/env node
|
||
/**
|
||
* scripts/cleanup-regions.mjs — removes the small holes left by coarsening.
|
||
*
|
||
* Run: npm run cleanup:regions (in place on public/regions.geo.json)
|
||
* node scripts/cleanup-regions.mjs [out] [src] [--bridge]
|
||
*
|
||
* The v0.14 coarsening produced a clean 1,796-region dataset, but small
|
||
* unfilled patches ("holes") remain inside regions — artifacts of the boolean
|
||
* dissolves and Natural Earth's own digitization. This fixes that in place,
|
||
* WITHOUT touching the dataset's contract (the region count, the `id`s, the
|
||
* `properties = { name, parent }`, the ~100 km floor and the island exemption
|
||
* are all preserved — only inner rings are dropped).
|
||
*
|
||
* ENCLOSED HOLES — every hole (inner ring) whose area is < HOLE_AREA_MAX_KM2
|
||
* is dropped, so the region fills in to be solid. Any larger hole (a genuine
|
||
* water body) is kept. NOTE: the big real lakes (Great Lakes, Caspian, Aral,
|
||
* Lake Superior, …) were already dissolved by the v0.14 coarsening, so they
|
||
* are not present as holes here; the only hole above the threshold in this
|
||
* dataset is a single medium one (~5,858 km2, kept). A 1,000 km2 threshold
|
||
* is comfortably above the 2,252 tiny artifact holes (all < 500 km2) and
|
||
* comfortably below any large real lake, so it fills exactly the artifacts.
|
||
*
|
||
* BRIDGING (opt-in, `--bridge`) — a MultiPolygon whose separate parts sit
|
||
* only a short distance apart (a thin gap across a strait, a nearby island)
|
||
* leaves a visible crack of background. This joins each such close pair with
|
||
* a small connector polygon and a boolean union, closing the crack while
|
||
* genuinely far parts (islands, separated territories) stay separate. It is
|
||
* OFF by default: the union itself can mint new (medium) holes, which is
|
||
* counter to the primary "fill holes" goal, so the default run is the
|
||
* clean, deterministic hole-fill only.
|
||
*
|
||
* Deterministic + idempotent: fixed threshold, deterministic planar geometry,
|
||
* 4-dp rounding. Re-running on an already-clean file rewrites the same file
|
||
* (no small holes remain to fill; the kept holes are above the threshold).
|
||
*/
|
||
|
||
import { readFileSync, writeFileSync } from 'node:fs';
|
||
import path from 'node:path';
|
||
import { fileURLToPath } from 'node:url';
|
||
import { union } from '@turf/union';
|
||
import { featureCollection } from '@turf/helpers';
|
||
|
||
const here = path.dirname(fileURLToPath(import.meta.url));
|
||
const root = path.resolve(here, '..');
|
||
const DEFAULT_FILE = path.join(root, 'public', 'regions.geo.json');
|
||
// positional args (out, src) plus named flags; keep the two orthogonal
|
||
const args = process.argv.slice(2).filter((a) => !a.startsWith('-'));
|
||
const flags = process.argv.slice(2).filter((a) => a.startsWith('-'));
|
||
const BRIDGE = flags.includes('--bridge');
|
||
const OUT = args[0] ? path.resolve(args[0]) : DEFAULT_FILE;
|
||
const SRC = args[1] ? path.resolve(args[1]) : DEFAULT_FILE;
|
||
|
||
// --- parameters ---------------------------------------------------------------
|
||
const HOLE_AREA_MAX_KM2 = 1000; // fill enclosed holes smaller than this (keep real lakes)
|
||
const BRIDGE_MAX_KM = 30; // bridge part-gaps smaller than this …
|
||
const BRIDGE_MIN_KM = 0.3; // … but only real gaps (a pair in (BRIDGE_MIN, BRIDGE_MAX))
|
||
const BRIDGE_EXT_KM = 2; // extend each bridge end into its part (ensures union overlap)
|
||
const BRIDGE_HW_MIN_KM = 2; // minimum half-width of a connector
|
||
const BRIDGE_HW_MAX_KM = 20; // cap so a connector never becomes a big blob
|
||
const R4 = 4; // coordinate precision (decimal places)
|
||
const SLIVER_EPS_KM2 = 1e-6; // drop only truly-degenerate (≈zero-area) parts
|
||
|
||
// --- area model (identical to coarsen-regions.mjs / regions.test.ts) -----------
|
||
const KM_LON = 111.32;
|
||
const KM_LAT = 110.574;
|
||
const rad = (deg) => (deg * Math.PI) / 180;
|
||
function ringAreaKm2(ring) {
|
||
let lat = 0;
|
||
for (const [, la] of ring) lat += la;
|
||
const kx = KM_LON * Math.cos(rad(ring.length ? lat / ring.length : 0));
|
||
let a = 0;
|
||
const n = ring.length;
|
||
for (let i = 0; i < n; i++) {
|
||
const [lon1, lat1] = ring[i];
|
||
const [lon2, lat2] = ring[(i + 1) % n];
|
||
a += lon1 * kx * (lat2 * KM_LAT) - lon2 * kx * (lat1 * KM_LAT);
|
||
}
|
||
return Math.abs(a / 2);
|
||
}
|
||
function geomAreaKm2(geom) {
|
||
if (!geom) return 0;
|
||
const polys = geom.type === 'Polygon' ? [geom.coordinates] : geom.coordinates;
|
||
let total = 0;
|
||
for (const poly of polys) total += ringAreaKm2(poly[0]);
|
||
return total;
|
||
}
|
||
|
||
// --- rounding / normalization (same conventions as coarsen-regions.mjs) --------
|
||
const rN = (v, N) => Math.round(v * 10 ** N) / 10 ** N;
|
||
function roundRing(ring, N = R4) {
|
||
const out = [];
|
||
for (const [x, y] of ring) {
|
||
const p = [rN(x, N), rN(y, N)];
|
||
const last = out[out.length - 1];
|
||
if (!last || last[0] !== p[0] || last[1] !== p[1]) out.push(p);
|
||
}
|
||
if (out.length > 1 && out[0][0] === out[out.length - 1][0] && out[0][1] === out[out.length - 1][1]) out.pop();
|
||
if (out.length < 3) return out;
|
||
return [...out, [...out[0]]];
|
||
}
|
||
function roundAll(geom) {
|
||
if (!geom || (geom.type !== 'Polygon' && geom.type !== 'MultiPolygon')) return null;
|
||
const polys = (geom.type === 'Polygon' ? [geom.coordinates] : geom.coordinates)
|
||
.map((poly) => poly.map((ring) => roundRing(ring)))
|
||
.filter((poly) => poly.length >= 1 && poly.every((r) => r.length >= 3));
|
||
if (!polys.length) return null;
|
||
if (polys.length === 1 && geom.type === 'Polygon') return { type: 'Polygon', coordinates: polys[0] };
|
||
return { type: 'MultiPolygon', coordinates: polys };
|
||
}
|
||
function normalizeGeom(geom) {
|
||
const rounded = roundAll(geom);
|
||
if (!rounded) return null;
|
||
const all = rounded.type === 'Polygon' ? [rounded.coordinates] : rounded.coordinates;
|
||
const kept = all.filter((poly) => ringAreaKm2(poly[0]) >= SLIVER_EPS_KM2);
|
||
const polys = kept.length ? kept : all;
|
||
if (polys.length === 1 && rounded.type === 'Polygon') return { type: 'Polygon', coordinates: polys[0] };
|
||
return { type: 'MultiPolygon', coordinates: polys };
|
||
}
|
||
|
||
// --- closest-point measurement (decimated planar point↔segment, ~km) ----------
|
||
function toXY(ring) {
|
||
let lat = 0;
|
||
for (const [, la] of ring) lat += la;
|
||
const kx = KM_LON * Math.cos(rad(ring.length ? lat / ring.length : 0));
|
||
return ring.map(([lon, lat]) => [lon * kx, lat * KM_LAT]);
|
||
}
|
||
function dec(pts, maxN = 200) {
|
||
if (pts.length <= maxN) return pts;
|
||
const step = Math.ceil(pts.length / maxN);
|
||
const out = [];
|
||
for (let i = 0; i < pts.length; i += step) out.push(pts[i]);
|
||
return out;
|
||
}
|
||
/**
|
||
* Closest point-pair between two planar rings: returns { d, p, q } with p on A,
|
||
* q on B and d = |p-q| in km (both on the decimated rings — a good-enough
|
||
* approximation for bridging; a tiny offset is corrected by the extension).
|
||
*/
|
||
function ringClosest(AX, BX) {
|
||
const A = dec(AX), B = dec(BX);
|
||
let best = Infinity, bestP = null, bestQ = null;
|
||
for (const p of A) {
|
||
for (let i = 0; i < B.length - 1; i++) {
|
||
const a = B[i], b = B[i + 1];
|
||
const dx = b[0] - a[0], dy = b[1] - a[1];
|
||
const l2 = dx * dx + dy * dy;
|
||
let t = l2 === 0 ? 0 : ((p[0] - a[0]) * dx + (p[1] - a[1]) * dy) / l2;
|
||
t = Math.max(0, Math.min(1, t));
|
||
const c = [a[0] + t * dx, a[1] + t * dy];
|
||
const dd = Math.hypot(p[0] - c[0], p[1] - c[1]);
|
||
if (dd < best) { best = dd; bestP = p; bestQ = c; }
|
||
}
|
||
}
|
||
for (const p of B) {
|
||
for (let i = 0; i < A.length - 1; i++) {
|
||
const a = A[i], b = A[i + 1];
|
||
const dx = b[0] - a[0], dy = b[1] - a[1];
|
||
const l2 = dx * dx + dy * dy;
|
||
let t = l2 === 0 ? 0 : ((p[0] - a[0]) * dx + (p[1] - a[1]) * dy) / l2;
|
||
t = Math.max(0, Math.min(1, t));
|
||
const c = [a[0] + t * dx, a[1] + t * dy];
|
||
const dd = Math.hypot(p[0] - c[0], p[1] - c[1]);
|
||
if (dd < best) { best = dd; bestP = c; bestQ = p; }
|
||
}
|
||
}
|
||
return { d: best, p: bestP, q: bestQ };
|
||
}
|
||
|
||
/**
|
||
* Build a thin connector (a quad in lon/lat) between two closest points p, q
|
||
* (both in planar km). The connector is extended a little into each part so the
|
||
* boolean union joins them, and its half-width is wide enough to span the local
|
||
* gap (capped so it never becomes a big blob).
|
||
*/
|
||
function makeBridgeRing(p, q) {
|
||
const gap = Math.hypot(q[0] - p[0], q[1] - p[1]);
|
||
if (gap < 1e-3) return null; // coincident → nothing to bridge
|
||
const midLat = (p[1] / KM_LAT + q[1] / KM_LAT) / 2;
|
||
const kx = KM_LON * Math.cos(rad(midLat));
|
||
const ux = (q[0] - p[0]) / gap, uy = (q[1] - p[1]) / gap; // unit p→q
|
||
const nx = -uy, ny = ux; // perpendicular
|
||
const hw = Math.min(BRIDGE_HW_MAX_KM, Math.max(BRIDGE_HW_MIN_KM, gap / 2));
|
||
const ext = BRIDGE_EXT_KM;
|
||
const pExt = [p[0] - ux * ext, p[1] - uy * ext];
|
||
const qExt = [q[0] + ux * ext, q[1] + uy * ext];
|
||
const toLL = ([x, y]) => [rN(x / kx, R4), rN(y / KM_LAT, R4)];
|
||
const c1 = toLL([pExt[0] + nx * hw, pExt[1] + ny * hw]);
|
||
const c2 = toLL([qExt[0] + nx * hw, qExt[1] + ny * hw]);
|
||
const c3 = toLL([qExt[0] - nx * hw, qExt[1] - ny * hw]);
|
||
const c4 = toLL([pExt[0] - nx * hw, pExt[1] - ny * hw]);
|
||
return [c1, c2, c3, c4, c1];
|
||
}
|
||
|
||
/** Bridge all close pairs of a MultiPolygon and union everything into one shape. */
|
||
function bridgeParts(geom) {
|
||
const parts = geom.coordinates;
|
||
const partsXY = parts.map((p) => toXY(p[0]));
|
||
const bridges = [];
|
||
let bridgedPairs = 0;
|
||
for (let i = 0; i < partsXY.length; i++) {
|
||
for (let j = i + 1; j < partsXY.length; j++) {
|
||
const { d, p, q } = ringClosest(partsXY[i], partsXY[j]);
|
||
if (d > BRIDGE_MIN_KM && d < BRIDGE_MAX_KM) {
|
||
const ring = makeBridgeRing(p, q);
|
||
if (ring) { bridges.push(ring); bridgedPairs++; }
|
||
}
|
||
}
|
||
}
|
||
if (!bridges.length) return { geom, bridgedPairs: 0, changed: false };
|
||
// union: original parts (as Polygons) + connector rings (as Polygons)
|
||
const feats = parts.map((poly) => ({ type: 'Feature', geometry: { type: 'Polygon', coordinates: poly } }));
|
||
for (const ring of bridges) feats.push({ type: 'Feature', geometry: { type: 'Polygon', coordinates: [ring] } });
|
||
let u;
|
||
try {
|
||
u = union(featureCollection(feats));
|
||
} catch (e) {
|
||
console.error(` [warn] bridge union failed: ${e.message}`);
|
||
return { geom, bridgedPairs: 0, changed: false };
|
||
}
|
||
const norm = u ? normalizeGeom(u.geometry) : null;
|
||
// Guard: a good bridge only adds a little area. Reject the union if it lost
|
||
// area (< 0.8×, a bad union that ate parts) or over-filled (> 1.3×, e.g. many
|
||
// overlapping bridges in a dense archipelago) — keep the original instead.
|
||
if (norm) {
|
||
const na = geomAreaKm2(norm);
|
||
const oa = geomAreaKm2(geom);
|
||
if (oa > 0 && (na < 0.8 * oa || na > 1.3 * oa)) {
|
||
return { geom, bridgedPairs: 0, changed: false };
|
||
}
|
||
}
|
||
return { geom: norm ?? geom, bridgedPairs, changed: true, bridgedParts: u ? (u.geometry.type === 'Polygon' ? 1 : u.geometry.coordinates.length) : parts.length };
|
||
}
|
||
|
||
// --- main ----------------------------------------------------------------------
|
||
function main() {
|
||
const t0 = Date.now();
|
||
const fc = JSON.parse(readFileSync(SRC, 'utf8'));
|
||
const inKB = (Buffer.byteLength(readFileSync(SRC)) / 1024).toFixed(0);
|
||
console.log(`input: ${fc.features.length} regions (${inKB} KB) from ${SRC}`);
|
||
console.log(`mode: holes <${HOLE_AREA_MAX_KM2} km2 filled, larger kept${BRIDGE ? `; + bridging part-gaps in (${BRIDGE_MIN_KM}, ${BRIDGE_MAX_KM}) km` : ''}`);
|
||
|
||
// diagnostics
|
||
let holesFilled = 0, holesFilledArea = 0, holesKept = 0, holesKeptArea = 0, featuresTouchedHole = 0;
|
||
let bridgedFeatures = 0, bridgedPairs = 0, partsBefore = 0, partsAfter = 0, areaBefore = 0, areaAfter = 0;
|
||
const bridgeSamples = [];
|
||
const keptHoleSamples = [];
|
||
|
||
const outFeatures = [];
|
||
for (const f of fc.features) {
|
||
let geom = f.geometry;
|
||
|
||
// ---- Opt-in: bridge close parts (additive connector + union) ----
|
||
if (BRIDGE && geom.type === 'MultiPolygon' && geom.coordinates.length >= 2) {
|
||
const before = { parts: geom.coordinates.length, area: geomAreaKm2(geom) };
|
||
const res = bridgeParts(geom);
|
||
if (res.changed) {
|
||
geom = res.geom;
|
||
bridgedFeatures++;
|
||
bridgedPairs += res.bridgedPairs;
|
||
partsBefore += before.parts;
|
||
partsAfter += geom.type === 'Polygon' ? 1 : geom.coordinates.length;
|
||
areaBefore += before.area;
|
||
areaAfter += geomAreaKm2(geom);
|
||
if (bridgeSamples.length < 20)
|
||
bridgeSamples.push(`${f.id} ${f.properties.name} (${f.properties.parent}): ${before.parts}p→${geom.type === 'Polygon' ? 1 : geom.coordinates.length}p, ${res.bridgedPairs} bridge(s), ${Math.round(before.area)}→${Math.round(geomAreaKm2(geom))} km2`);
|
||
}
|
||
}
|
||
|
||
// ---- Fill small enclosed holes: drop inner rings < threshold, keep larger ----
|
||
const isMulti = geom.type === 'MultiPolygon';
|
||
const polys = isMulti ? geom.coordinates : [geom.coordinates];
|
||
let touchedHole = false;
|
||
const newPolys = polys.map((poly) => {
|
||
const kept = [poly[0]]; // outer ring always kept
|
||
for (let i = 1; i < poly.length; i++) {
|
||
const area = ringAreaKm2(poly[i]);
|
||
if (area < 1e-6) continue; // degenerate — drop
|
||
if (area < HOLE_AREA_MAX_KM2) {
|
||
holesFilled++; holesFilledArea += area; touchedHole = true; // fill (drop)
|
||
} else {
|
||
holesKept++; holesKeptArea += area;
|
||
kept.push(poly[i]); // a genuine water body — keep the hole
|
||
if (keptHoleSamples.length < 20) keptHoleSamples.push(`${f.id}(${f.properties.parent}) ${Math.round(area)} km2`);
|
||
}
|
||
}
|
||
return kept;
|
||
});
|
||
if (touchedHole) featuresTouchedHole++;
|
||
|
||
const rebuilt = isMulti
|
||
? { type: 'MultiPolygon', coordinates: newPolys }
|
||
: { type: 'Polygon', coordinates: newPolys[0] };
|
||
// normalize (round 4dp, drop degenerate). NEVER drop a feature: if nothing
|
||
// usable survives (a degenerate 3-pt 0-area micro-dot), keep the
|
||
// hole-filled (rebuilt) geometry, else the untouched original.
|
||
const norm = normalizeGeom(rebuilt);
|
||
outFeatures.push({ type: 'Feature', id: f.id, properties: { name: f.properties.name, parent: f.properties.parent }, geometry: norm ?? rebuilt ?? geom });
|
||
}
|
||
|
||
const out = { type: 'FeatureCollection', features: outFeatures };
|
||
const json = JSON.stringify(out);
|
||
const outKB = (Buffer.byteLength(json) / 1024).toFixed(0);
|
||
writeFileSync(OUT, json + '\n');
|
||
|
||
console.log(`\n--- RESULTS ---`);
|
||
console.log(`regions: ${fc.features.length} → ${out.features.length} (should be unchanged)`);
|
||
console.log(`file: ${inKB} KB → ${outKB} KB (${(outKB / 1024).toFixed(2)} MB)`);
|
||
console.log(`\nENCLOSED HOLES:`);
|
||
console.log(` filled ${holesFilled} holes (${Math.round(holesFilledArea).toLocaleString('en-US')} km2) in ${featuresTouchedHole} features`);
|
||
console.log(` kept ${holesKept} hole(s) (${Math.round(holesKeptArea).toLocaleString('en-US')} km2) as real water`);
|
||
if (keptHoleSamples.length) console.log(` kept-hole list: ${keptHoleSamples.join(', ')}`);
|
||
if (BRIDGE) {
|
||
console.log(`\nBRIDGING (part-gaps):`);
|
||
console.log(` bridged ${bridgedFeatures} features (${bridgedPairs} close pairs): parts ${partsBefore} → ${partsAfter}, area ${Math.round(areaBefore).toLocaleString('en-US')} → ${Math.round(areaAfter).toLocaleString('en-US')} km2`);
|
||
for (const s of bridgeSamples) console.log(` ${s}`);
|
||
} else {
|
||
console.log(`\nbridging: skipped (pass --bridge to enable)`);
|
||
}
|
||
console.log(`\nwrote ${OUT}; elapsed ${((Date.now() - t0) / 1000).toFixed(1)}s`);
|
||
}
|
||
|
||
main();
|