rustial_engine/math/
mercator.rs

1//! Web Mercator (EPSG:3857) projection.
2//!
3//! # Notes
4//!
5//! Web Mercator is the de-facto projection used by slippy-map tile systems.
6//! It is conformal (preserves local angles) but strongly distorts area and
7//! scale toward the poles.
8//!
9//! Valid latitude range is approximately +/-85.051129 degrees. Inputs outside
10//! this range produce very large values, `inf`, or `NaN` with the raw
11//! projection formula.
12//!
13//! # Altitude behavior
14//!
15//! This projection operates on latitude/longitude only. Altitude (`z`) is
16//! passed through unchanged in both projection directions.
17
18use crate::bounds::GeoBounds;
19use crate::coord::{GeoCoord, WorldCoord, MAX_MERCATOR_LAT};
20use crate::ellipsoid::Ellipsoid;
21use crate::projection::Projection;
22use std::f64::consts::PI;
23
24/// Semi-major axis of the WGS-84 ellipsoid in meters.
25///
26/// Web Mercator (EPSG:3857) treats the Earth as a sphere with this radius.
27const EARTH_RADIUS: f64 = Ellipsoid::WGS84.a;
28
29/// Web Mercator projection utilities.
30///
31/// `project`/`unproject` are fast, allocation-free wrappers over the
32/// [`Projection`] trait implementation.
33pub struct WebMercator;
34
35impl WebMercator {
36    /// Project a geographic coordinate to Web Mercator world coordinates (meters).
37    ///
38    /// Input is assumed to already be in valid Web Mercator range.
39    #[inline]
40    pub fn project(geo: &GeoCoord) -> WorldCoord {
41        <Self as Projection>::project(&Self, geo)
42    }
43
44    /// Checked variant of [`project`](Self::project).
45    ///
46    /// Returns `None` when `geo` is outside Web Mercator valid range.
47    #[inline]
48    pub fn project_checked(geo: &GeoCoord) -> Option<WorldCoord> {
49        if !geo.is_web_mercator_valid() {
50            return None;
51        }
52        Some(Self::project(geo))
53    }
54
55    /// Project with automatic Mercator clamping/wrapping.
56    ///
57    /// Latitude is clamped to +/-85.051129 and longitude wrapped to
58    /// `[-180, 180]` before projection.
59    #[inline]
60    pub fn project_clamped(geo: &GeoCoord) -> WorldCoord {
61        Self::project(&geo.clamped_mercator())
62    }
63
64    /// Inverse-project Web Mercator coordinates back to geographic coordinates.
65    #[inline]
66    pub fn unproject(world: &WorldCoord) -> GeoCoord {
67        <Self as Projection>::unproject(&Self, world)
68    }
69
70    /// The maximum extent of Web Mercator along one axis, in meters.
71    ///
72    /// Equals `R * PI` where `R` is the projection sphere radius.
73    #[inline]
74    pub fn max_extent() -> f64 {
75        EARTH_RADIUS * PI
76    }
77
78    /// Full width/height of the Web Mercator world square, in meters.
79    #[inline]
80    pub fn world_size() -> f64 {
81        2.0 * Self::max_extent()
82    }
83}
84
85impl Projection for WebMercator {
86    fn project(&self, geo: &GeoCoord) -> WorldCoord {
87        // EPSG:3857 forward equations:
88        // x = R * lon_rad
89        // y = R * ln(tan(pi/4 + lat_rad/2))
90        let x = EARTH_RADIUS * geo.lon.to_radians();
91        let lat_rad = geo.lat.to_radians();
92        let y = EARTH_RADIUS * ((PI / 4.0 + lat_rad / 2.0).tan()).ln();
93        WorldCoord::new(x, y, geo.alt)
94    }
95
96    fn unproject(&self, world: &WorldCoord) -> GeoCoord {
97        // EPSG:3857 inverse equations:
98        // lon = x / R
99        // lat = 2 * atan(exp(y / R)) - pi/2
100        let mut lon = (world.position.x / EARTH_RADIUS).to_degrees();
101        lon = ((lon + 180.0).rem_euclid(360.0)) - 180.0;
102        let lat = (2.0 * (world.position.y / EARTH_RADIUS).exp().atan() - PI / 2.0).to_degrees();
103        GeoCoord::new(lat, lon, world.position.z)
104    }
105
106    /// Scale factor for Web Mercator: `sec(lat) = 1 / cos(lat)`.
107    ///
108    /// Approaches infinity near the poles, which is why the projection
109    /// is limited to approximately 85.06 degrees latitude.
110    fn scale_factor(&self, geo: &GeoCoord) -> f64 {
111        1.0 / geo.lat.to_radians().cos()
112    }
113
114    /// Web Mercator is valid within approximately 85.06 degrees latitude
115    /// and 180 degrees longitude.
116    fn projection_bounds(&self) -> GeoBounds {
117        GeoBounds::new(
118            GeoCoord::from_lat_lon(-MAX_MERCATOR_LAT, -180.0),
119            GeoCoord::from_lat_lon(MAX_MERCATOR_LAT, 180.0),
120        )
121    }
122}
123
124#[cfg(test)]
125mod tests {
126    use super::*;
127
128    #[test]
129    fn roundtrip_origin() {
130        let geo = GeoCoord::from_lat_lon(0.0, 0.0);
131        let world = WebMercator::project(&geo);
132        let back = WebMercator::unproject(&world);
133        assert!((back.lat - geo.lat).abs() < 1e-10);
134        assert!((back.lon - geo.lon).abs() < 1e-10);
135    }
136
137    #[test]
138    fn roundtrip_nonzero() {
139        let geo = GeoCoord::from_lat_lon(51.09916, 17.03664);
140        let world = WebMercator::project(&geo);
141        let back = WebMercator::unproject(&world);
142        assert!((back.lat - geo.lat).abs() < 1e-8);
143        assert!((back.lon - geo.lon).abs() < 1e-8);
144    }
145
146    #[test]
147    fn project_checked_rejects_invalid_lat() {
148        let geo = GeoCoord::from_lat_lon(89.0, 0.0);
149        assert!(WebMercator::project_checked(&geo).is_none());
150    }
151
152    #[test]
153    fn project_clamped_accepts_invalid_lat() {
154        let geo = GeoCoord::from_lat_lon(89.0, 0.0);
155        let world = WebMercator::project_clamped(&geo);
156        assert!(world.position.y.is_finite());
157    }
158
159    #[test]
160    fn world_size_is_double_extent() {
161        assert!((WebMercator::world_size() - 2.0 * WebMercator::max_extent()).abs() < 1e-10);
162    }
163
164    #[test]
165    fn scale_factor_equator() {
166        let sf = WebMercator.scale_factor(&GeoCoord::from_lat_lon(0.0, 0.0));
167        assert!((sf - 1.0).abs() < 1e-10);
168    }
169
170    #[test]
171    fn scale_factor_60_degrees() {
172        // sec(60 deg) = 2.0
173        let sf = WebMercator.scale_factor(&GeoCoord::from_lat_lon(60.0, 0.0));
174        assert!((sf - 2.0).abs() < 1e-10);
175    }
176
177    #[test]
178    fn projection_bounds_mercator() {
179        let bounds = WebMercator.projection_bounds();
180        assert!((bounds.sw().lat - (-MAX_MERCATOR_LAT)).abs() < 1e-10);
181        assert!((bounds.ne().lat - MAX_MERCATOR_LAT).abs() < 1e-10);
182    }
183
184    #[test]
185    fn project_clamped_wraps_longitude() {
186        let a = GeoCoord {
187            lat: 0.0,
188            lon: 190.0,
189            alt: 0.0,
190        };
191        let b = GeoCoord::from_lat_lon(0.0, -170.0);
192        let wa = WebMercator::project_clamped(&a);
193        let wb = WebMercator::project_clamped(&b);
194        assert!((wa.position.x - wb.position.x).abs() < 1e-10);
195    }
196
197    #[test]
198    fn altitude_passthrough() {
199        let geo = GeoCoord::new(45.0, 12.0, 1234.5);
200        let world = WebMercator::project(&geo);
201        let back = WebMercator::unproject(&world);
202        assert!((world.position.z - 1234.5).abs() < 1e-12);
203        assert!((back.alt - 1234.5).abs() < 1e-12);
204    }
205
206    #[test]
207    fn unproject_wraps_longitude_back_into_valid_range() {
208        let extent = WebMercator::max_extent();
209        let world = WorldCoord::new(extent + 1000.0, 0.0, 0.0);
210        let geo = WebMercator::unproject(&world);
211        assert!(geo.lon >= -180.0 && geo.lon <= 180.0);
212    }
213}
rustial_engine/math/mercator.rs

rustial_engine/math/
mercator.rs
