Coarsen sub-national regions to ~100km floor (v0.14.0)

Part of #16

scripts/regions.test.ts

@@ -1,12 +1,16 @@
 /**
  * scripts/regions.test.ts — validation for the committed `public/regions.geo.json`.
  *
- * The dataset is now REAL named sub-national administrative regions (Natural
- * Earth admin-1: states, provinces, governorates, oblasts, …) rather than the
- * old geometric grid — one feature per named region, each carrying
- * `properties = { name, parent }` where `parent` is the country's ISO_A3 (the
- * RegionMap key the map resolves ownership through). Tiny or uncovered
- * countries collapse to a single whole-country region (`<ISO_A3>-01`).
+ * The dataset is REAL named sub-national administrative regions (Natural
+ * Earth admin-1: states, provinces, governorates, oblasts, …) that the
+ * coarsening step (`scripts/coarsen-regions.mjs`, v0.14) has merged down to a
+ * ~100×100 km (10,000 km²) floor: every small non-island region was dissolved
+ * into its nearest same-country neighbor, so the only regions left below the
+ * floor are islands (no same-country sibling nearby). One feature per named
+ * region, each carrying `properties = { name, parent }` where `parent` is the
+ * country's ISO_A3 (the RegionMap key the map resolves ownership through).
+ * Tiny or uncovered countries collapse to a single whole-country region
+ * (`<ISO_A3>-01`).
  *
  * This test is self-contained (no import of the generator): it re-implements
  * the same planar area model the generator/analysis use, so the area-based
@@ -109,6 +113,99 @@ const countByParent = new Map<string, number>();
 for (const f of regions.features)
   countByParent.set(f.properties.parent, (countByParent.get(f.properties.parent) ?? 0) + 1);
 
+// --- island/close-neighbor helpers (mirror scripts/coarsen-regions.mjs) ------
+
+interface LonInt {
+  a: number;
+  b: number;
+  full: boolean;
+}
+interface RegionBBox {
+  latMin: number;
+  latMax: number;
+  lon: LonInt;
+}
+
+const COARSE_FLOOR_KM2 = 10_000; // the coarsener's ~100×100 km floor
+const CLOSE_DEG = 0.2; // "spatially close" bbox gap (both lat and lon ≤ this)
+
+/**
+ * Unwrap a ring's longitudes to a continuous line (np.unwrap shortest-diff
+ * scheme, anchored at the first vertex's [0,360) value). `full` when the
+ * unwrapped span is ≥ 300° (a full-360° antimeridian band, e.g. Antarctica).
+ * Must stay in sync with the coarsening script.
+ */
+function unwrapRing(ring: Ring) {
+  const n = ring.length;
+  const norm = (x: number) => ((x % 360) + 360) % 360;
+  const u: number[] = new Array(n);
+  u[0] = norm(ring[0][0]);
+  for (let i = 1; i < n; i++) {
+    // the previous UNWRAPPED value may lie outside [0,360); normalize it
+    // before measuring the shortest signed difference, or a bogus ~180° jump
+    // is inserted and the ring's span balloons.
+    const a = norm(u[i - 1]);
+    const b = norm(ring[i][0]);
+    const d = (((b - a) % 360) + 540) % 360 - 180;
+    u[i] = u[i - 1] + d;
+  }
+  let min = Infinity;
+  let max = -Infinity;
+  for (const v of u) {
+    if (v < min) min = v;
+    if (v > max) max = v;
+  }
+  return { lons: u, full: max - min >= 300 };
+}
+
+/**
+ * Cyclic 1-D gap (degrees) between two tight longitude intervals in their own
+ * unwrapped spaces; a `full` interval is close to everything.
+ */
+function lonGap(intA: LonInt, intB: LonInt): number {
+  if (intA.full || intB.full) return 0;
+  let best = Infinity;
+  for (const k of [-1, 0, 1]) {
+    const a2 = intB.a + 360 * k;
+    const b2 = intB.b + 360 * k;
+    const gap = intA.b < a2 ? a2 - intA.b : b2 < intA.a ? intA.a - b2 : 0;
+    if (gap < best) best = gap;
+  }
+  return best;
+}
+
+// Per-region bbox in antimeridian-safe unwrapped longitude space.
+const regionBBox = new Map<string, RegionBBox>();
+for (const f of regions.features) {
+  const polys = f.geometry.type === 'Polygon' ? [f.geometry.coordinates] : f.geometry.coordinates;
+  let latMin = Infinity;
+  let latMax = -Infinity;
+  let a = Infinity;
+  let b = -Infinity;
+  let full = false;
+  for (const poly of polys) {
+    const outer = poly[0];
+    for (const [, la] of outer) {
+      if (la < latMin) latMin = la;
+      if (la > latMax) latMax = la;
+    }
+    const { lons, full: fl } = unwrapRing(outer);
+    if (fl) full = true;
+    for (const v of lons) {
+      if (v < a) a = v;
+      if (v > b) b = v;
+    }
+  }
+  regionBBox.set(f.id, { latMin, latMax, lon: { a, b, full } });
+}
+
+/** Bounding-box closeness (same definition the coarsening script uses). */
+function bboxesClose(x: RegionBBox, y: RegionBBox): boolean {
+  const latGap = x.latMax < y.latMin ? y.latMin - x.latMax : y.latMax < x.latMin ? x.latMin - y.latMax : 0;
+  if (latGap > CLOSE_DEG) return false;
+  return lonGap(x.lon, y.lon) <= CLOSE_DEG;
+}
+
 describe('regions.geo.json', () => {
   it('is a valid GeoJSON FeatureCollection of well-formed features', () => {
     expect(regions.type).toBe('FeatureCollection');
@@ -188,6 +285,47 @@ describe('regions.geo.json', () => {
     );
   });
 
+  it('leaves only islands below the ~10,000 km2 coarsening floor (non-island small regions merged)', () => {
+    // v0.14 coarsening: every region smaller than the floor that has a
+    // same-country sibling within CLOSE_DEG of its bbox is dissolved into its
+    // nearest such neighbor, so a below-floor region may only survive if it
+    // is an island. Count the non-island residuals below the floor and allow
+    // a small tolerance band (strictly < 1% of all regions) — the coarsener's
+    // fixed point leaves none, and this catches a regression in the merge
+    // logic (e.g. a broken neighbor test that exempts land regions) without
+    // being brittle to one or two borderline cases.
+    const smallNonIslands: string[] = [];
+    let smallIslands = 0;
+    let smallTotal = 0;
+    for (const f of regions.features) {
+      const a = regionArea.get(f.id) ?? 0;
+      if (a >= COARSE_FLOOR_KM2) continue;
+      smallTotal++;
+      const bb = regionBBox.get(f.id)!;
+      let hasNeighbor = false;
+      for (const g of regions.features) {
+        if (g.id === f.id || g.properties.parent !== f.properties.parent) continue;
+        if (bboxesClose(bb, regionBBox.get(g.id)!)) {
+          hasNeighbor = true;
+          break;
+        }
+      }
+      if (hasNeighbor) {
+        smallNonIslands.push(`${f.id} (${Math.round(a)} km2)`);
+      } else smallIslands++;
+    }
+    const frac = regions.features.length ? smallNonIslands.length / regions.features.length : 0;
+    expect(
+      frac < 0.01,
+      `${smallNonIslands.length} non-island region(s) below the ${COARSE_FLOOR_KM2} km2 floor ` +
+        `(${(100 * frac).toFixed(2)}% of all regions): ${smallNonIslands.slice(0, 10).join(', ')}`,
+    ).toBe(true);
+    console.log(
+      `\n  below-floor regions: ${smallTotal} — ${smallIslands} islands (exempt, no same-country neighbor within ${CLOSE_DEG}°)` +
+        ` / ${smallNonIslands.length} non-island residuals (must stay < 1% of all regions)`,
+    );
+  });
+
   it('sums to within [40%, 175%] of each country area (median in [70%, 130%])', () => {
     // The band is a planar-approximation sanity check. It is wider than a naive
     // [40%, 110%] because Natural Earth's admin-1 polygons are digitized