ezu_core/
coord.rs

1//! Web-Mercator tile / world coordinate utilities.
2//!
3//! World coordinates are in the unit square `[0, 1] x [0, 1]` covering the whole
4//! Web-Mercator projection. This avoids zoom-dependent units and makes
5//! deterministic seeding zoom-stable.
6
7/// A Web-Mercator XYZ tile identifier.
8#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
9pub struct TileId {
10    pub z: u8,
11    pub x: u32,
12    pub y: u32,
13}
14
15impl TileId {
16    pub const fn new(z: u8, x: u32, y: u32) -> Self {
17        Self { z, x, y }
18    }
19
20    /// Number of tiles along one axis at this zoom level.
21    #[inline]
22    pub fn axis_tiles(self) -> u32 {
23        1u32 << self.z
24    }
25
26    /// The tile one zoom level up that contains this tile, or `None`
27    /// at zoom 0 (which has no parent).
28    #[inline]
29    pub fn parent(self) -> Option<TileId> {
30        if self.z == 0 {
31            None
32        } else {
33            Some(TileId::new(self.z - 1, self.x >> 1, self.y >> 1))
34        }
35    }
36
37    /// The ancestor tile at the given (lower) zoom, or `None` if `z`
38    /// is greater than or equal to this tile's zoom.
39    #[inline]
40    pub fn ancestor_at(self, z: u8) -> Option<TileId> {
41        if z >= self.z {
42            return None;
43        }
44        let dz = self.z - z;
45        Some(TileId::new(z, self.x >> dz, self.y >> dz))
46    }
47
48    /// `true` iff `other` lies inside this tile's spatial bounds at a
49    /// deeper zoom. A tile is *not* considered its own ancestor.
50    pub fn is_ancestor_of(self, other: TileId) -> bool {
51        if other.z <= self.z {
52            return false;
53        }
54        let dz = other.z - self.z;
55        other.x >> dz == self.x && other.y >> dz == self.y
56    }
57}
58
59/// A position in the global Web-Mercator unit square `[0, 1] x [0, 1]`.
60#[derive(Debug, Clone, Copy, PartialEq)]
61pub struct WorldPos {
62    pub x: f64,
63    pub y: f64,
64}
65
66impl WorldPos {
67    pub const fn new(x: f64, y: f64) -> Self {
68        Self { x, y }
69    }
70}
71
72/// Convert a tile-local position (in `[0, extent]`) to world unit-square coordinates.
73#[inline]
74pub fn tile_to_world(tile: TileId, tx: f64, ty: f64, extent: f64) -> WorldPos {
75    let n = tile.axis_tiles() as f64;
76    WorldPos {
77        x: (tile.x as f64 + tx / extent) / n,
78        y: (tile.y as f64 + ty / extent) / n,
79    }
80}
81
82/// Length of the equator on the WGS-84 ellipsoid, in metres. One world
83/// unit of x spans this at the equator.
84pub const EARTH_CIRCUMFERENCE_M: f64 = 40_075_016.685_578_5;
85
86/// The latitude Web Mercator runs out at, in degrees. The projection
87/// sends the poles to infinity, so the world square stops here.
88pub const MERCATOR_MAX_LAT: f64 = 85.051_128_779_8;
89
90/// World x of a longitude in degrees. `-180` is `0.0`, `180` is `1.0`.
91#[inline]
92pub fn lon_to_world_x(lon_deg: f64) -> f64 {
93    (lon_deg + 180.0) / 360.0
94}
95
96/// World y of a latitude in degrees, clamped to the Mercator domain.
97/// North is `0.0`, south is `1.0`, matching the y-down tile grid.
98#[inline]
99pub fn lat_to_world_y(lat_deg: f64) -> f64 {
100    let lat = lat_deg
101        .clamp(-MERCATOR_MAX_LAT, MERCATOR_MAX_LAT)
102        .to_radians();
103    (1.0 - lat.tan().asinh() / std::f64::consts::PI) / 2.0
104}
105
106/// Longitude in degrees of a world x. Inverse of [`lon_to_world_x`].
107#[inline]
108pub fn world_x_to_lon(wx: f64) -> f64 {
109    wx * 360.0 - 180.0
110}
111
112/// Latitude in degrees of a world y. Inverse of [`lat_to_world_y`].
113#[inline]
114pub fn world_y_to_lat(wy: f64) -> f64 {
115    (std::f64::consts::PI * (1.0 - 2.0 * wy))
116        .sinh()
117        .atan()
118        .to_degrees()
119}
120
121/// Ground metres per world unit at world y `wy`.
122///
123/// Web Mercator inflates distances away from the equator by `1 / cos(lat)`,
124/// and is conformal, so the same factor applies along both axes: a world
125/// unit square at `wy` covers `metres_per_world_unit(wy).powi(2)` square
126/// metres of ground. Callers converting a real-world density (people per
127/// km², say) into one expressed in world or tile-pixel units need this.
128#[inline]
129pub fn metres_per_world_unit(wy: f64) -> f64 {
130    // lat = atan(sinh(pi (1 - 2 wy))), and cos(atan(sinh(u))) = 1/cosh(u),
131    // so the cosine falls out without the round trip through a latitude.
132    let u = std::f64::consts::PI * (1.0 - 2.0 * wy);
133    EARTH_CIRCUMFERENCE_M / u.cosh()
134}
135
136#[cfg(test)]
137mod tests {
138    use super::*;
139
140    #[test]
141    fn parent_walks_one_level() {
142        assert_eq!(TileId::new(0, 0, 0).parent(), None);
143        assert_eq!(TileId::new(3, 5, 6).parent(), Some(TileId::new(2, 2, 3)));
144    }
145
146    #[test]
147    fn ancestor_at_handles_invalid() {
148        let t = TileId::new(5, 10, 20);
149        assert_eq!(t.ancestor_at(5), None); // same zoom → not an ancestor
150        assert_eq!(t.ancestor_at(6), None); // deeper zoom
151        assert_eq!(t.ancestor_at(3), Some(TileId::new(3, 2, 5)));
152        assert_eq!(t.ancestor_at(0), Some(TileId::new(0, 0, 0)));
153    }
154
155    #[test]
156    fn is_ancestor_of() {
157        let parent = TileId::new(5, 10, 20);
158        assert!(parent.is_ancestor_of(TileId::new(7, 41, 81))); // inside
159        assert!(!parent.is_ancestor_of(TileId::new(7, 44, 81))); // wrong x branch
160        assert!(!parent.is_ancestor_of(parent)); // not self
161        assert!(!parent.is_ancestor_of(TileId::new(4, 5, 10))); // ancestor, not descendant
162    }
163
164    #[test]
165    fn lon_lat_round_trip_through_world_coords() {
166        for lon in [-180.0, -74.0, 0.0, 139.7, 180.0] {
167            let wx = lon_to_world_x(lon);
168            assert!((world_x_to_lon(wx) - lon).abs() < 1e-9, "lon {lon}");
169        }
170        for lat in [-80.0, -35.7, 0.0, 51.5, 80.0] {
171            let wy = lat_to_world_y(lat);
172            assert!((world_y_to_lat(wy) - lat).abs() < 1e-9, "lat {lat}");
173        }
174        // North is y = 0 and the equator is the middle of the square.
175        assert!(lat_to_world_y(80.0) < lat_to_world_y(-80.0));
176        assert!((lat_to_world_y(0.0) - 0.5).abs() < 1e-12);
177        // Beyond the projection's domain the value saturates.
178        assert_eq!(lat_to_world_y(89.0), lat_to_world_y(MERCATOR_MAX_LAT));
179    }
180
181    #[test]
182    fn metres_per_world_unit_matches_mercator_scale() {
183        // The equator is the unscaled reference.
184        assert!((metres_per_world_unit(0.5) - EARTH_CIRCUMFERENCE_M).abs() < 1.0);
185        // 60°N sits at wy where cos(lat) = 1/2, so the scale halves.
186        let wy_60n = (1.0 - 60f64.to_radians().tan().asinh() / std::f64::consts::PI) / 2.0;
187        let ratio = metres_per_world_unit(wy_60n) / EARTH_CIRCUMFERENCE_M;
188        assert!((ratio - 0.5).abs() < 1e-9, "ratio {ratio}");
189        // Symmetric about the equator.
190        assert!((metres_per_world_unit(0.2) - metres_per_world_unit(0.8)).abs() < 1e-6);
191    }
192}
ezu_core/coord.rs

ezu_core/
coord.rs
