umsh_regiondb/
morton.rs

1//! The fixed lookup grid.
2//!
3//! Longitude and latitude are each mapped onto 16 bits and interleaved into a
4//! 32-bit Z-order key, so that a quadtree cell at any depth is one contiguous
5//! range of keys and a lookup is a single indexed query.
6//!
7//! Every operation here has an exact counterpart in the Python builder
8//! (`tools/regiondb-build/src/regiondb_build/morton.py`). Both are written as
9//! the same two floating-point operations rather than in whichever idiom each
10//! language finds natural, because a position that lands in different cells on
11//! different platforms would be a lookup that disagrees with itself.
12
13/// Depth of the grid: 16 bits per axis, so cells are roughly 600 m by 300 m at
14/// the equator.
15pub const MAX_DEPTH: u32 = 16;
16
17/// Cells per axis at maximum depth.
18pub const GRID: f64 = 65536.0;
19
20/// A coordinate that cannot be placed on the grid.
21#[derive(Clone, Copy, Debug, PartialEq, Eq)]
22pub enum MortonError {
23    /// Longitude or latitude was NaN or infinite.
24    NotFinite,
25    /// Latitude was outside `[-90, 90]`.
26    LatitudeOutOfRange,
27}
28
29impl core::fmt::Display for MortonError {
30    fn fmt(&self, formatter: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
31        match self {
32            Self::NotFinite => write!(formatter, "coordinate is not finite"),
33            Self::LatitudeOutOfRange => write!(formatter, "latitude is outside [-90, 90]"),
34        }
35    }
36}
37
38impl std::error::Error for MortonError {}
39
40/// Wrap longitude into `[-180, 180)`.
41///
42/// Exactly +180 wraps to -180: the two name the same meridian, and picking one
43/// keeps the grid a partition rather than an overlap.
44pub fn normalize_longitude(longitude: f64) -> Result<f64, MortonError> {
45    if !longitude.is_finite() {
46        return Err(MortonError::NotFinite);
47    }
48    Ok(longitude - 360.0 * ((longitude + 180.0) / 360.0).floor())
49}
50
51/// Validate latitude, which is clamped into the grid rather than wrapped.
52pub fn check_latitude(latitude: f64) -> Result<f64, MortonError> {
53    if !latitude.is_finite() {
54        return Err(MortonError::NotFinite);
55    }
56    if !(-90.0..=90.0).contains(&latitude) {
57        return Err(MortonError::LatitudeOutOfRange);
58    }
59    Ok(latitude)
60}
61
62/// Grid coordinates of the cell containing a position.
63pub fn cell_xy(latitude: f64, longitude: f64) -> Result<(u32, u32), MortonError> {
64    let longitude = normalize_longitude(longitude)?;
65    let latitude = check_latitude(latitude)?;
66    let x = ((longitude + 180.0) / 360.0 * GRID).floor();
67    let y = ((latitude + 90.0) / 180.0 * GRID).floor();
68    Ok((clamp_axis(x), clamp_axis(y)))
69}
70
71fn clamp_axis(value: f64) -> u32 {
72    if value <= 0.0 {
73        0
74    } else if value >= GRID - 1.0 {
75        GRID as u32 - 1
76    } else {
77        value as u32
78    }
79}
80
81/// Interleave two 16-bit values, `x` into the even bits.
82pub fn interleave(x: u32, y: u32) -> u32 {
83    spread(x) | (spread(y) << 1)
84}
85
86fn spread(value: u32) -> u32 {
87    let mut value = value & 0xFFFF;
88    value = (value | (value << 8)) & 0x00FF_00FF;
89    value = (value | (value << 4)) & 0x0F0F_0F0F;
90    value = (value | (value << 2)) & 0x3333_3333;
91    value = (value | (value << 1)) & 0x5555_5555;
92    value
93}
94
95/// The maximum-depth Morton key for a position.
96pub fn key(latitude: f64, longitude: f64) -> Result<u32, MortonError> {
97    let (x, y) = cell_xy(latitude, longitude)?;
98    Ok(interleave(x, y))
99}
100
101#[cfg(test)]
102mod tests {
103    use super::*;
104
105    #[test]
106    fn longitude_180_names_the_same_meridian_as_minus_180() {
107        assert_eq!(normalize_longitude(180.0).unwrap(), -180.0);
108        assert_eq!(key(0.0, 180.0).unwrap(), key(0.0, -180.0).unwrap());
109    }
110
111    #[test]
112    fn longitudes_wrap_rather_than_clamp() {
113        assert!((normalize_longitude(181.0).unwrap() - -179.0).abs() < 1e-9);
114        assert!((normalize_longitude(-181.0).unwrap() - 179.0).abs() < 1e-9);
115        assert!((normalize_longitude(540.0).unwrap() - -180.0).abs() < 1e-9);
116    }
117
118    #[test]
119    fn latitude_90_lands_in_the_last_row() {
120        let (_, y) = cell_xy(90.0, 0.0).unwrap();
121        assert_eq!(y, GRID as u32 - 1);
122    }
123
124    #[test]
125    fn rejects_coordinates_that_are_not_positions() {
126        assert_eq!(check_latitude(90.5), Err(MortonError::LatitudeOutOfRange));
127        assert_eq!(check_latitude(f64::NAN), Err(MortonError::NotFinite));
128        assert_eq!(
129            normalize_longitude(f64::INFINITY),
130            Err(MortonError::NotFinite)
131        );
132    }
133
134    #[test]
135    fn interleave_places_longitude_in_the_even_bits() {
136        assert_eq!(interleave(1, 0), 0b01);
137        assert_eq!(interleave(0, 1), 0b10);
138        assert_eq!(interleave(0xFFFF, 0), 0x5555_5555);
139        assert_eq!(interleave(0, 0xFFFF), 0xAAAA_AAAA);
140    }
141
142    #[test]
143    fn corners_of_a_cell_share_its_key_prefix() {
144        // Two positions in the same maximum-depth cell must produce the same
145        // key, or the cache would answer them from different ranges.
146        let first = key(37.5119, -122.2495).unwrap();
147        let second = key(37.511_900_1, -122.249_500_1).unwrap();
148        assert_eq!(first, second);
149    }
150}
umsh_regiondb/morton.rs

umsh_regiondb/
morton.rs
