History/scripts/split-arkhangel-band.mjs

#!/usr/bin/env node
/**
 * scripts/split-arkhangel-band.mjs — split the RU-ARK Arctic-coast band into
 * two new regions, RU-PEO ("Pechora", the SW group) and RU-YAM ("Yamal", the
 * NE group) (issue #41, a follow-up to #31).
 *
 * The #31 split (`scripts/split-arkhangel.mjs`) left `RU-ARK` as a 33-part
 * MultiPolygon (~393,256 km²). Measured on that feature (centroid latitude /
 * planar area), the parts divide cleanly by centroid latitude:
 *
 *   lat < 68            4 parts  ~305,000 km²  (Arkhangel'sk proper + 3 slivers)
 *   66 ≤ lat < 68       1 part   ~47 km²
 *   lat ≥ 68            28 parts ~88,000 km²   (the "Arctic-coast band")
 *
 * The band is really *other* Arctic oblasts (Nenets / Pechora / Yamal / Taimyr
 * side) mis-lumped into Arkhangel'sk during the v0.14 coarsening. This stage
 * re-coarsens it: the 28 band parts are split into two new regions by a
 * Voronoi of the two band parts with the largest planar area (the SW anchor
 * vs the NE anchor), while the 5 lat < 68 parts stay in `RU-ARK`.
 *
 *   RU-ARK   (keeps its id and name)  — the 5 lat < 68 parts (original order)
 *                                          → still the `arkhangel-sk-lhbk` key
 *   RU-PEO   "Pechora"                — the band's SW Voronoi cell (13 parts)
 *   RU-YAM   "Yamal"                  — the band's NE Voronoi cell (15 parts)
 *
 * Deterministic, turf-free threshold spec (all data-derived — no magic
 * coordinates; the same planar model the tests/adjacency build use):
 *   1. centroid of a part = the mean of its outer-ring vertices (verified to
 *      give the same band assignment as turf's area-weighted centroid, 0 flips).
 *   2. band = the parts whose centroid lat ≥ 68 (asserted to be exactly 28).
 *   3. anchors = the 2 band parts with the largest planar area (≈ 47,223 and
 *      ≈ 32,907 km², indices 2 and 1 in the 33-part feature); anchorW = the
 *      anchor with the smaller centroid longitude.
 *   4. each non-anchor band part → the nearer anchor by km-weighted centroid
 *      distance, `dist(a,b) = hypot(Δlon·111.32·cos(midLat), Δlat·110.574)`;
 *      a tie → west (anchorW). Each anchor → its own group. This yields
 *      PEO = {1,4,5,9,11,12,13,14,16,17,18,20,22} and
 *      YAM = {2,6,15,19,21,23,24,25,26,27,28,29,30,31,32} (0-based indices).
 *
 * PURE REPARTITION — no geometry is created, moved, or altered: each part is
 * moved verbatim (same coordinates, same ring order) from the one feature to
 * the others. The union of the three output geometries is byte-identical to
 * the input feature's parts, so no border moves and no new overlap is
 * possible (the same guarantee as #31). `RU-ARK` keeps its stable feature id,
 * so its region key (`arkhangel-sk-lhbk`, FNV-1a of the id) and every existing
 * `Nation.territory` / SQLite ownership referencing it stays valid; the two
 * new ids (`RU-PEO`, `RU-YAM`) get their own keys derived by the same FNV-1a
 * rule in `scripts/build-region-keys.mjs` (`pechora-qwhs`, `yamal-gosw` —
 * run it after this stage).
 *
 * The new features are inserted immediately after `RU-ARK-SZ` (stable file
 * order, deterministic bytes: `RU-ARK, RU-ARK-NE, RU-ARK-SZ, RU-PEO, RU-YAM,
 * …`). Each new feature's `properties` copies `RU-ARK`'s (only `name` and
 * `parent` exist in this dataset) with `name` overridden; `parent` stays RUS.
 * Each feature's geometry is a MultiPolygon of its assigned parts (in
 * original order), or a single Polygon if it ended up with one part (it won't
 * here, but the builder handles it).
 *
 * Pipeline position: AFTER `split-arkhangel` and BEFORE `build-region-keys`
 * (the manifests must be rebuilt after any geo change):
 *
 *   generate → coarsen → cleanup → smooth → fill-gaps → split-arkhangel
 *   → **split-arkhangel-band** → build-region-keys → build-region-adjacency
 *
 * Idempotent: when RU-PEO and RU-YAM already exist the file is left untouched
 * (no write at all), so a full re-run of the pipeline reproduces byte-identical
 * output.
 *
 * Usage:  node scripts/split-arkhangel-band.mjs [inFile [outFile]]
 *          (or `npm run split:regions:band`)
 *   inFile  default: public/regions.geo.json
 *   outFile default: inFile (edit in place)
 *
 * Deliberately turf-free: no `@turf/*` import (turf is only ever imported by
 * the geo-generator `.mjs` stages, never by this re-partition or by any `.ts`).
 */
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, '..');
const inFile = process.argv[2] ?? path.join(root, 'public', 'regions.geo.json');
const outFile = process.argv[3] ?? inFile;

// A part whose centroid lat is ≥ this is part of the Arctic-coast band; the
// lat < 68 core stays in RU-ARK (issue #41).
const BAND_MIN_LAT = 68;

// The post-#31 mainland feature is a 33-part MultiPolygon; any other shape
// means the file is not the state this stage expects → fail loud.
const EXPECTED_PARTS = 33;
// The Arctic-coast band is exactly 28 of those parts (lat ≥ 68).
const EXPECTED_BAND = 28;

// Planar km model (matches the test's area model and the adjacency build):
// x = lon · 111.32 · cos(midLat), y = lat · 110.574.
const KM_LON = 111.32;
const KM_LAT = 110.574;
const rad = (deg) => (deg * Math.PI) / 180;

/**
 * Mean of a part's outer-ring vertices — a turf-free centroid. Verified to give
 * the same band assignment as turf's area-weighted centroid (0 flips).
 */
function centroidOf(part) {
  const ring = part[0];
  let lon = 0;
  let lat = 0;
  for (const [lo, la] of ring) {
    lon += lo;
    lat += la;
  }
  return [lon / ring.length, lat / ring.length];
}

/** Planar area of a part (km²), shoelace over its outer ring at its mean latitude. */
function areaKm2(part) {
  const ring = part[0];
  let lat = 0;
  for (const [, la] of ring) lat += la;
  const midLat = lat / ring.length;
  const kx = KM_LON * Math.cos(rad(midLat));
  let a = 0;
  for (let i = 0; i < ring.length; i++) {
    const [lon1, lat1] = ring[i];
    const [lon2, lat2] = ring[(i + 1) % ring.length];
    a += lon1 * kx * (lat2 * KM_LAT) - lon2 * kx * (lat1 * KM_LAT);
  }
  return Math.abs(a / 2);
}

/** Planar km distance between two centroids (midpoint-latitude weighted). */
function kmDist(a, b) {
  const dLon = b[0] - a[0];
  const dLat = b[1] - a[1];
  const midLat = (a[1] + b[1]) / 2;
  return Math.hypot(dLon * KM_LON * Math.cos(rad(midLat)), dLat * KM_LAT);
}

// --- load --------------------------------------------------------------------
const geo = JSON.parse(readFileSync(inFile, 'utf8'));
if (geo?.type !== 'FeatureCollection' || !Array.isArray(geo.features)) {
  throw new Error(`${inFile} is not a FeatureCollection`);
}
const features = geo.features;

// --- idempotency: already split? ----------------------------------------------
const ids = new Set(features.map((f) => String(f.id ?? '')));
if (ids.has('RU-PEO') && ids.has('RU-YAM')) {
  console.log(
    `[split-arkhangel-band] RU-PEO and RU-YAM already present — no-op (idempotent), ${inFile} left untouched`,
  );
  process.exit(0);
}

// --- load + validate RU-ARK ---------------------------------------------------
const idx = features.findIndex((f) => f.id === 'RU-ARK');
if (idx < 0) throw new Error('RU-ARK not found in ' + inFile);
const ark = features[idx];
if (ark.geometry?.type !== 'MultiPolygon') {
  throw new Error(
    `RU-ARK geometry is ${ark.geometry?.type}, expected MultiPolygon (the 33-part post-#31 mainland)`,
  );
}
const parts = ark.geometry.coordinates;
if (parts.length !== EXPECTED_PARTS) {
  throw new Error(`RU-ARK has ${parts.length} parts, expected ${EXPECTED_PARTS} (the post-#31 mainland)`);
}

// --- band selection (centroid lat ≥ 68) and the lat < 68 core -----------------
const partsByIndex = parts.map((p, i) => ({ i, c: centroidOf(p) }));
const band = partsByIndex.filter((o) => o.c[1] >= BAND_MIN_LAT);
if (band.length !== EXPECTED_BAND) {
  throw new Error(`expected exactly ${EXPECTED_BAND} Arctic-coast band parts (lat ≥ ${BAND_MIN_LAT}), found ${band.length}`);
}
const core = partsByIndex.filter((o) => o.c[1] < BAND_MIN_LAT);

// --- anchors: the two band parts with the largest planar area -----------------
// (≈ 47,223 and ≈ 32,907 km², indices 2 and 1 in the 33-part feature).
const byArea = [...band].sort((a, b) => areaKm2(parts[b.i]) - areaKm2(parts[a.i]));
const big = byArea[0];
const small = byArea[1];
const anchorW = big.c[0] <= small.c[0] ? big : small; // the west (SW) anchor
const anchorE = big.c[0] <= small.c[0] ? small : big; // the east (NE) anchor

// --- assignment (Voronoi by km-weighted centroid distance; tie → west) --------
const peo = new Set([anchorW.i]);
const yam = new Set([anchorE.i]);
for (const o of band) {
  if (o.i === anchorW.i || o.i === anchorE.i) continue;
  const dW = kmDist(o.c, anchorW.c);
  const dE = kmDist(o.c, anchorE.c);
  (dW <= dE ? peo : yam).add(o.i); // a tie goes to the west (anchorW) group
}

// --- rebuild geometries (parts verbatim, original order; MultiPolygon/Polygon) --
const geom = (ps) =>
  ps.length > 1 ? { type: 'MultiPolygon', coordinates: ps } : { type: 'Polygon', coordinates: ps[0] };

const coreParts = [...core].sort((a, b) => a.i - b.i).map((o) => parts[o.i]);
const peoParts = [...peo].sort((a, b) => a - b).map((i) => parts[i]);
const yamParts = [...yam].sort((a, b) => a - b).map((i) => parts[i]);

ark.geometry = geom(coreParts);

const mkFeature = (id, name, ps) => ({
  type: 'Feature',
  id,
  // copy every property RU-ARK has (currently just { name, parent }), overriding name
  properties: { ...ark.properties, name },
  geometry: geom(ps),
});
const peoF = mkFeature('RU-PEO', 'Pechora', peoParts);
const yamF = mkFeature('RU-YAM', 'Yamal', yamParts);

// --- report -------------------------------------------------------------------
const km2 = (ps) => ps.reduce((s, p) => s + areaKm2(p), 0);
console.log(`[split-arkhangel-band] input: ${inFile}`);
console.log(
  `  RU-ARK  (core, lat < 68)   ${coreParts.length} parts (of ${parts.length}), ≈ ${Math.round(km2(coreParts)).toLocaleString('en-US')} km²`,
);
console.log(
  `  RU-PEO  (Pechora, SW)      ${peoParts.length} parts, ≈ ${Math.round(km2(peoParts)).toLocaleString('en-US')} km²  [${[...peo].sort((a, b) => a - b).join(',')}]`,
);
console.log(
  `  RU-YAM  (Yamal, NE)        ${yamParts.length} parts, ≈ ${Math.round(km2(yamParts)).toLocaleString('en-US')} km²  [${[...yam].sort((a, b) => a - b).join(',')}]`,
);
if (peoParts.length === 0 || yamParts.length === 0) {
  throw new Error(
    `threshold mis-split: Pechora=${peoParts.length} parts, Yamal=${yamParts.length} parts — check the band before shipping`,
  );
}

// --- insert at a stable position (immediately after RU-ARK-SZ) and write -------
// RU-ARK-SZ must exist (the #31 split ran first); if it doesn't, this stage is
// being run out of order → fail loud.
const szIdx = features.findIndex((f) => f.id === 'RU-ARK-SZ');
if (szIdx < 0) {
  throw new Error('RU-ARK-SZ not found — the #31 split (scripts/split-arkhangel.mjs) must run before this stage');
}
features.splice(szIdx + 1, 0, peoF, yamF);
writeFileSync(outFile, JSON.stringify(geo));
console.log(
  `[split-arkhangel-band] wrote ${outFile} — features ${features.length - 2} → ${features.length}` +
    ` (order …, RU-ARK, RU-ARK-NE, RU-ARK-SZ, RU-PEO, RU-YAM, …)`,
);