facett-map3d-flight/src/sweden.rs
//! **Sweden extent + a Web-Mercator slippy tile grid** — the spatial index that
//! turns "fly the whole of Sweden" into a bounded, prioritisable set of tiles.
//!
//! ## What this is
//! - [`SwedenExtent`] — the WGS-84 bounding box of Sweden (lon ~10.9–24.2,
//! lat ~55.3–69.1) and the set of slippy tiles that cover it at a zoom `z`.
//! - [`SpatialKey`] — a slippy / Web-Mercator tile `{z, x, y}` **consistent with
//! [`crate::origin::WorldPos::from_lonlat`]**: a key's tile is derived from the
//! SAME spherical Web-Mercator projection the floating origin uses, so a streamed
//! tile's geometry lands exactly where the camera expects it (no half-pixel skew
//! between the tile grid and the render frame).
//! - [`tiles_in_view`] — given the camera-snapped [`FloatingOrigin`], a zoom, and a
//! world-metre radius, the in-view tiles **priority-ordered nearest-first** by
//! [`WorldPos::ground_distance`] — the frustum/priority seam the §5 streamer loads
//! against.
//!
//! Pure `f64` math (no IO, no GPU, no clock) → every function here is directly
//! inject-assert testable (FC-7).
use crate::origin::{EARTH_RADIUS_M, FloatingOrigin, WorldPos};
/// Half the Web-Mercator world span in metres: `π·R`. The projected world is the
/// square `[-WORLD_HALF, WORLD_HALF]²` (x east, y north), so a zoom-`z` tile spans
/// `2·WORLD_HALF / 2^z` metres on a side. Matches [`WorldPos::from_lonlat`]'s frame.
pub const WORLD_HALF_M: f64 = std::f64::consts::PI * EARTH_RADIUS_M;
/// A **slippy / Web-Mercator tile key** `{z, x, y}`. `x` grows east (0 at the
/// antimeridian), `y` grows south (0 at the north Mercator limit) — the standard
/// OSM/XYZ convention, so these keys name the SAME tiles an `osm2geoparquet` /
/// katana-osm tiler writes.
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash, PartialOrd, Ord)]
pub struct SpatialKey {
/// Zoom level (0 = whole world in one tile).
pub z: u8,
/// Tile column, `0..2^z`, growing east.
pub x: u32,
/// Tile row, `0..2^z`, growing south.
pub y: u32,
}
impl SpatialKey {
#[must_use]
pub const fn new(z: u8, x: u32, y: u32) -> Self {
Self { z, x, y }
}
/// Number of tiles per axis at this zoom (`2^z`), as `f64` for the projections.
#[must_use]
fn n(z: u8) -> f64 {
(1u64 << z) as f64
}
/// The tile whose Web-Mercator cell contains a `WorldPos` (x,y in metres). This
/// is the single projection all the other constructors go through, so the grid
/// is consistent with the floating-origin world frame **by construction**.
#[must_use]
pub fn from_world(w: WorldPos, z: u8) -> Self {
let n = Self::n(z);
// Normalise x∈[-WORLD_HALF,WORLD_HALF] → u∈[0,1] (east), y → v∈[0,1] (south).
let u = (w.x + WORLD_HALF_M) / (2.0 * WORLD_HALF_M);
let v = (WORLD_HALF_M - w.y) / (2.0 * WORLD_HALF_M);
let clampi = |t: f64| (t * n).floor().clamp(0.0, n - 1.0) as u32;
Self::new(z, clampi(u), clampi(v))
}
/// The tile containing a WGS-84 `(lon, lat)` at zoom `z`, via the same
/// [`WorldPos::from_lonlat`] projection the renderer uses (lat clamped to the
/// Mercator limit inside `from_lonlat`).
#[must_use]
pub fn from_lonlat(lon_deg: f64, lat_deg: f64, z: u8) -> Self {
Self::from_world(WorldPos::from_lonlat(lon_deg, lat_deg, 0.0), z)
}
/// The tile **centre** in the `f64` world frame (altitude 0) — the point the
/// priority metric measures to.
#[must_use]
pub fn center_world(self) -> WorldPos {
let n = Self::n(self.z);
let u = (self.x as f64 + 0.5) / n;
let v = (self.y as f64 + 0.5) / n;
let x = u * 2.0 * WORLD_HALF_M - WORLD_HALF_M;
let y = WORLD_HALF_M - v * 2.0 * WORLD_HALF_M;
WorldPos::new(x, y, 0.0)
}
/// The tile centre as WGS-84 `(lon, lat)` degrees (inverse spherical Mercator).
#[must_use]
pub fn center_lonlat(self) -> (f64, f64) {
let n = Self::n(self.z);
let lon = (self.x as f64 + 0.5) / n * 360.0 - 180.0;
let lat_rad = (std::f64::consts::PI * (1.0 - 2.0 * (self.y as f64 + 0.5) / n)).sinh().atan();
(lon, lat_rad.to_degrees())
}
/// Tile side length in world (Mercator) metres at this zoom.
#[must_use]
pub fn size_metres(self) -> f64 {
2.0 * WORLD_HALF_M / Self::n(self.z)
}
}
/// **Sweden's WGS-84 extent** — the bounding box the streamer flies over. The
/// bounds are the country's mainland envelope (a touch generous so coastal tiles
/// are covered): lon `10.9 .. 24.2`, lat `55.3 .. 69.1`.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct SwedenExtent;
impl SwedenExtent {
/// West / minimum longitude (°). Öresund / west coast.
pub const MIN_LON: f64 = 10.9;
/// East / maximum longitude (°). Torne valley.
pub const MAX_LON: f64 = 24.2;
/// South / minimum latitude (°). Smygehuk.
pub const MIN_LAT: f64 = 55.3;
/// North / maximum latitude (°). Treriksröset.
pub const MAX_LAT: f64 = 69.1;
/// `(min_lon, min_lat, max_lon, max_lat)` in degrees.
#[must_use]
pub const fn bbox() -> (f64, f64, f64, f64) {
(Self::MIN_LON, Self::MIN_LAT, Self::MAX_LON, Self::MAX_LAT)
}
/// Whether a WGS-84 `(lon, lat)` lies inside the Sweden extent.
#[must_use]
pub fn contains(lon_deg: f64, lat_deg: f64) -> bool {
(Self::MIN_LON..=Self::MAX_LON).contains(&lon_deg) && (Self::MIN_LAT..=Self::MAX_LAT).contains(&lat_deg)
}
/// The inclusive tile-index range `((x0, y0), (x1, y1))` covering the extent at
/// zoom `z`. NB `y0` (north) ≤ `y1` (south) because tile-y grows southward.
#[must_use]
pub fn tile_range(z: u8) -> ((u32, u32), (u32, u32)) {
// Corners: NW = (min_lon, max_lat), SE = (max_lon, min_lat).
let nw = SpatialKey::from_lonlat(Self::MIN_LON, Self::MAX_LAT, z);
let se = SpatialKey::from_lonlat(Self::MAX_LON, Self::MIN_LAT, z);
((nw.x, nw.y), (se.x, se.y))
}
/// Every tile that covers the Sweden extent at zoom `z`, row-major. This is the
/// full working set the streamer pages through as the camera flies the country.
#[must_use]
pub fn tiles(z: u8) -> Vec<SpatialKey> {
let ((x0, y0), (x1, y1)) = Self::tile_range(z);
let mut out = Vec::new();
for y in y0..=y1 {
for x in x0..=x1 {
out.push(SpatialKey::new(z, x, y));
}
}
out
}
/// The extent's centre as a `WorldPos` at altitude 0 (a handy camera anchor).
#[must_use]
pub fn center_world(altitude_m: f64) -> WorldPos {
let lon = (Self::MIN_LON + Self::MAX_LON) * 0.5;
let lat = (Self::MIN_LAT + Self::MAX_LAT) * 0.5;
WorldPos::from_lonlat(lon, lat, altitude_m)
}
}
/// **Frustum / priority tile selection.** Given the camera-snapped
/// [`FloatingOrigin`] (its anchor is the query centre — the camera ground point),
/// a zoom `z`, and a `radius_m` in world (Mercator) metres, return every Sweden
/// tile whose centre lies within `radius_m` of the camera, **ordered nearest-first**
/// by [`WorldPos::ground_distance`] (ties broken by `(y, x)` for determinism).
///
/// This is the seam the §5 out-of-core streamer loads against: request the front of
/// this list first (nearest tiles), evict tiles that fall off its tail. Pure +
/// deterministic — a scripted camera path yields a bit-reproducible tile order.
#[must_use]
pub fn tiles_in_view(origin: FloatingOrigin, z: u8, radius_m: f64) -> Vec<SpatialKey> {
let center = origin.anchor;
let radius = radius_m.max(0.0);
// Bound the candidate window: a world-metre box of ±radius around the centre,
// converted to a tile-index window, intersected with Sweden's tile range.
let lo = SpatialKey::from_world(WorldPos::new(center.x - radius, center.y + radius, 0.0), z); // NW
let hi = SpatialKey::from_world(WorldPos::new(center.x + radius, center.y - radius, 0.0), z); // SE
let ((sx0, sy0), (sx1, sy1)) = SwedenExtent::tile_range(z);
let x0 = lo.x.max(sx0);
let x1 = hi.x.min(sx1);
let y0 = lo.y.max(sy0);
let y1 = hi.y.min(sy1);
let r2 = radius * radius;
let mut scored: Vec<(f64, SpatialKey)> = Vec::new();
if x0 <= x1 && y0 <= y1 {
for y in y0..=y1 {
for x in x0..=x1 {
let key = SpatialKey::new(z, x, y);
let d2 = {
let c = key.center_world();
let dx = c.x - center.x;
let dy = c.y - center.y;
dx * dx + dy * dy
};
if d2 <= r2 {
scored.push((d2, key));
}
}
}
}
// Nearest-first; deterministic tie-break on the key (z,x,y are Ord).
scored.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal).then(a.1.cmp(&b.1)));
scored.into_iter().map(|(_, k)| k).collect()
}
#[cfg(test)]
mod tests {
use super::*;
const STHLM: (f64, f64) = (18.07, 59.33); // Stockholm
/// **INJECT-ASSERT (the grid is consistent with the render frame):** projecting a
/// `(lon,lat)` to a tile and back to the tile centre stays within half a tile of
/// the input, and the round-trip through `WorldPos` agrees — the slippy grid and
/// the floating-origin world frame describe the SAME Mercator.
#[test]
fn tile_grid_round_trips_through_worldpos() {
for &(lon, lat) in &[STHLM, (13.0, 55.6), (20.3, 67.85), (11.97, 57.7)] {
for z in [6u8, 10, 13] {
let key = SpatialKey::from_lonlat(lon, lat, z);
// Re-deriving the tile from the tile's own centre is a fixed point.
let re = SpatialKey::from_world(key.center_world(), z);
assert_eq!(re, key, "tile→centre→tile is a fixed point at z={z}");
// The tile the point falls in contains it: centre within ~1 tile.
let (clon, clat) = key.center_lonlat();
let w = WorldPos::from_lonlat(lon, lat, 0.0);
let c = key.center_world();
let half = key.size_metres() * 0.5 + 1.0;
assert!((w.x - c.x).abs() <= half && (w.y - c.y).abs() <= half, "point sits inside its tile ({clon},{clat})");
}
}
}
/// **INJECT-ASSERT (the grid IS the OSM/XYZ slippy grid — authoritative tile refs):**
/// [`SpatialKey::from_lonlat`] must reproduce the canonical slippy tile numbers a
/// standard `z/x/y` tiler / OSM tile server assigns, so GeoParquet tiles named the
/// OSM way line up with the render frame. Checked two ways: (1) against the
/// closed-form OSM formula `x=⌊(λ+180)/360·2ᶻ⌋`,
/// `y=⌊(1−asinh(tan φ)/π)/2·2ᶻ⌋` recomputed *independently* here, and (2) against
/// hand-verified integer references for Swedish landmarks.
#[test]
fn spatial_key_matches_canonical_osm_slippy_tiles() {
// The canonical slippy-tile algorithm (OSM wiki), computed independently of the
// metre-based `from_world` path so this is a real cross-check, not a tautology.
fn canonical(lon: f64, lat: f64, z: u8) -> (u32, u32) {
let n = (1u64 << z) as f64;
let x = ((lon + 180.0) / 360.0 * n).floor();
let lat_r = lat.to_radians();
let y = ((1.0 - lat_r.tan().asinh() / std::f64::consts::PI) / 2.0 * n).floor();
(x as u32, y as u32)
}
let landmarks = [
("Stockholm", 18.0686, 59.3293),
("Kiruna", 20.2253, 67.8558),
("Sigtuna", 17.7206, 59.6169),
("Göteborg", 11.9746, 57.7089),
("Malmö", 13.0038, 55.6050),
];
for &(name, lon, lat) in &landmarks {
for z in [0u8, 6, 8, 10, 12, 13, 14, 16] {
let key = SpatialKey::from_lonlat(lon, lat, z);
let (cx, cy) = canonical(lon, lat, z);
assert_eq!((key.x, key.y), (cx, cy), "{name} z{z}: SpatialKey must equal the OSM slippy tile");
}
}
// Hand-verified integer references (recomputed offline from the OSM formula) —
// pins the numbers so a future refactor can't drift them silently.
let refs = [
("Stockholm", 18.0686, 59.3293, 10u8, 563u32, 301u32),
("Stockholm", 18.0686, 59.3293, 13, 4507, 2409),
("Kiruna", 20.2253, 67.8558, 10, 569, 246),
("Kiruna", 20.2253, 67.8558, 13, 4556, 1969),
("Sigtuna", 17.7206, 59.6169, 12, 2249, 1198),
("Göteborg", 11.9746, 57.7089, 12, 2184, 1239),
("Malmö", 13.0038, 55.6050, 8, 137, 80),
];
for &(name, lon, lat, z, x, y) in &refs {
let key = SpatialKey::from_lonlat(lon, lat, z);
assert_eq!((key.x, key.y), (x, y), "{name} z{z} tile is the known reference {x},{y}");
}
}
/// **INJECT-ASSERT (tile-centre inverse round-trips + tile geometry is exact):** a
/// tile's `center_lonlat` (inverse spherical Mercator) lands back in the same tile,
/// its `center_world` agrees with projecting that lon/lat, the side length is the
/// exact `2·π·R/2ᶻ`, and neighbouring tile centres are exactly one side apart.
#[test]
fn tile_centre_inverse_and_size_are_exact() {
for z in [6u8, 10, 13] {
let side = 2.0 * WORLD_HALF_M / (1u64 << z) as f64;
for &(lon, lat) in &[STHLM, (20.2253, 67.8558), (13.0038, 55.6050)] {
let key = SpatialKey::from_lonlat(lon, lat, z);
// side length is the exact Mercator tile span.
assert!((key.size_metres() - side).abs() < 1e-6, "z{z} side is 2πR/2^z");
// center_lonlat is the inverse-Mercator of the tile centre → same tile.
let (clon, clat) = key.center_lonlat();
assert_eq!(SpatialKey::from_lonlat(clon, clat, z), key, "tile centre falls in its own tile");
// center_world and projecting center_lonlat describe the same point.
let cw = key.center_world();
let pw = WorldPos::from_lonlat(clon, clat, 0.0);
assert!((cw.x - pw.x).abs() < 1e-3 && (cw.y - pw.y).abs() < 1e-3, "center_world == project(center_lonlat)");
// The east neighbour's centre is exactly one tile side east.
let east = SpatialKey::new(z, key.x + 1, key.y);
assert!((east.center_world().x - cw.x - side).abs() < 1e-6, "east neighbour is one side over");
}
}
}
/// **INJECT-ASSERT (Sweden bbox covers all its tiles):** every tile
/// `SwedenExtent::tiles(z)` returns has its centre inside the extent's tile
/// window, the tile of every corner + Stockholm is present, and the count is the
/// full `width × height` grid — so the streamer's working set spans the country.
#[test]
fn sweden_bbox_covers_all_its_tiles() {
let z = 8u8;
let tiles = SwedenExtent::tiles(z);
let ((x0, y0), (x1, y1)) = SwedenExtent::tile_range(z);
let w = (x1 - x0 + 1) as usize;
let h = (y1 - y0 + 1) as usize;
assert_eq!(tiles.len(), w * h, "tiles() is the full {w}×{h} covering grid");
// Every returned key is inside the covering window.
assert!(tiles.iter().all(|k| (x0..=x1).contains(&k.x) && (y0..=y1).contains(&k.y)), "all tiles inside the window");
// Corners + Stockholm land in the set.
for &(lon, lat) in &[STHLM, (SwedenExtent::MIN_LON + 0.1, SwedenExtent::MIN_LAT + 0.1), (SwedenExtent::MAX_LON - 0.1, SwedenExtent::MAX_LAT - 0.1)] {
assert!(SwedenExtent::contains(lon, lat), "corner/city inside extent");
let key = SpatialKey::from_lonlat(lon, lat, z);
assert!(tiles.contains(&key), "tiles() contains {key:?} for {lon},{lat}");
}
// A point OUTSIDE Sweden is not claimed as inside.
assert!(!SwedenExtent::contains(2.35, 48.85), "Paris is not in the Sweden extent");
}
/// **INJECT-ASSERT (tiles_in_view selects Stockholm-area tiles nearest-first):** a
/// camera over Stockholm gets the Stockholm tile FIRST, every returned tile is
/// within the radius, the distances are non-decreasing (priority order), and a
/// far city (Kiruna) is excluded by a modest radius.
#[test]
fn tiles_in_view_priority_is_nearest_first_over_stockholm() {
let z = 11u8;
let cam = WorldPos::from_lonlat(STHLM.0, STHLM.1, 4000.0);
let origin = FloatingOrigin::at(cam);
let radius = 30_000.0; // 30 km of Mercator metres
let view = tiles_in_view(origin, z, radius);
assert!(!view.is_empty(), "the camera over Stockholm sees tiles");
// The nearest tile is the one the camera sits in.
let here = SpatialKey::from_world(cam, z);
assert_eq!(view[0], here, "nearest-first: the camera's own tile leads");
// Priority order: ground_distance is non-decreasing down the list.
let mut prev = -1.0;
for k in &view {
let d = k.center_world().ground_distance(cam);
assert!(d <= radius + k.size_metres(), "tile within the view radius");
assert!(d + 1e-6 >= prev, "distances are non-decreasing (priority order)");
prev = d;
}
// A distant city is outside a 30 km radius (Kiruna ~1000 km north).
let kiruna = SpatialKey::from_lonlat(20.22, 67.85, z);
assert!(!view.contains(&kiruna), "Kiruna is not in a 30 km Stockholm view");
// Widening the radius only adds tiles (monotone) and keeps nearest-first.
let wide = tiles_in_view(origin, z, radius * 3.0);
assert!(wide.len() >= view.len(), "a wider radius sees at least as many tiles");
assert_eq!(wide[0], here, "still nearest-first when widened");
#[cfg(feature = "testmatrix")]
facett_core::testmatrix::emit(
"facett-osm3d-sweden",
"tiles_in_view_priority",
view[0] == here && !view.contains(&kiruna),
&format!("z={z} tiles={} first_is_camera_tile=true excludes_far_city=true", view.len()),
);
}
/// INJECT-ASSERT: `tiles_in_view` clips to the Sweden extent — a camera parked in
/// the sea west of the country still only ever returns Sweden tiles.
#[test]
fn tiles_in_view_is_clipped_to_sweden() {
let z = 9u8;
// Out in the North Sea, west of Sweden.
let cam = WorldPos::from_lonlat(6.0, 58.0, 2000.0);
let view = tiles_in_view(FloatingOrigin::at(cam), z, 200_000.0);
let ((x0, y0), (x1, y1)) = SwedenExtent::tile_range(z);
assert!(view.iter().all(|k| (x0..=x1).contains(&k.x) && (y0..=y1).contains(&k.y)), "every returned tile is a Sweden tile");
}
}