rustial_engine/math/
equirectangular.rs

1//! Equirectangular (Plate Carree) projection (EPSG:4326 scaled to meters).
2//!
3//! # Overview
4//!
5//! The Equirectangular projection is the simplest cylindrical projection:
6//! longitude and latitude are mapped linearly to X and Y, then scaled by
7//! the WGS-84 semi-major axis so that output units are **meters** at the
8//! equator.
9//!
10//! ```text
11//! x = R * lon_rad
12//! y = R * lat_rad
13//! ```
14//!
15//! # Distortion
16//!
17//! - **Equator:** no distortion (scale factor = 1.0).
18//! - **Away from equator:** X distances are progressively stretched by
19//!   `sec(lat)`.  Y distances remain undistorted at all latitudes.
20//! - **Poles:** the projection is valid up to +/-90 degrees latitude,
21//!   but X distortion is infinite at the poles.
22//!
23//! # When to use
24//!
25//! Equirectangular is useful as a simple baseline, for data interchange
26//! (many raster datasets are stored in EPSG:4326 pixel grids), and for
27//! server-side tile pre-processing where rendering fidelity is not
28//! critical.  For map display, prefer [`WebMercator`](crate::WebMercator).
29//!
30//! # Altitude
31//!
32//! Altitude (`alt` / `z`) is passed through unchanged in both directions.
33
34use crate::coord::{GeoCoord, WorldCoord};
35use crate::ellipsoid::Ellipsoid;
36use crate::projection::Projection;
37
38/// Equirectangular projection: trivial linear mapping of lon/lat to meters.
39///
40/// Uses the WGS-84 semi-major axis (`R = 6 378 137 m`) as the scaling
41/// factor so that output units are meters at the equator.  This is a
42/// zero-size type -- construction is free.
43pub struct Equirectangular;
44
45impl Projection for Equirectangular {
46    /// Forward projection: geographic degrees to world meters.
47    ///
48    /// `x = R * lon_rad`, `y = R * lat_rad`, `z = alt` (passthrough).
49    fn project(&self, geo: &GeoCoord) -> WorldCoord {
50        let a = Ellipsoid::WGS84.a;
51        WorldCoord::new(a * geo.lon.to_radians(), a * geo.lat.to_radians(), geo.alt)
52    }
53
54    /// Inverse projection: world meters back to geographic degrees.
55    ///
56    /// `lat = degrees(y / R)`, `lon = degrees(x / R)`, `alt = z`.
57    fn unproject(&self, world: &WorldCoord) -> GeoCoord {
58        let a = Ellipsoid::WGS84.a;
59        let lat = (world.position.y / a).to_degrees().clamp(-90.0, 90.0);
60        let lon = ((world.position.x / a).to_degrees() + 180.0).rem_euclid(360.0) - 180.0;
61        GeoCoord::new(
62            lat,
63            lon,
64            world.position.z,
65        )
66    }
67
68    /// Scale factor for Equirectangular: `sec(lat) = 1 / cos(lat)`.
69    ///
70    /// The projection maps longitude linearly, so X distances are
71    /// stretched relative to true distances by `sec(lat)`, the same
72    /// as the meridian convergence factor.  Along Y the scale is 1.0,
73    /// so this returns the *maximum* linear distortion at the point.
74    fn scale_factor(&self, geo: &GeoCoord) -> f64 {
75        1.0 / geo.lat.to_radians().cos()
76    }
77
78    // `projection_bounds` is not overridden: the default full-globe
79    // range (+/-90 lat, +/-180 lon) is correct for Equirectangular.
80}
81
82#[cfg(test)]
83mod tests {
84    use super::*;
85
86    // -- Roundtrip ---------------------------------------------------------
87
88    #[test]
89    fn roundtrip_origin() {
90        let geo = GeoCoord::from_lat_lon(0.0, 0.0);
91        let world = Equirectangular.project(&geo);
92        let back = Equirectangular.unproject(&world);
93        assert!((back.lat - geo.lat).abs() < 1e-10);
94        assert!((back.lon - geo.lon).abs() < 1e-10);
95    }
96
97    #[test]
98    fn roundtrip_nonzero() {
99        let geo = GeoCoord::from_lat_lon(45.0, 90.0);
100        let world = Equirectangular.project(&geo);
101        let back = Equirectangular.unproject(&world);
102        assert!((back.lat - geo.lat).abs() < 1e-8);
103        assert!((back.lon - geo.lon).abs() < 1e-8);
104    }
105
106    #[test]
107    fn poles_roundtrip() {
108        let north = GeoCoord::from_lat_lon(90.0, 0.0);
109        let north_w = Equirectangular.project(&north);
110        let north_back = Equirectangular.unproject(&north_w);
111        assert!((north_back.lat - 90.0).abs() < 1e-9);
112
113        let south = GeoCoord::from_lat_lon(-90.0, 0.0);
114        let south_w = Equirectangular.project(&south);
115        let south_back = Equirectangular.unproject(&south_w);
116        assert!((south_back.lat + 90.0).abs() < 1e-9);
117    }
118
119    #[test]
120    fn anti_meridian_roundtrip() {
121        // +180 and -180 represent the same meridian; the normalisation
122        // in `unproject` canonicalises both to the (-180, 180] range.
123        let east = GeoCoord::from_lat_lon(10.0, 180.0);
124        let east_back = Equirectangular.unproject(&Equirectangular.project(&east));
125        assert!((east_back.lon.abs() - 180.0).abs() < 1e-9);
126
127        let west = GeoCoord::from_lat_lon(10.0, -180.0);
128        let west_back = Equirectangular.unproject(&Equirectangular.project(&west));
129        assert!((west_back.lon.abs() - 180.0).abs() < 1e-9);
130    }
131
132    #[test]
133    fn extreme_latitudes_stable() {
134        for lat in [89.999, -89.999, 89.9, -89.9] {
135            let geo = GeoCoord::from_lat_lon(lat, 45.0);
136            let world = Equirectangular.project(&geo);
137            let back = Equirectangular.unproject(&world);
138            assert!((back.lat - lat).abs() < 1e-6, "lat roundtrip failed for {lat}");
139        }
140    }
141
142    // -- Metric output at the equator -------------------------------------
143
144    #[test]
145    fn equator_scale() {
146        let a = GeoCoord::from_lat_lon(0.0, 0.0);
147        let b = GeoCoord::from_lat_lon(0.0, 1.0);
148        let wa = Equirectangular.project(&a);
149        let wb = Equirectangular.project(&b);
150        let dx = wb.position.x - wa.position.x;
151        // 1 degree of longitude at the equator ~ 111,319.49 m.
152        assert!((dx - 111_319.49).abs() < 1.0);
153    }
154
155    // -- Scale factor ------------------------------------------------------
156
157    #[test]
158    fn scale_factor_equator() {
159        let sf = Equirectangular.scale_factor(&GeoCoord::from_lat_lon(0.0, 0.0));
160        assert!((sf - 1.0).abs() < 1e-10);
161    }
162
163    #[test]
164    fn scale_factor_45_degrees() {
165        // sec(45 deg) = sqrt(2)
166        let sf = Equirectangular.scale_factor(&GeoCoord::from_lat_lon(45.0, 0.0));
167        assert!((sf - std::f64::consts::SQRT_2).abs() < 1e-10);
168    }
169
170    // -- Altitude passthrough ----------------------------------------------
171
172    #[test]
173    fn altitude_passthrough() {
174        let geo = GeoCoord::new(30.0, 60.0, 1234.5);
175        let world = Equirectangular.project(&geo);
176        let back = Equirectangular.unproject(&world);
177        assert!((world.position.z - 1234.5).abs() < 1e-12);
178        assert!((back.alt - 1234.5).abs() < 1e-12);
179    }
180
181    // -- Projection bounds default to full globe --------------------------
182
183    #[test]
184    fn projection_bounds_full_globe() {
185        let bounds = Equirectangular.projection_bounds();
186        assert!((bounds.sw().lat - (-90.0)).abs() < 1e-10);
187        assert!((bounds.ne().lat - 90.0).abs() < 1e-10);
188        assert!((bounds.sw().lon - (-180.0)).abs() < 1e-10);
189        assert!((bounds.ne().lon - 180.0).abs() < 1e-10);
190    }
191}
rustial_engine/math/equirectangular.rs

rustial_engine/math/
equirectangular.rs
