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