diff --git a/cockpit/public/body.manifest.json b/cockpit/public/body.manifest.json index 2b20d4b87..7ae58e69e 100644 --- a/cockpit/public/body.manifest.json +++ b/cockpit/public/body.manifest.json @@ -1,7 +1,8 @@ { "helix_latest": "body.20260629c.v6helix.soa.gz", "osm_latest": "berlin.helix.soa.gz", - "iceland_latest": "iceland.helix.soa.gz", + "iceland_latest": "iceland_dem.helix.soa.gz", + "iceland_note": "iceland_dem.helix.soa.gz = real Terrarium DEM heightfield (z9 keyless terrain-RGB, 1152x896 grid = 1,032,192 verts / 2,060,290 tris, equirectangular, true-scale elevation 0..0.0067 in [-1,1], sea level -> y=0). Solid landmass, replaces the sparse iceland.helix.soa.gz scatter (kept on disk, unreferenced). Baked by geo/src/bin/iceland_dem.rs via bso2::encode_mesh_bso2; DEM fetched by scripts/fetch_iceland_dem.py.", "note": "20260629c re-emit from soa_v2 (geometry identical to 20260629b): 39 connective structures (ligaments / tendons / interosseous membranes / fascia / retinacula / iliotibial tract) reclassified out of the ORGAN and SKIN layers into the now-live CONNECTIVE layer 7 — they were FMA-filed under /viscera/solid_organ/ligament_organ, so the is_a walk tagged them viscus->organ and they floated in the organ view as tan limb-shaped strays (interosseous membrane of leg/forearm, calcaneal tendon, long plantar ligament). Carries the 20260629b fixes (teeth->skeleton, per-vessel slicer-fill diameter). BSO2 ver 6 = F16 pos + Signed360 NORMAL + HXFL trailer; Gouraud per-vertex shading. Published to fma-body-soa-v3-v1; Dockerfile pulls same-origin.", "verts": 4283525 } diff --git a/cockpit/public/iceland_dem.helix.soa.gz b/cockpit/public/iceland_dem.helix.soa.gz new file mode 100644 index 000000000..a9848139e Binary files /dev/null and b/cockpit/public/iceland_dem.helix.soa.gz differ diff --git a/cockpit/src/BodyHelix.tsx b/cockpit/src/BodyHelix.tsx index cb6a6b09b..7a660a304 100644 --- a/cockpit/src/BodyHelix.tsx +++ b/cockpit/src/BodyHelix.tsx @@ -606,7 +606,15 @@ export default function BodyHelix() { background: active ? '#1c2738' : '#0e1219', color: active ? '#cdd9e5' : '#6b7686', font: '12px ui-monospace, monospace', }); const q = query.trim().toLowerCase(); - const groups = LAYERS.map((l) => ({ + // Geo scenes: the terrain bake stamps a single layer — present it as "terrain" and drop the + // empty anatomy layers, so a map never reads as skin/muscle/skeleton (the layer *id* is kept, so + // the show/hide toggle still filters the real geometry). Anatomy keeps the full LAYERS taxonomy. + const geoUI = Boolean(new URLSearchParams(window.location.search).get('scene') ?? pathScene()); + const activeLayers = + geoUI && d + ? LAYERS.filter((l) => d.conceptList.some((c) => c.layer === l.id)).map((l) => ({ ...l, name: 'terrain', color: '#7c8f5c' })) + : LAYERS; + const groups = activeLayers.map((l) => ({ l, items: d ? d.conceptList.filter((c) => c.layer === l.id && (!q || c.name.toLowerCase().includes(q))) : [], })).filter((g) => g.items.length > 0 || !q); @@ -637,7 +645,7 @@ export default function BodyHelix() {
- {LAYERS.map((l) => ( + {activeLayers.map((l) => ( diff --git a/geo/Cargo.toml b/geo/Cargo.toml index 1293b88b7..018480ed6 100644 --- a/geo/Cargo.toml +++ b/geo/Cargo.toml @@ -53,6 +53,12 @@ required-features = ["osm"] name = "osm_helix" required-features = ["osm", "helix"] +# iceland_dem → BSO2 /helix wire from a Terrarium DEM heightfield (no osm reader +# at runtime; needs `helix` for bso2 + ndarray + lance-graph-contract). +[[bin]] +name = "iceland_dem" +required-features = ["helix"] + [profile.release] opt-level = 3 codegen-units = 1 diff --git a/geo/src/bin/iceland_dem.rs b/geo/src/bin/iceland_dem.rs new file mode 100644 index 000000000..483d73e0b --- /dev/null +++ b/geo/src/bin/iceland_dem.rs @@ -0,0 +1,230 @@ +//! `iceland_dem` — bake a real Iceland DEM heightfield into the BSO2 `/helix` +//! wire so `/ice` renders SOLID terrain (not the sparse scatter). +//! +//! Input is a `.demgrid` produced by `scripts/fetch_iceland_dem.py` (keyless AWS +//! Terrarium terrain-RGB tiles → stitched, downsampled elevation raster). We turn +//! the W×H elevation grid into a shared-vertex heightfield mesh (two triangles +//! per cell, per-vertex normals from the terrain gradient), normalize it to the +//! viewer's `[-1,1]` frame TRUE-SCALE (elevation in the same metric units as the +//! horizontal extent — the shader applies `uExag`, so we do NOT pre-exaggerate), +//! and run it through the SAME `bso2::encode_mesh_bso2` path the OSM bake uses +//! (F16 pos + Signed360 normals + HXFL trailer). Sea level (elev ≤ 0) → y = 0, +//! which the shader paints as ocean. +//! +//! ```text +//! usage: iceland_dem (writes .soa + .blocks) +//! ``` +//! Gzip `.soa` afterwards (BodyHelix fetches a gzip'd artifact). + +use geo_hhtl::bso2::{encode_mesh_bso2, MeshConcept, CLASSID_GEO_V3}; +use geo_hhtl::hhtl::point_to_hhtl4; +use geo_hhtl::osm_read::M_PER_DEG; + +use lance_graph_contract::canonical_node::{NodeGuid, TailVariant}; + +/// Deep key zoom for the 4-tier HHTL cascade (matches the OSM bake's KEY_ZOOM). +const KEY_ZOOM: u32 = 32; +/// Terrain basin tier (family) — a natural-area basin, distinct from the OSM +/// building/road/area families (1/2/3); 4 = "terrain". +const FAMILY_TERRAIN: u32 = 4; + +struct Dem { + w: usize, + h: usize, + west: f64, + east: f64, + lats: Vec, // len h, row 0 = north + elev: Vec, // len w*h, row-major, row 0 = north, col 0 = west (metres) +} + +fn read_demgrid(path: &str) -> Result> { + let b = std::fs::read(path)?; + if b.len() < 24 || &b[0..4] != b"DEMG" { + return Err("not a DEMG file".into()); + } + let u32_at = |o: usize| u32::from_le_bytes(b[o..o + 4].try_into().unwrap()); + let f64_at = |o: usize| f64::from_le_bytes(b[o..o + 8].try_into().unwrap()); + let ver = u32_at(4); + if ver != 1 { + return Err(format!("unsupported DEMG version {ver}").into()); + } + let w = u32_at(8) as usize; + let h = u32_at(12) as usize; + let west = f64_at(16); + let east = f64_at(24); + let mut o = 32; + let mut lats = Vec::with_capacity(h); + for _ in 0..h { + lats.push(f64_at(o)); + o += 8; + } + let n = w * h; + let mut elev = Vec::with_capacity(n); + for _ in 0..n { + elev.push(f32::from_le_bytes(b[o..o + 4].try_into().unwrap())); + o += 4; + } + Ok(Dem { w, h, west, east, lats, elev }) +} + +fn main() { + let mut args = std::env::args().skip(1); + let (Some(input), Some(output)) = (args.next(), args.next()) else { + eprintln!("usage: iceland_dem "); + std::process::exit(2); + }; + let dem = match read_demgrid(&input) { + Ok(d) => d, + Err(e) => { + eprintln!("iceland_dem: {e}"); + std::process::exit(1); + } + }; + let (w, h) = (dem.w, dem.h); + eprintln!( + "DEM {w}x{h} verts · lon {:.4}..{:.4} · lat {:.4}..{:.4}", + dem.west, dem.east, dem.lats[0], dem.lats[h - 1] + ); + + // ── Equirectangular projection about the grid centre (matches osm_read). ── + let lon0 = (dem.west + dem.east) * 0.5; + let lat0 = (dem.lats[0] + dem.lats[h - 1]) * 0.5; + let cos_lat0 = lat0.to_radians().cos(); + let lon_at = |c: usize| dem.west + (dem.east - dem.west) * c as f64 / (w - 1).max(1) as f64; + // metric (x = east, z = north) in metres, y = elevation (sea clamped to 0). + let mut mx = vec![0.0f32; w * h]; + let mut mz = vec![0.0f32; w * h]; + let mut my = vec![0.0f32; w * h]; + for r in 0..h { + let z = ((dem.lats[r] - lat0) * M_PER_DEG) as f32; + for c in 0..w { + let x = ((lon_at(c) - lon0) * cos_lat0 * M_PER_DEG) as f32; + let i = r * w + c; + mx[i] = x; + mz[i] = z; + my[i] = dem.elev[i].max(0.0); // sea level & bathymetry → 0 + } + } + + // ── Normalize to [-1,1] by horizontal extent; elevation TRUE-SCALE (÷half). ── + let (mut lox, mut hix, mut loz, mut hiz) = (f32::MAX, f32::MIN, f32::MAX, f32::MIN); + for i in 0..w * h { + lox = lox.min(mx[i]); + hix = hix.max(mx[i]); + loz = loz.min(mz[i]); + hiz = hiz.max(mz[i]); + } + let cx = (lox + hix) * 0.5; + let cz = (loz + hiz) * 0.5; + let half = ((hix - lox).max(hiz - loz) * 0.5).max(1.0); + let inv = 1.0 / half; + // DISPLAY frame (Dx, Dy=up=elevation, Dz). No vertical exaggeration — the + // shader raises the true-scale island by uExag=10 + the Kurvenlineal. + let mut pos = vec![[0.0f32; 3]; w * h]; + for i in 0..w * h { + pos[i] = [(mx[i] - cx) * inv, my[i] * inv, (mz[i] - cz) * inv]; + } + + // ── Per-vertex normals from the terrain gradient (display frame). ── + // du ≈ +east (col+), dv ≈ +north (row-, since row 0 is north). cross(dv,du) + // points +Dy on flat ground and tilts with the slope. One-sided at edges. + let idx = |r: usize, c: usize| r * w + c; + let sub = |a: [f32; 3], b: [f32; 3]| [a[0] - b[0], a[1] - b[1], a[2] - b[2]]; + let cross = |a: [f32; 3], b: [f32; 3]| { + [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]] + }; + let mut nrm = vec![[0.0f32, 1.0, 0.0]; w * h]; + for r in 0..h { + let rn = r.saturating_sub(1); // north neighbour (smaller row) + let rs = (r + 1).min(h - 1); // south neighbour + for c in 0..w { + let cw = c.saturating_sub(1); + let ce = (c + 1).min(w - 1); + let du = sub(pos[idx(r, ce)], pos[idx(r, cw)]); // east tangent + let dv = sub(pos[idx(rn, c)], pos[idx(rs, c)]); // north tangent + let mut n = cross(dv, du); + let m = (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt(); + if m > 1e-12 { + n = [n[0] / m, n[1] / m, n[2] / m]; + if n[1] < 0.0 { + n = [-n[0], -n[1], -n[2]]; // keep the up-facing hemisphere + } + } else { + n = [0.0, 1.0, 0.0]; + } + nrm[idx(r, c)] = n; + } + } + + // ── Triangles: two per grid cell, shared vertices (watertight surface). ── + let mut tris: Vec<[u32; 3]> = Vec::with_capacity((w - 1) * (h - 1) * 2); + for r in 0..h - 1 { + for c in 0..w - 1 { + let a = idx(r, c) as u32; + let b = idx(r, c + 1) as u32; + let d = idx(r + 1, c) as u32; + let e = idx(r + 1, c + 1) as u32; + // Wound so the THREE.js face normal (v1-v0)×(v2-v0) points +Dy (up): + // the terrain top is the FrontSide face an overhead camera sees. + tris.push([a, b, d]); + tris.push([b, e, d]); + } + } + + // ── Concepts = grid rows (contiguous vertex ranges; the per-concept + // max-diameter clamp stays a thin full-width strip → drops nothing). ── + let rows: Vec = (0..h).flat_map(|r| std::iter::repeat(r as u32).take(w)).collect(); + let mid = w / 2; + let mut concepts = Vec::with_capacity(h); + for r in 0..h { + let v_start = (r * w) as u32; + let mut c = [0.0f32; 3]; + for i in r * w..(r + 1) * w { + c[0] += pos[i][0]; + c[1] += pos[i][1]; + c[2] += pos[i][2]; + } + let invn = 1.0 / w as f32; + let centroid = [c[0] * invn, c[1] * invn, c[2] * invn]; + let key4 = point_to_hhtl4(lon_at(mid), dem.lats[r], KEY_ZOOM); + let key = NodeGuid::mint_for( + TailVariant::V3, + CLASSID_GEO_V3, + key4.heel, + key4.hip, + key4.twig, + key4.leaf, + FAMILY_TERRAIN, + r as u32, + ); + concepts.push(MeshConcept { + key, + layer: 4, // default-on geo toggle group + label: 0, // names[0] = "iceland-terrain" + centroid, + v_start, + v_count: w as u32, + }); + } + + let labels = br#"{"names":["iceland-terrain"]}"#; + let (soa, blocks) = encode_mesh_bso2(&pos, &nrm, &rows, &tris, &concepts, labels); + + let nc = u32::from_le_bytes(soa[6..10].try_into().unwrap()); + let nv = u32::from_le_bytes(soa[10..14].try_into().unwrap()); + let nt = u32::from_le_bytes(soa[14..18].try_into().unwrap()); + eprintln!( + "BSO2: {nc} concepts · {nv} verts · {nt} tris · {} B soa · {} B blocks", + soa.len(), + blocks.len() + ); + + let blocks_path = format!("{}.blocks", output.strip_suffix(".soa").unwrap_or(&output)); + if let Err(e) = + std::fs::write(&output, &soa).and_then(|()| std::fs::write(&blocks_path, &blocks)) + { + eprintln!("iceland_dem: write: {e}"); + std::process::exit(1); + } + eprintln!("wrote {output} + {blocks_path}"); +} diff --git a/geo/src/bso2.rs b/geo/src/bso2.rs index 0ace32048..fa09b7864 100644 --- a/geo/src/bso2.rs +++ b/geo/src/bso2.rs @@ -325,6 +325,128 @@ pub fn encode_bso2(scene: &GeoScene) -> (Vec, Vec) { (o, blocks) } +/// One concept (structure) for the generic mesh encoder [`encode_mesh_bso2`]: +/// its already-minted key, toggle/colour byte, label index, DISPLAY-frame +/// centroid `(Dx, Dy=up, Dz)`, and its **contiguous** vertex range. +pub struct MeshConcept { + /// The minted `NodeGuid` (geo classid + HHTL tiers). + pub key: NodeGuid, + /// Layer toggle byte (`4` = the default-on geo group, as the OSM bake uses). + pub layer: u8, + /// Label index into `labels_json`'s `names[]`. + pub label: u32, + /// Concept centroid in the DISPLAY frame `(Dx, Dy=up, Dz)`. + pub centroid: [f32; 3], + /// First vertex index owned by this concept. + pub v_start: u32, + /// Vertex count (the concept owns `[v_start, v_start + v_count)`). + pub v_count: u32, +} + +/// Encode an **arbitrary display-frame triangle mesh** into the exact BSO2 ver-6 +/// `/helix` wire (F16 pos + `Signed360` normals + HXFL trailer) that +/// `BodyHelix.tsx` decodes — the generalization of [`encode_bso2`] for meshes +/// that are NOT OSM footprints (e.g. a DEM heightfield grid). The byte layout is +/// identical to [`encode_bso2`]'s so the same reader renders it unchanged. +/// +/// `pos`/`nrm` are already the normalized DISPLAY frame `(Dx, Dy=up, Dz)`; +/// positions are written as `srcPos = (-Dx, Dz, Dy)` (F16) and normals via +/// [`signed360`], the same remap [`encode_bso2`] applies, so BodyHelix's +/// `(-x, z, y)` decode lands the mesh upright. `rows[i]` is vertex `i`'s concept +/// (the client keys its per-vertex layer/clamp off this array, NOT `v_range`). +/// Returns `(bso2_bytes, blocks_sidecar_bytes)`. +#[must_use] +pub fn encode_mesh_bso2( + pos: &[[f32; 3]], + nrm: &[[f32; 3]], + rows: &[u32], + tris: &[[u32; 3]], + concepts: &[MeshConcept], + labels_json: &[u8], +) -> (Vec, Vec) { + let nv = pos.len(); + let nt = tris.len(); + let nc = concepts.len(); + assert_eq!(nrm.len(), nv, "nrm len must match pos"); + assert_eq!(rows.len(), nv, "rows len must match pos"); + + let mut o = Vec::with_capacity(nc * 40 + nv * 22 + nt * 12 + labels_json.len() + 64); + o.extend_from_slice(b"BSO2"); + o.extend_from_slice(&6u16.to_le_bytes()); + o.extend_from_slice(&(nc as u32).to_le_bytes()); + o.extend_from_slice(&(nv as u32).to_le_bytes()); + o.extend_from_slice(&(nt as u32).to_le_bytes()); + // per-concept: guid | material | layer | label | centroid | (v_start,v_count) + for c in concepts { + o.extend_from_slice(c.key.as_bytes()); + } + o.extend(concepts.iter().map(|_| 4u8)); // material (doppler default) + o.extend(concepts.iter().map(|c| c.layer)); // layer (toggle/colour byte) + for c in concepts { + o.extend_from_slice(&c.label.to_le_bytes()); + } + for c in concepts { + // Same (-Dx, Dz, Dy) pre-remap as the vertex positions so focus targets + // land upright (BodyHelix applies (-x, z, y) to this column too). + for v in [-c.centroid[0], c.centroid[2], c.centroid[1]] { + o.extend_from_slice(&v.to_le_bytes()); + } + } + for c in concepts { + o.extend_from_slice(&c.v_start.to_le_bytes()); + o.extend_from_slice(&c.v_count.to_le_bytes()); + } + // per-vertex: pos 3×F16 (srcPos = -Dx, Dz, Dy) | helix 6B (Signed360) | row u32. + for p in pos { + for &v in &[-p[0], p[2], p[1]] { + o.extend_from_slice(&F16::from_f32(v).0.to_le_bytes()); + } + } + for n in nrm { + o.extend_from_slice(&signed360(*n)); + } + for r in rows { + o.extend_from_slice(&r.to_le_bytes()); + } + for t in tris { + for idx in t { + o.extend_from_slice(&idx.to_le_bytes()); + } + } + // labels_json (names[]) + materials_json (unused by the reader). + o.extend_from_slice(&(labels_json.len() as u32).to_le_bytes()); + o.extend_from_slice(labels_json); + let materials = b"{}"; + o.extend_from_slice(&(materials.len() as u32).to_le_bytes()); + o.extend_from_slice(materials); + // HXFL trailer (12 B, must be the final bytes): the rim-dequant floor. We take + // the analytic `end = 255` Signed360 branch, so these are unused, but the + // well-formed ver-6 artifact carries the tag BodyHelix reads from the tail. + o.extend_from_slice(b"HXFL"); + o.extend_from_slice(&(-2.256_794_5f32).to_le_bytes()); + o.extend_from_slice(&11.535_854f32.to_le_bytes()); + + // .blocks sidecar: per-concept (center3, radius) in DISPLAY space (unused by + // geo scenes — LOD is disabled — but emitted to match the OSM path). + let mut blocks = Vec::with_capacity(nc * 16); + for c in concepts { + let mut rad = 0.0f32; + let end = (c.v_start + c.v_count) as usize; + for v in &pos[c.v_start as usize..end.min(nv)] { + let d = (v[0] - c.centroid[0]).powi(2) + + (v[1] - c.centroid[1]).powi(2) + + (v[2] - c.centroid[2]).powi(2); + rad = rad.max(d); + } + for v in c.centroid { + blocks.extend_from_slice(&v.to_le_bytes()); + } + blocks.extend_from_slice(&rad.sqrt().max(1e-4).to_le_bytes()); + } + + (o, blocks) +} + #[cfg(test)] mod tests { use super::*; diff --git a/scripts/fetch_iceland_dem.py b/scripts/fetch_iceland_dem.py new file mode 100755 index 000000000..7a3052c22 --- /dev/null +++ b/scripts/fetch_iceland_dem.py @@ -0,0 +1,153 @@ +#!/usr/bin/env python3 +"""Fetch a keyless Iceland DEM heightfield from AWS Terrarium terrain-RGB tiles. + +Terrarium tiles (`elevation-tiles-prod`, public, no key) encode elevation in the +PNG's RGB: ``elev_m = R*256 + G + B/256 - 32768``. We fetch the slippy-tile grid +covering Iceland at a chosen zoom, stitch into one elevation raster, block-mean +downsample to a tractable heightfield, and emit a compact ``.demgrid`` binary the +`iceland_dem` Rust baker reads (no HTTP dep in the geo crate). + +Grid file layout (little-endian): + "DEMG" 4 bytes magic + version u32 = 1 + W u32 columns (west->east, col 0 = west edge) + H u32 rows (north->south, row 0 = north edge) + west_lon f64 + east_lon f64 (lon is LINEAR across columns — WebMercator x is linear in lon) + lat[H] f64 latitude of each row centre (row 0 = north); non-linear in row + elev[H*W] f32 metres, row-major, row 0 = north, col 0 = west + +usage: + python3 scripts/fetch_iceland_dem.py OUT.demgrid [--zoom 9] [--downsample 4] +""" +import io +import math +import struct +import sys +import urllib.request +from concurrent.futures import ThreadPoolExecutor + +import numpy as np +from PIL import Image + +# Iceland bounding box (slightly padded so the whole island + a sea margin is in). +LON_W, LON_E = -25.0, -13.0 +LAT_S, LAT_N = 63.0, 67.0 + +TILE = 256 +BASE = "https://s3.amazonaws.com/elevation-tiles-prod/terrarium/{z}/{x}/{y}.png" +UA = {"User-Agent": "q2-iceland-dem/1.0 (keyless terrarium bake)"} + + +def lon_to_tilex(lon, z): + return (lon + 180.0) / 360.0 * (1 << z) + + +def lat_to_tiley(lat, z): + r = math.radians(lat) + return (1.0 - math.log(math.tan(r) + 1.0 / math.cos(r)) / math.pi) / 2.0 * (1 << z) + + +def tiley_to_lat(ty, z): + """Inverse WebMercator: fractional tile-y (at zoom z) -> latitude (deg).""" + n = math.pi * (1.0 - 2.0 * ty / (1 << z)) + return math.degrees(math.atan(math.sinh(n))) + + +def fetch_tile(z, x, y, tries=4): + url = BASE.format(z=z, x=x, y=y) + last = None + for _ in range(tries): + try: + req = urllib.request.Request(url, headers=UA) + data = urllib.request.urlopen(req, timeout=60).read() + im = Image.open(io.BytesIO(data)).convert("RGB") + return x, y, np.asarray(im, dtype=np.float32) + except Exception as e: # noqa: BLE001 + last = e + print(f" tile {z}/{x}/{y} FAILED: {last!r}", file=sys.stderr) + return x, y, None + + +def main(): + if len(sys.argv) < 2: + print(__doc__) + sys.exit(2) + out = sys.argv[1] + zoom = 9 + down = 4 + for i, a in enumerate(sys.argv): + if a == "--zoom": + zoom = int(sys.argv[i + 1]) + if a == "--downsample": + down = int(sys.argv[i + 1]) + + x0 = int(math.floor(lon_to_tilex(LON_W, zoom))) + x1 = int(math.floor(lon_to_tilex(LON_E, zoom))) + # tile-y grows south, so LAT_N (north) -> smaller y. + y0 = int(math.floor(lat_to_tiley(LAT_N, zoom))) + y1 = int(math.floor(lat_to_tiley(LAT_S, zoom))) + xs = list(range(x0, x1 + 1)) + ys = list(range(y0, y1 + 1)) + print(f"zoom {zoom}: x {x0}..{x1} ({len(xs)} tiles), y {y0}..{y1} ({len(ys)} tiles) " + f"= {len(xs) * len(ys)} tiles", file=sys.stderr) + + Wtiles, Htiles = len(xs), len(ys) + stitched = np.zeros((Htiles * TILE, Wtiles * TILE), dtype=np.float32) + + jobs = [(zoom, x, y) for y in ys for x in xs] + got, missing = 0, 0 + with ThreadPoolExecutor(max_workers=16) as ex: + for x, y, rgb in ex.map(lambda t: fetch_tile(*t), jobs): + gx = x - x0 + gy = y - y0 + if rgb is None: + missing += 1 + continue + elev = rgb[:, :, 0] * 256.0 + rgb[:, :, 1] + rgb[:, :, 2] / 256.0 - 32768.0 + stitched[gy * TILE:(gy + 1) * TILE, gx * TILE:(gx + 1) * TILE] = elev + got += 1 + print(f"fetched {got} tiles, {missing} missing", file=sys.stderr) + if missing > len(jobs) // 10: + print("too many missing tiles — aborting", file=sys.stderr) + sys.exit(1) + + # Geographic extent of the stitched tile grid (tile edges, exact). + west = xs[0] / (1 << zoom) * 360.0 - 180.0 + east = (xs[-1] + 1) / (1 << zoom) * 360.0 - 180.0 + Hpx, Wpx = stitched.shape + + # Block-mean downsample by `down` (crop to a multiple first). + Hc = (Hpx // down) * down + Wc = (Wpx // down) * down + st = stitched[:Hc, :Wc].reshape(Hc // down, down, Wc // down, down).mean(axis=(1, 3)) + Hout, Wout = st.shape + print(f"stitched {Wpx}x{Hpx} px -> downsampled {Wout}x{Hout} verts " + f"({Wout * Hout:,} verts)", file=sys.stderr) + + # Per-row latitude: row r centre sits at tile-y = y0 + (r+0.5)*down/TILE. + lats = np.empty(Hout, dtype=np.float64) + for r in range(Hout): + ty = y0 + (r + 0.5) * down / TILE + lats[r] = tiley_to_lat(ty, zoom) + + emin, emax = float(st.min()), float(st.max()) + land = int((st > 0.5).sum()) + print(f"elev range {emin:.1f}..{emax:.1f} m; land verts {land:,} " + f"({100 * land / st.size:.1f}%)", file=sys.stderr) + + with open(out, "wb") as f: + f.write(b"DEMG") + f.write(struct.pack("