History/scripts/generate-regions.mjs

#!/usr/bin/env node
/**
 * scripts/generate-regions.mjs — dependency-free generator for `public/regions.geo.json`.
 *
 * Subdivides every country in `public/countries.geo.json` into sub-national
 * regions sized ~100–200 km on a side (adaptive: a regular lat/lon grid
 * clipped to the country, so vast/sparse countries — Siberia, Sahara,
 * Antarctica — end up with bigger regions while Europe gets smaller ones).
 * Each country gets a few to tens of regions (capped at 36); tiny countries
 * (and tiny antimeridian islet groups) stay a single region.
 *
 * Model: a country's lat/lon bbox is tiled with a cols × rows grid of
 * ~square-in-km cells. Each grid cell yields at most ONE region feature:
 * every sub-polygon of the country is clipped against the cell with
 * Sutherland–Hodgman (edges in order: left, bottom, right, top), clipped
 * pieces smaller than 300 km² are dropped, and the kept pieces become the
 * rings of that cell's region (a MultiPolygon when there is more than one).
 * Region ids are `<ISO_A3>-<nn>`, nn = 1..k in row-major order (rows
 * top→bottom by descending lat, left→right by ascending lon); names are
 * `<Country> R<nn>`. Coordinates are rounded to 4 decimal places (≈11 m).
 * Because a cell produces at most one region, the per-country count can
 * never exceed the (capped) cell count ≤ 36, even for archipelagos.
 *
 * Area: local planar approximation, x = lon·111.32·cos(midLat), y = lat·
 * 110.574 (km), shoelace per sub-polygon ring, summed in abs. Target region
 * area T = 12,000 km² (≈110 km side); target count n = clamp(round(area /
 * T), 1, 36). S = sqrt(bboxAreaKm2 / n); cols = max(1, round(wKm/S)),
 * rows = max(1, round(hKm/S)); the grid is then reduced if cols·rows > 36.
 *
 * Antimeridian handling (three cases):
 *  - Small antimeridian island groups (area ≤ 50,000 km² — Fiji, Kiribati,
 *    Tuvalu, …): one region = the whole country (`<ISO_A3>-01`, plain name).
 *  - Large antimeridian countries made of ordinary sub-polygons (Russia):
 *    the sub-polygons are shifted into longitude space [0, 360) (lons < 0 →
 *    +360) so the country sits in one contiguous band, tiled there, and
 *    output rings that straddle ±180° are split at the antimeridian into a
 *    MultiPolygon with each ring fully on one side — standard −180..180
 *    GeoJSON that renders without cross-map lines.
 *  - Full-360° "polar" wrap rings (Antarctica's main ring, whose only
 *    longitude jump sits at the pole and which crosses the prime meridian
 *    once): such a ring self-intersects in any single planar representation,
 *    so it is replaced by two simple "wedge" rings — the eastern and western
 *    hemispheric halves, each closed along the prime meridian, the pole
 *    (lat −90), and the antimeridian. The wedges tile −180..180 space
 *    without overlap and are tiled like any other ring (no shifting needed).
 *
 * Deterministic: same input → byte-identical output.
 * Run: `npm run generate:regions`
 */

import { readFileSync, writeFileSync } from 'node:fs';
import path from 'node:path';
import { fileURLToPath } from 'node:url';

const here = path.dirname(fileURLToPath(import.meta.url));
const root = path.resolve(here, '..');

// --- parameters ------------------------------------------------------------------
const T = 12_000; // target region area, km² (≈110 km side)
const MAX_REGIONS = 36; // per-country cap ("tens")
const SLIVER_KM2 = 300; // drop clipped pieces smaller than this
const BIG_KM2 = 50_000; // antimeridian countries at/above this are tiled, not 1 region
const KM_LON = 111.32; // km per degree of longitude (at the equator)
const KM_LAT = 110.574; // km per degree of latitude

// --- geo helpers (same planar approximation as scripts/regions.test.ts) -----------
const rad = (deg) => (deg * Math.PI) / 180;
const r4 = (v) => Math.round(v * 1e4) / 1e4;

function polygons(geom) {
  return geom.type === 'Polygon' ? [geom.coordinates] : geom.coordinates;
}

function allCoords(geom) {
  const out = [];
  for (const poly of polygons(geom)) for (const ring of poly) for (const c of ring) out.push(c);
  return out;
}

function meanLat(geom) {
  const cs = allCoords(geom);
  let s = 0;
  for (const c of cs) s += c[1];
  return cs.length ? s / cs.length : 0;
}

function bbox(rings) {
  let minLon = Infinity,
    maxLon = -Infinity,
    minLat = Infinity,
    maxLat = -Infinity;
  for (const ring of rings)
    for (const [lo, la] of ring) {
      if (lo < minLon) minLon = lo;
      if (lo > maxLon) maxLon = lo;
      if (la < minLat) minLat = la;
      if (la > maxLat) maxLat = la;
    }
  return { minLon, maxLon, minLat, maxLat };
}

/** |shoelace| of a (possibly open) ring in km², using one mid-latitude for x. */
function ringAreaKm2(ring, midLat) {
  const kx = KM_LON * Math.cos(rad(midLat));
  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);
}

/** Remove consecutive duplicate points and a trailing closing point equal to the first. */
function dedupeRing(ring) {
  const out = [];
  for (const p of ring) {
    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();
  return out;
}

/** Close an open ring by appending a copy of the first point. */
function closeRing(ring) {
  return ring.length ? [...ring, [...ring[0]]] : ring;
}

/**
 * The tiling/area subject rings of a feature, in standard −180..180 space:
 * one open ring per sub-polygon (outer ring only — holes are ignored by
 * design: Sutherland–Hodgman clips simple rings and any enclosed water is
 * at most a lake, which the area-sanity tolerance in the test absorbs).
 * A full-360° wrap ring is replaced by its two hemispheric wedges (see
 * makeWedges); if that is not possible the ring is returned as-is (the
 * caller may still shift it).
 */
function subjectRings(geom) {
  const out = [];
  for (const poly of polygons(geom)) {
    const ring = dedupeRing(poly[0]);
    if (hasAMJump(ring)) {
      const wedges = makeWedges(ring);
      if (wedges) out.push(...wedges);
      else out.push(ring);
    } else {
      out.push(ring);
    }
  }
  return out;
}

// --- Sutherland–Hodgman: clip an open subject ring against one half-plane --------
function clipEdge(poly, inside, intersect) {
  const out = [];
  for (let i = 0; i < poly.length; i++) {
    const a = poly[i];
    const b = poly[(i + 1) % poly.length];
    const aIn = inside(a);
    const bIn = inside(b);
    if (aIn && bIn) out.push(b);
    else if (aIn && !bIn) out.push(intersect(a, b));
    else if (!aIn && bIn) {
      out.push(intersect(a, b));
      out.push(b);
    }
  }
  return out;
}

/**
 * Clip an open subject ring against an axis-aligned (lon/lat) rectangle.
 * Edges are processed in order: left, bottom, right, top. Returns the clipped
 * open ring, or [] if empty.
 */
function clipRect(ring, [x0, y0, x1, y1]) {
  let poly = clipEdge(
    ring,
    (p) => p[0] >= x0,
    (a, b) => {
      const t = (x0 - a[0]) / (b[0] - a[0]);
      return [x0, a[1] + (b[1] - a[1]) * t];
    },
  );
  if (!poly.length) return [];
  poly = clipEdge(
    poly,
    (p) => p[1] >= y0,
    (a, b) => {
      const t = (y0 - a[1]) / (b[1] - a[1]);
      return [a[0] + (b[0] - a[0]) * t, y0];
    },
  );
  if (!poly.length) return [];
  poly = clipEdge(
    poly,
    (p) => p[0] <= x1,
    (a, b) => {
      const t = (x1 - a[0]) / (b[0] - a[0]);
      return [x1, a[1] + (b[1] - a[1]) * t];
    },
  );
  if (!poly.length) return [];
  poly = clipEdge(
    poly,
    (p) => p[1] <= y1,
    (a, b) => {
      const t = (y1 - a[1]) / (b[1] - a[1]);
      return [a[0] + (b[0] - a[0]) * t, y1];
    },
  );
  return poly;
}

// --- antimeridian handling ---------------------------------------------------------

/** True if the ring has an edge with |Δlon| > 180 (an antimeridian jump). */
function hasAMJump(ring) {
  const n = ring.length;
  for (let i = 0; i < n; i++)
    if (Math.abs(ring[i][0] - ring[(i + 1) % n][0]) > 180) return true;
  return false;
}

function crossesAntimeridian(geom) {
  let minLon = Infinity,
    maxLon = -Infinity;
  for (const [lo] of allCoords(geom)) {
    if (lo < minLon) minLon = lo;
    if (lo > maxLon) maxLon = lo;
  }
  return maxLon > 0 && minLon < 0 && maxLon - minLon > 180;
}

/**
 * Replace a full-360° wrap ring (one antimeridian jump, one prime-meridian
 * crossing) by two simple hemispheric "wedge" rings, each closed along the
 * prime meridian (x = 0), the pole (lat = −90), and the antimeridian
 * (x = ±180). Both wedges lie in standard −180..180 space, are disjoint, and
 * together cover exactly the original ring's enclosed region. Returns
 * [westRing, eastRing] (open rings), or null if the ring has an unexpected
 * shape (no wrap handling possible).
 */
function makeWedges(ring) {
  const n = ring.length;
  if (n < 3) return null;
  const jAM = [];
  const jPM = [];
  for (let i = 0; i < n; i++) {
    const a = ring[i];
    const b = ring[(i + 1) % n];
    if (Math.abs(a[0] - b[0]) > 180) jAM.push(i);
    else if (a[0] !== 0 && b[0] !== 0 && (a[0] < 0) !== (b[0] < 0)) jPM.push(i);
  }
  if (jAM.length !== 1 || jPM.length !== 1) return null;

  const iAM = jAM[0];
  const a = ring[iAM]; // the +180 end of the jump
  const b = ring[(iAM + 1) % n]; // the −180 end of the jump
  if (Math.abs(a[0] - 180) > 1e-6 || Math.abs(b[0] + 180) > 1e-6) return null;

  const iPM = jPM[0];
  const p = ring[iPM];
  const q = ring[(iPM + 1) % n];
  const t = (0 - p[0]) / (q[0] - p[0]);
  const P = [0, p[1] + (q[1] - p[1]) * t]; // the exact prime-meridian crossing

  // Chain from b (−180) forward to P: everything in between must have lon ≤ 0.
  const west = [];
  for (let i = (iAM + 1) % n; ; i = (i + 1) % n) {
    const pt = ring[i];
    if (pt[0] > 0) return null;
    west.push(pt);
    if (i === iPM) break;
  }
  west.push(P);
  // Chain from P forward to a (+180): everything in between must have lon ≥ 0.
  const east = [P];
  for (let i = (iPM + 1) % n; ; i = (i + 1) % n) {
    const pt = ring[i];
    if (pt[0] < 0) return null;
    east.push(pt);
    if (i === iAM) break;
  }

  const pole = (x) => [x, -90];
  // West wedge: P → pole along x=0 → pole along the pole → up the −180 side →
  // the western coast back to P. East wedge: the mirror image on x ≥ 0.
  const westWedge = dedupeRing([P, pole(0), pole(-180), b, ...west.slice(1)]);
  // The east chain q→…→a runs in *increasing* lon, so to trace the coast back
  // from a to q (keeping the wedge simple) we walk it in reverse.
  const eastWedge = dedupeRing([P, pole(0), pole(180), a, ...east.slice(1, -1).reverse()]);
  if (westWedge.length < 3 || eastWedge.length < 3) return null;
  return [westWedge, eastWedge];
}

/**
 * Split an open ring (in shifted [0,360) space) at the antimeridian (x = 180).
 *
 * Instead of walking the ring between crossings (which breaks on jagged
 * coastlines that cross the antimeridian many times), clip the ring against
 * the two half-planes x ≤ 180 and x ≥ 180 with Sutherland–Hodgman. Each piece
 * is a proper closed-on-the-split-line polygon, the two pieces share only the
 * x = 180 boundary, and the split is area-preserving (the boundary has zero
 * area). Both pieces are returned in shifted space; the caller remaps any
 * x > 180 piece back into −180..0.
 */
function splitAtAM(ring) {
  if (ring.length < 3) return [];
  const cross = (a, b) => {
    const t = (180 - a[0]) / (b[0] - a[0]);
    return [180, a[1] + (b[1] - a[1]) * t];
  };
  const left = clipEdge(ring, (p) => p[0] <= 180, cross); // x in [0, 180]
  // Right piece: x in [180, 360). Remap EVERY vertex by -360 (including the
  // shared x = 180 boundary vertex → -180) so the piece sits fully on the west
  // side of the antimeridian with no ±180 jump inside the ring.
  const right = clipEdge(ring, (p) => p[0] >= 180, cross).map(
    ([lo, la]) => [lo - 360, la],
  );
  const pieces = [];
  if (left.length >= 3) pieces.push(left);
  if (right.length >= 3) pieces.push(right);
  return pieces;
}

// --- feature processing ------------------------------------------------------------

function makeFeature(id, name, parent, rings) {
  // Each ring is an INDEPENDENT piece (a separate island or a separate clipped
  // cell fragment), never a hole. GeoJSON: a MultiPolygon is an array of
  // polygons, each polygon an array of rings — so wrap every ring in its own
  // (single-ring) polygon. ([rings] would instead make ring 0 the outer and
  // every other ring a hole, swallowing e.g. Fiji's other islands.)
  return {
    type: 'Feature',
    id,
    properties: { name, parent },
    geometry:
      rings.length === 1
        ? { type: 'Polygon', coordinates: [rings[0]] }
        : { type: 'MultiPolygon', coordinates: rings.map((r) => [r]) },
  };
}

/**
 * Whole-country single region (no tiling, no sliver drops): all sub-polygons
 * in one feature, coordinates rounded. Used for tiny countries (n = 1) and
 * small antimeridian-crossing island groups.
 */
function wholeCountryFeature(feature) {
  const id = feature.id;
  const name = feature.properties.name;
  const rings = polygons(feature.geometry).map((poly) =>
    closeRing(dedupeRing(poly[0]).map(([lo, la]) => [r4(lo), r4(la)])),
  );
  return [makeFeature(`${id}-01`, name, id, rings)];
}

/**
 * Tile the country's bbox with a cols × rows grid; each non-empty cell yields
 * one region feature (its kept clipped pieces as a Polygon/MultiPolygon).
 * `shifted` = true for large non-wrapping antimeridian countries (Russia):
 * work in [0, 360) longitude space, then split output rings at ±180°.
 */
function tileFeature(feature, rings, shifted) {
  const id = feature.id;
  const name = feature.properties.name;
  const midLat = meanLat(feature.geometry);
  const area = rings.reduce((s, r) => s + ringAreaKm2(r, midLat), 0);

  const n = Math.min(MAX_REGIONS, Math.max(1, Math.round(area / T)));
  if (n === 1) return wholeCountryFeature(feature);

  const shiftLon = (lo) => (shifted && lo < 0 ? lo + 360 : lo);
  const work = rings.map((r) => r.map(([lo, la]) => [shiftLon(lo), la]));
  const bb = bbox(work);

  const kx = KM_LON * Math.cos(rad(midLat));
  const wKm = (bb.maxLon - bb.minLon) * kx;
  const hKm = (bb.maxLat - bb.minLat) * KM_LAT;
  const S = Math.sqrt((wKm * hKm) / n);
  let cols = Math.max(1, Math.round(wKm / S));
  let rows = Math.max(1, Math.round(hKm / S));
  while (cols * rows > MAX_REGIONS) {
    if (cols >= rows) cols = Math.max(1, Math.floor(MAX_REGIONS / rows));
    else rows = Math.max(1, Math.floor(MAX_REGIONS / cols));
  }

  const dLon = (bb.maxLon - bb.minLon) / cols;
  const dLat = (bb.maxLat - bb.minLat) / rows;

  const out = [];
  let nn = 0;
  for (let r = 0; r < rows; r++) {
    const yTop = bb.maxLat - r * dLat;
    const yBot = bb.maxLat - (r + 1) * dLat;
    for (let c = 0; c < cols; c++) {
      const xLeft = bb.minLon + c * dLon;
      const xRight = bb.minLon + (c + 1) * dLon;
      // One region per cell: clip every sub-ring, keep the pieces that clear
      // the sliver threshold, and merge them into a single feature.
      const cellPieces = [];
      for (const ring of work) {
        const clipped = clipRect(ring, [xLeft, yBot, xRight, yTop]);
        if (clipped.length < 3) continue;
        if (ringAreaKm2(clipped, midLat) < SLIVER_KM2) continue;
        // splitAtAM returns pieces already in standard −180..180 space (the
        // x > 180 half remapped), so no further per-vertex remap is needed.
        const pieces = shifted
          ? splitAtAM(dedupeRing(clipped)).map((pc) => dedupeRing(pc))
          : [dedupeRing(clipped)];
        for (const pc of pieces)
          if (pc.length >= 3) cellPieces.push(pc.map(([lo, la]) => [r4(lo), r4(la)]));
      }
      if (!cellPieces.length) continue;
      const closed = cellPieces
        .map(dedupeRing)
        .filter((pc) => pc.length >= 3)
        .map(closeRing);
      if (!closed.length) continue;
      nn += 1;
      const nnStr = String(nn).padStart(2, '0');
      out.push(makeFeature(`${id}-${nnStr}`, `${name} R${nnStr}`, id, closed));
    }
  }
  return out;
}

function processFeature(feature) {
  const geom = feature.geometry;
  const am = crossesAntimeridian(geom);
  if (am && featureAreaKm2(feature) <= BIG_KM2) {
    // Small antimeridian island group (Fiji, Kiribati, Tuvalu, …): one region.
    return wholeCountryFeature(feature);
  }
  // Wrapping rings (Antarctica) are already in standard space as wedges — no
  // shifting; non-wrapping antimeridian countries (Russia) are shifted.
  const rings = subjectRings(geom);
  const hasWrap = polygons(geom).some((p) => hasAMJump(dedupeRing(p[0])));
  const shifted = am && !hasWrap;
  return tileFeature(feature, rings, shifted);
}

/** Total planar area of a feature's tiling rings (same convention as tileFeature). */
function featureAreaKm2(feature) {
  const midLat = meanLat(feature.geometry);
  return subjectRings(feature.geometry).reduce((s, r) => s + ringAreaKm2(r, midLat), 0);
}

// --- main ---------------------------------------------------------------------------

function main() {
  const src = readFileSync(path.join(root, 'public', 'countries.geo.json'), 'utf8');
  const countries = JSON.parse(src);

  const features = [];
  const amSingletons = [];
  const amTiled = [];
  for (const f of countries.features) {
    if (crossesAntimeridian(f.geometry)) {
      (featureAreaKm2(f) <= BIG_KM2 ? amSingletons : amTiled).push(f.id);
    }
    features.push(...processFeature(f));
  }
  // Sorted by (parent, nn) — the processing order already guarantees this, but
  // sort explicitly so the contract holds regardless of traversal changes.
  features.sort((a, b) => {
    if (a.properties.parent !== b.properties.parent)
      return a.properties.parent < b.properties.parent ? -1 : 1;
    return a.id.localeCompare(b.id, 'en', { numeric: true });
  });

  const out = { type: 'FeatureCollection', features };
  const json = JSON.stringify(out);
  writeFileSync(path.join(root, 'public', 'regions.geo.json'), json + '\n');

  // --- summary ----------------------------------------------------------------------
  const counts = new Map();
  for (const f of features)
    counts.set(f.properties.parent, (counts.get(f.properties.parent) ?? 0) + 1);
  const sortedCounts = [...counts.values()].sort((a, b) => a - b);
  const median =
    sortedCounts.length % 2
      ? sortedCounts[Math.floor(sortedCounts.length / 2)]
      : (sortedCounts[sortedCounts.length / 2 - 1] + sortedCounts[sortedCounts.length / 2]) / 2;

  // Max region area (planar approximation, same mid-latitude convention as the
  // tiling: the parent country's mean latitude).
  const midLatByParent = new Map(countries.features.map((f) => [f.id, meanLat(f.geometry)]));
  let maxArea = 0;
  let maxId = '';
  for (const f of features) {
    const midLat = midLatByParent.get(f.properties.parent);
    let a = 0;
    for (const poly of polygons(f.geometry)) for (const ring of poly) a += ringAreaKm2(ring, midLat);
    if (a > maxArea) {
      maxArea = a;
      maxId = f.id;
    }
  }

  console.log(`regions: ${features.length} total (${countries.features.length} countries)`);
  console.log(`file:    public/regions.geo.json — ${(Buffer.byteLength(json + '\n') / 1024).toFixed(1)} KB`);
  console.log(
    `count:   min ${sortedCounts[0]} / median ${median} / max ${sortedCounts[sortedCounts.length - 1]} per country`,
  );
  console.log(`max region: ${maxId} ≈ ${Math.round(maxArea).toLocaleString('en-US')} km²`);
  if (amSingletons.length) console.log(`AM → 1 region: ${amSingletons.join(', ')}`);
  if (amTiled.length) console.log(`AM → tiled (shifted/wedges): ${amTiled.join(', ')}`);
}

main();