Re-coarsen RU-ARK Arctic-coast band into RU-PEO + RU-YAM (#41)

#41: re-partition the RU-ARK mainland's 28-part lat>=68 Arctic-coast band into two new regions via a new turf-free, idempotent pipeline stage scripts/split-arkhangel-band.mjs (npm run split:regions:band), mirroring the #31 split: - RU-PEO (Pechora, 13 parts, ~37,678 km2) -> key pechora-qwhs - RU-YAM (Yamal, 15 parts, ~51,722 km2) -> key yamal-gosw - RU-ARK keeps its id/name/key arkhangel-sk-lhbk (ownership-safe); the 5-part lat<68 core (~305,229 km2) stays RU-ARK. Pure re-partition: no geometry moved (byte-identical union). The two new features are inserted immediately after RU-ARK-SZ. Manifests regenerated with the real scripts (1,798 -> 1,800 keys + adjacency; adjacency pairs 4,471 -> 4,473). Fixed-point dataset unchanged except feature count (11,794 parts / 11,795 rings / 164,801 verts / 1 KZ-KZY hole). Four dataset/manifest tests updated (RUS 82 -> 84 sub-regions). Version 0.25.0 -> 0.26.0. Fixes #41
2026-10-14 19:26:49 +00:00
parent 867b2eb0a1
commit f6ad554823
12 changed files with 382 additions and 70 deletions
@@ -0,0 +1,247 @@
+#!/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, …)`,
+);