facett-map3d/src/terrain.rs
//! **Real terrain elevation (DEM) for the orbit map** — ported from the flight
//! engine ([`facett_map3d_flight::terrain`]) so the lightweight, deck-embeddable
//! orbit [`crate::Map3D`] can break its flat `y = 0` ground plane and show **real
//! topography** under the extruded city, not a billiard-table quad.
//!
//! `facett-map3d-flight` depends on `facett-map3d` (not the other way round), so the
//! elevation maths cannot be *imported* here without a dependency cycle — it is
//! **ported** (the same SRTM `.hgt` loader + bilinear sampling + procedural fallback),
//! kept byte-for-byte behaviour-compatible so the two crates agree on "how high is the
//! ground at this `(lon, lat)`". The flight crate builds its mesh in the `f64`
//! floating-origin world frame; the orbit map works in the normalised `[-1, 1]`
//! **design box**, so the mesh *builder* lives in [`crate::extrude`] and this module
//! carries only the source-agnostic [`Elevation`] query (`Dem` / `ProceduralTerrain`
//! / `TerrainField`), all pure CPU math, headless + deterministically testable (FC-7).
use std::fs;
use std::io;
use std::path::Path;
/// SRTM `.hgt` **void** sentinel (no-data); replaced by `0.0` on load, counted in
/// [`Dem::void_count`].
pub const HGT_VOID: i16 = -32_768;
/// **Anything that answers "how high is the ground at this `(lon, lat)`?"** — the one
/// query the relief mesh builder ([`crate::extrude::relief_terrain_into`]) goes
/// through, so a real [`Dem`], the [`ProceduralTerrain`] fallback and the
/// [`TerrainField`] enum are interchangeable.
pub trait Elevation {
/// Ground elevation in **metres** at WGS-84 `(lon, lat)` (degrees). Must be total
/// (finite for any input) so the mesh builder never emits a NaN vertex — sampling
/// off a tile clamps to the border, it does not panic.
fn elevation_at(&self, lon_deg: f64, lat_deg: f64) -> f32;
}
// ───────────────────────────────────────────────────────────────────────────────
// Dem — a regular geographic elevation grid + bilinear sampling.
// ───────────────────────────────────────────────────────────────────────────────
/// A **Digital Elevation Model**: a regular `cols × rows` grid of elevation samples
/// (metres) over a geographic rectangle. Row 0 is the **north** edge, column 0 the
/// **west** edge (SRTM `.hgt` convention), stored row-major. Sampled with
/// [`Dem::elevation_at`] (bilinear, edge-clamped).
#[derive(Clone, Debug, PartialEq)]
pub struct Dem {
/// Longitude (°) of column 0 — the **west** edge.
pub west_lon: f64,
/// Latitude (°) of row 0 — the **north** edge.
pub north_lat: f64,
/// Degrees of longitude between adjacent columns (east +).
pub d_lon: f64,
/// Degrees of latitude between adjacent rows, going **south** (positive).
pub d_lat: f64,
/// Number of columns (samples west→east).
pub cols: usize,
/// Number of rows (samples north→south).
pub rows: usize,
/// Row-major heights (metres); `heights[row * cols + col]`, row 0 = north.
pub heights: Vec<f32>,
/// How many samples were the `.hgt` void sentinel (filled with `0.0` on load).
pub void_count: usize,
}
impl Dem {
/// Raw sample at integer `(col, row)`, indices clamped to the grid (edge-extend).
#[must_use]
pub fn sample(&self, col: usize, row: usize) -> f32 {
let c = col.min(self.cols - 1);
let r = row.min(self.rows - 1);
self.heights[r * self.cols + c]
}
/// The DEM's geographic bounds `(min_lon, min_lat, max_lon, max_lat)` in degrees.
#[must_use]
pub fn bounds(&self) -> (f64, f64, f64, f64) {
let max_lon = self.west_lon + self.d_lon * (self.cols.saturating_sub(1)) as f64;
let min_lat = self.north_lat - self.d_lat * (self.rows.saturating_sub(1)) as f64;
(self.west_lon, min_lat, max_lon, self.north_lat)
}
/// Build a DEM from an explicit height grid (row 0 = north, row-major). Panics if
/// `heights.len() != cols * rows`.
#[must_use]
pub fn from_grid(west_lon: f64, north_lat: f64, d_lon: f64, d_lat: f64, cols: usize, rows: usize, heights: Vec<f32>) -> Self {
assert_eq!(heights.len(), cols * rows, "height grid must be cols*rows");
assert!(cols >= 1 && rows >= 1, "grid needs at least one sample per axis");
Self { west_lon, north_lat, d_lon, d_lat, cols, rows, heights, void_count: 0 }
}
/// A **synthetic** DEM: evaluate `f(lon, lat) -> metres` on a `cols × rows` grid
/// spanning `[west_lon, west_lon + span_lon] × [north_lat - span_lat, north_lat]`.
/// Deterministic (no RNG).
#[must_use]
pub fn synthetic(
west_lon: f64,
north_lat: f64,
span_lon: f64,
span_lat: f64,
cols: usize,
rows: usize,
f: impl Fn(f64, f64) -> f32,
) -> Self {
assert!(cols >= 2 && rows >= 2, "a synthetic grid needs ≥2×2 samples");
let d_lon = span_lon / (cols - 1) as f64;
let d_lat = span_lat / (rows - 1) as f64;
let mut heights = Vec::with_capacity(cols * rows);
for r in 0..rows {
let lat = north_lat - d_lat * r as f64;
for c in 0..cols {
let lon = west_lon + d_lon * c as f64;
heights.push(f(lon, lat));
}
}
Self { west_lon, north_lat, d_lon, d_lat, cols, rows, heights, void_count: 0 }
}
/// Parse an **SRTM `.hgt`** byte buffer: `N×N` big-endian `i16` metres over the
/// 1°×1° cell whose **south-west corner** is `(sw_lon, sw_lat)`. `N` is inferred
/// from the length (`2·N²`) so both 1-arc-second (3601) and 3-arc-second (1201)
/// tiles load. Void samples ([`HGT_VOID`]) become `0.0` and are counted.
///
/// # Errors
/// If the buffer length is not `2·N²` for an integer `N ≥ 2`.
pub fn from_hgt_bytes(bytes: &[u8], sw_lon: f64, sw_lat: f64) -> io::Result<Self> {
if bytes.len() < 8 || !bytes.len().is_multiple_of(2) {
return Err(io::Error::new(io::ErrorKind::InvalidData, "hgt: odd/too-short byte length"));
}
let samples = bytes.len() / 2;
let n = (samples as f64).sqrt().round() as usize;
if n < 2 || n * n != samples {
return Err(io::Error::new(io::ErrorKind::InvalidData, format!("hgt: {samples} samples is not a square N×N grid")));
}
let mut heights = Vec::with_capacity(samples);
let mut void_count = 0usize;
for chunk in bytes.as_chunks::<2>().0 {
let v = i16::from_be_bytes([chunk[0], chunk[1]]);
if v == HGT_VOID {
void_count += 1;
heights.push(0.0);
} else {
heights.push(f32::from(v));
}
}
let d = 1.0 / (n - 1) as f64; // the cell is exactly 1°×1°
Ok(Self { west_lon: sw_lon, north_lat: sw_lat + 1.0, d_lon: d, d_lat: d, cols: n, rows: n, heights, void_count })
}
/// Load an SRTM `.hgt` file from disk, taking the SW corner from its **filename**
/// (`N59E017.hgt`, `S12W034.hgt`) via [`parse_hgt_name`].
///
/// # Errors
/// If the filename is not a valid `.hgt` corner name, the file cannot be read, or
/// its bytes are not a square grid.
pub fn load_hgt(path: impl AsRef<Path>) -> io::Result<Self> {
let path = path.as_ref();
let name = path
.file_name()
.and_then(|s| s.to_str())
.ok_or_else(|| io::Error::new(io::ErrorKind::InvalidInput, "hgt: no filename"))?;
let (sw_lon, sw_lat) =
parse_hgt_name(name).ok_or_else(|| io::Error::new(io::ErrorKind::InvalidInput, format!("hgt: bad corner name {name:?}")))?;
let bytes = fs::read(path)?;
Self::from_hgt_bytes(&bytes, sw_lon, sw_lat)
}
}
impl Elevation for Dem {
/// **Bilinear** elevation at `(lon, lat)`, edge-clamped (continuous, no
/// stair-steps; exact for a linear height field).
fn elevation_at(&self, lon_deg: f64, lat_deg: f64) -> f32 {
let fc = ((lon_deg - self.west_lon) / self.d_lon).clamp(0.0, (self.cols - 1) as f64);
let fr = ((self.north_lat - lat_deg) / self.d_lat).clamp(0.0, (self.rows - 1) as f64);
let c0 = fc.floor() as usize;
let r0 = fr.floor() as usize;
let c1 = (c0 + 1).min(self.cols - 1);
let r1 = (r0 + 1).min(self.rows - 1);
let tx = (fc - c0 as f64) as f32;
let ty = (fr - r0 as f64) as f32;
let h00 = self.sample(c0, r0);
let h10 = self.sample(c1, r0);
let h01 = self.sample(c0, r1);
let h11 = self.sample(c1, r1);
let top = h00 + (h10 - h00) * tx;
let bot = h01 + (h11 - h01) * tx;
top + (bot - top) * ty
}
}
/// Parse an SRTM `.hgt` **corner filename** into the SW-corner `(lon, lat)` in degrees:
/// `"N59E017.hgt"` → `(17.0, 59.0)`, `"S12W034"` → `(-34.0, -12.0)`. The `.hgt`
/// extension is optional. `None` if the pattern does not match.
#[must_use]
pub fn parse_hgt_name(name: &str) -> Option<(f64, f64)> {
let stem = name.strip_suffix(".hgt").or_else(|| name.strip_suffix(".HGT")).unwrap_or(name);
let b = stem.as_bytes();
if b.len() < 7 {
return None;
}
let lat_sign = match b[0].to_ascii_uppercase() {
b'N' => 1.0,
b'S' => -1.0,
_ => return None,
};
let lat_deg: f64 = stem.get(1..3)?.parse().ok()?;
let lon_sign = match b[3].to_ascii_uppercase() {
b'E' => 1.0,
b'W' => -1.0,
_ => return None,
};
let lon_deg: f64 = stem.get(4..7)?.parse().ok()?;
Some((lon_sign * lon_deg, lat_sign * lat_deg))
}
// ───────────────────────────────────────────────────────────────────────────────
// ProceduralTerrain — the deterministic synthetic fallback.
// ───────────────────────────────────────────────────────────────────────────────
/// A **procedural** heightfield: smooth, deterministic, analytic topography — the
/// synthetic fallback for CI / when no DEM tile is present. Sums a few sine/cosine
/// octaves so it has hills + ridges at multiple scales; bounded in
/// `[base_m − amplitude_m, base_m + amplitude_m]`, continuous, no clock / RNG (FC-7).
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct ProceduralTerrain {
/// Mean ground height (metres) — the field oscillates around this.
pub base_m: f32,
/// Peak deviation from `base_m` (metres); the field stays within `base ± amplitude`.
pub amplitude_m: f32,
/// Angular frequency of the coarsest octave (radians per degree).
pub freq: f64,
/// Phase seed — different seeds give different (but each deterministic) terrain.
pub seed: u64,
}
impl Default for ProceduralTerrain {
/// A Sweden-flavoured default: gentle ~600 m relief around a 300 m base.
fn default() -> Self {
Self { base_m: 300.0, amplitude_m: 600.0, freq: 0.9, seed: 0x5E_7A_11 }
}
}
impl ProceduralTerrain {
/// A procedural field with a given relief around a base height.
#[must_use]
pub fn new(base_m: f32, amplitude_m: f32, freq: f64, seed: u64) -> Self {
Self { base_m, amplitude_m, freq, seed }
}
}
impl Elevation for ProceduralTerrain {
fn elevation_at(&self, lon_deg: f64, lat_deg: f64) -> f32 {
// Four octaves of a sine/cosine field, each octave sampled in its OWN rotated
// frame; amplitude halving + frequency doubling per octave. `h` normalised to
// [-1, 1] by the amplitude sum → strictly bounded, deterministic, always finite.
//
// The per-octave ROTATION is load-bearing, not decoration. The field used to be
// separable in world axes — `sin(lon·f + φ) · cos(lat·f·1.3 + 0.7φ)` — and a
// separable field degenerates over any viewport whose lon span is small against
// 1/f: the sin factor goes near-constant, `h ≈ C · cos(lat·…)`, and the relief
// collapses into PARALLEL EAST-WEST RIDGES. On korp's flight pane (0.35° span)
// that rendered as a handful of horizontal lines marching up the screen —
// reported live by the owner as "skärmen har 5 horisontella linjer fördelade i
// y-led". Rotating each octave's sample frame by its own phase keeps the field
// deterministic and bounded while leaving no world axis privileged, so no
// camera-aligned banding survives the octave sum.
let mut h = 0.0f64;
let mut amp = 1.0f64;
let mut sum = 0.0f64;
let mut f = self.freq;
for oct in 0..4u64 {
let mix = self.seed.wrapping_add(oct).wrapping_mul(0x9E37_79B9_7F4A_7C15);
let phase = (mix >> 40) as f64 / (1u64 << 24) as f64 * std::f64::consts::TAU;
let (s, c) = phase.sin_cos();
let u = lon_deg * c - lat_deg * s;
let v = lon_deg * s + lat_deg * c;
h += amp * (u * f + phase).sin() * (v * f * 1.3 + phase * 0.7).cos();
sum += amp;
amp *= 0.5;
f *= 2.0;
}
(f64::from(self.base_m) + f64::from(self.amplitude_m) * (h / sum)) as f32
}
}
// ───────────────────────────────────────────────────────────────────────────────
// TerrainField — real DEM or synthetic fallback behind one query.
// ───────────────────────────────────────────────────────────────────────────────
/// The **elevation source** the orbit relief binds to: either a real loaded [`Dem`]
/// or the [`ProceduralTerrain`] fallback, behind one [`Elevation`] query.
#[derive(Clone, Debug, PartialEq)]
pub enum TerrainField {
/// A real Digital Elevation Model (e.g. an SRTM `.hgt` tile).
Dem(Box<Dem>),
/// The deterministic synthetic fallback.
Synthetic(ProceduralTerrain),
}
impl TerrainField {
/// Load a real SRTM `.hgt` tile as the elevation source.
///
/// # Errors
/// Propagates [`Dem::load_hgt`] failures.
pub fn load_hgt(path: impl AsRef<Path>) -> io::Result<Self> {
Dem::load_hgt(path).map(|d| Self::Dem(Box::new(d)))
}
/// The default deterministic synthetic fallback.
#[must_use]
pub fn synthetic() -> Self {
Self::Synthetic(ProceduralTerrain::default())
}
/// **Real-source-or-fallback**: try to load the `.hgt` at `path`; if it is missing
/// or unreadable, fall back to the deterministic [`ProceduralTerrain`]. Returns
/// whether the real tile loaded.
#[must_use]
pub fn load_hgt_or_synthetic(path: impl AsRef<Path>) -> (Self, bool) {
match Self::load_hgt(path) {
Ok(field) => (field, true),
Err(_) => (Self::synthetic(), false),
}
}
/// **The orbit map's terrain source, resolved from the environment.** If
/// `$FACETT_HGT` points at a readable SRTM `.hgt` tile it becomes the real DEM;
/// otherwise the deterministic [`ProceduralTerrain`] fallback (so a fixture-free
/// run still shows real relief, honestly labelled by [`Self::is_real`]).
#[must_use]
pub fn from_env() -> (Self, bool) {
match std::env::var("FACETT_HGT") {
Ok(path) if !path.is_empty() => Self::load_hgt_or_synthetic(path),
_ => (Self::synthetic(), false),
}
}
/// Whether this is a real DEM (vs the synthetic fallback).
#[must_use]
pub fn is_real(&self) -> bool {
matches!(self, Self::Dem(_))
}
}
impl Elevation for TerrainField {
fn elevation_at(&self, lon_deg: f64, lat_deg: f64) -> f32 {
match self {
Self::Dem(d) => d.elevation_at(lon_deg, lat_deg),
Self::Synthetic(p) => p.elevation_at(lon_deg, lat_deg),
}
}
}
#[cfg(test)]
mod tests {
use super::*;
/// Bilinear sampling is exact on a linear ramp (pins the interpolation maths).
#[test]
fn bilinear_is_exact_on_a_linear_ramp() {
let (cols, rows) = (11usize, 7usize);
let (a, b) = (10.0f32, -3.0f32);
let mut heights = vec![0.0f32; cols * rows];
for r in 0..rows {
for c in 0..cols {
heights[r * cols + c] = a * c as f32 + b * r as f32;
}
}
let dem = Dem::from_grid(10.0, 60.0, 0.1, 0.05, cols, rows, heights);
for &(fc, fr) in &[(2.5f64, 3.0f64), (0.0, 0.0), (9.9, 5.5), (7.25, 1.75)] {
let lon = 10.0 + fc * 0.1;
let lat = 60.0 - fr * 0.05;
let want = a as f64 * fc + b as f64 * fr;
let got = dem.elevation_at(lon, lat);
assert!((f64::from(got) - want).abs() < 1e-3, "bilinear exact at ({fc},{fr}): want {want}, got {got}");
}
}
/// Sampling outside the grid clamps to the border (never panics / NaNs).
#[test]
fn sampling_clamps_outside_the_grid() {
let dem = Dem::from_grid(0.0, 1.0, 1.0, 1.0, 2, 2, vec![0.0, 5.0, 10.0, 15.0]);
assert_eq!(dem.elevation_at(-100.0, 100.0), 0.0);
assert_eq!(dem.elevation_at(100.0, -100.0), 15.0);
assert!(dem.elevation_at(0.5, 0.5).is_finite());
}
/// The real `.hgt` byte loader round-trips a known grid + a void sentinel.
#[test]
fn hgt_loader_round_trips_a_known_tile() {
let n = 4usize;
let mut bytes = Vec::with_capacity(n * n * 2);
for r in 0..n {
for c in 0..n {
let v: i16 = if (r, c) == (1, 2) { HGT_VOID } else { (r * 100 + c) as i16 };
bytes.extend_from_slice(&v.to_be_bytes());
}
}
let dem = Dem::from_hgt_bytes(&bytes, 17.0, 59.0).expect("valid hgt");
assert_eq!((dem.cols, dem.rows), (4, 4));
assert_eq!(dem.void_count, 1);
assert!((dem.d_lon - 1.0 / 3.0).abs() < 1e-12);
assert_eq!(dem.elevation_at(17.0, 60.0), 0.0);
assert_eq!(dem.sample(3, 2), 203.0);
assert_eq!(dem.sample(2, 1), 0.0);
assert!(Dem::from_hgt_bytes(&bytes[..bytes.len() - 2], 17.0, 59.0).is_err());
}
/// `.hgt` corner-name parsing handles both hemispheres + the optional extension.
#[test]
fn hgt_names_parse_both_hemispheres() {
assert_eq!(parse_hgt_name("N59E017.hgt"), Some((17.0, 59.0)));
assert_eq!(parse_hgt_name("S12W034"), Some((-34.0, -12.0)));
assert_eq!(parse_hgt_name("N00E000.hgt"), Some((0.0, 0.0)));
assert_eq!(parse_hgt_name("garbage"), None);
assert_eq!(parse_hgt_name("X59E017.hgt"), None);
}
/// The procedural fallback is bounded in `base ± amp`, deterministic, and has relief.
#[test]
fn procedural_terrain_is_bounded_and_deterministic() {
let p = ProceduralTerrain::new(300.0, 600.0, 0.9, 42);
let mut lo = f32::INFINITY;
let mut hi = f32::NEG_INFINITY;
for i in 0..200 {
let lon = 10.0 + i as f64 * 0.07;
let lat = 55.0 + i as f64 * 0.06;
let h = p.elevation_at(lon, lat);
assert_eq!(h, p.elevation_at(lon, lat), "deterministic");
assert!(h.is_finite());
lo = lo.min(h);
hi = hi.max(h);
}
assert!(lo >= 300.0 - 600.0 - 1e-3 && hi <= 300.0 + 600.0 + 1e-3, "bounded: [{lo},{hi}]");
assert!(hi - lo > 100.0, "procedural terrain has relief (Δ={})", hi - lo);
}
/// **★ The procedural field carries no world-axis banding over a flight-pane
/// window.** The RED this pins: the pre-rotation field was separable in world
/// axes, so over korp's 0.35°-wide flight window the lon factor went
/// near-constant and the relief collapsed to a pure function of LATITUDE —
/// horizontal ridge bands on screen ("5 horisontella linjer", live report).
///
/// The probe is the degeneracy itself, not a pixel: over that window, fit the
/// rank-poorest model there is — "elevation depends only on lat" (per-row means)
/// — and measure how much variance it explains. The OLD field scores ~0.97
/// (measured red before the rotation landed); a field with no privileged axis
/// must leave real variance across longitude too. Deleting the per-octave
/// rotation in `elevation_at` turns this red again.
#[test]
fn the_synthetic_field_is_not_a_function_of_latitude_alone() {
let p = ProceduralTerrain::default();
const N: usize = 64;
let (lon0, lat0, dlon, dlat) = (17.0f64, 59.0f64, 0.35f64, 0.385f64);
let mut grid = [[0.0f64; N]; N];
let mut total_mean = 0.0f64;
for (r, row) in grid.iter_mut().enumerate() {
for (c, cell) in row.iter_mut().enumerate() {
let lon = lon0 + dlon * c as f64 / (N - 1) as f64;
let lat = lat0 + dlat * r as f64 / (N - 1) as f64;
*cell = f64::from(p.elevation_at(lon, lat));
total_mean += *cell;
}
}
total_mean /= (N * N) as f64;
let mut var_total = 0.0f64;
let mut var_lat_only = 0.0f64; // variance the "lat alone" model explains
for row in &grid {
let row_mean = row.iter().sum::<f64>() / N as f64;
var_lat_only += (row_mean - total_mean).powi(2) * N as f64;
for h in row {
var_total += (h - total_mean).powi(2);
}
}
assert!(var_total > 0.0, "precondition: the window has relief at all");
let explained = var_lat_only / var_total;
assert!(
explained < 0.85,
"the 'elevation is lat alone' model explains {:.1} % of the window's \
variance — that is the horizontal-band degeneracy, not terrain",
explained * 100.0,
);
}
/// `load_hgt_or_synthetic` / `from_env` fall back to synthetic when no tile is present.
#[test]
fn missing_tile_falls_back_to_synthetic() {
let (field, real) = TerrainField::load_hgt_or_synthetic("/no/such/tile/N59E017.hgt");
assert!(!real);
assert!(!field.is_real());
assert!(field.elevation_at(17.0, 59.0).is_finite());
}
/// `TerrainField` forwards `elevation_at` to a loaded real DEM.
#[test]
fn terrain_field_dispatches_to_the_real_dem() {
let dem = Dem::from_grid(0.0, 1.0, 1.0, 1.0, 2, 2, vec![10.0, 20.0, 30.0, 40.0]);
let field = TerrainField::Dem(Box::new(dem.clone()));
assert!(field.is_real());
assert_eq!(field.elevation_at(0.25, 0.75), dem.elevation_at(0.25, 0.75));
}
}