rustial-engine 0.0.1

//! Web Mercator (EPSG:3857) projection.
//!
//! # Notes
//!
//! Web Mercator is the de-facto projection used by slippy-map tile systems.
//! It is conformal (preserves local angles) but strongly distorts area and
//! scale toward the poles.
//!
//! Valid latitude range is approximately +/-85.051129 degrees. Inputs outside
//! this range produce very large values, `inf`, or `NaN` with the raw
//! projection formula.
//!
//! # Altitude behavior
//!
//! This projection operates on latitude/longitude only. Altitude (`z`) is
//! passed through unchanged in both projection directions.

use crate::bounds::GeoBounds;
use crate::coord::{GeoCoord, WorldCoord, MAX_MERCATOR_LAT};
use crate::ellipsoid::Ellipsoid;
use crate::projection::Projection;
use std::f64::consts::PI;

/// Semi-major axis of the WGS-84 ellipsoid in meters.
///
/// Web Mercator (EPSG:3857) treats the Earth as a sphere with this radius.
const EARTH_RADIUS: f64 = Ellipsoid::WGS84.a;

/// Web Mercator projection utilities.
///
/// `project`/`unproject` are fast, allocation-free wrappers over the
/// [`Projection`] trait implementation.
pub struct WebMercator;

impl WebMercator {
    /// Project a geographic coordinate to Web Mercator world coordinates (meters).
    ///
    /// Input is assumed to already be in valid Web Mercator range.
    #[inline]
    pub fn project(geo: &GeoCoord) -> WorldCoord {
        <Self as Projection>::project(&Self, geo)
    }

    /// Checked variant of [`project`](Self::project).
    ///
    /// Returns `None` when `geo` is outside Web Mercator valid range.
    #[inline]
    pub fn project_checked(geo: &GeoCoord) -> Option<WorldCoord> {
        if !geo.is_web_mercator_valid() {
            return None;
        }
        Some(Self::project(geo))
    }

    /// Project with automatic Mercator clamping/wrapping.
    ///
    /// Latitude is clamped to +/-85.051129 and longitude wrapped to
    /// `[-180, 180]` before projection.
    #[inline]
    pub fn project_clamped(geo: &GeoCoord) -> WorldCoord {
        Self::project(&geo.clamped_mercator())
    }

    /// Inverse-project Web Mercator coordinates back to geographic coordinates.
    #[inline]
    pub fn unproject(world: &WorldCoord) -> GeoCoord {
        <Self as Projection>::unproject(&Self, world)
    }

    /// The maximum extent of Web Mercator along one axis, in meters.
    ///
    /// Equals `R * PI` where `R` is the projection sphere radius.
    #[inline]
    pub fn max_extent() -> f64 {
        EARTH_RADIUS * PI
    }

    /// Full width/height of the Web Mercator world square, in meters.
    #[inline]
    pub fn world_size() -> f64 {
        2.0 * Self::max_extent()
    }
}

impl Projection for WebMercator {
    fn project(&self, geo: &GeoCoord) -> WorldCoord {
        // EPSG:3857 forward equations:
        // x = R * lon_rad
        // y = R * ln(tan(pi/4 + lat_rad/2))
        let x = EARTH_RADIUS * geo.lon.to_radians();
        let lat_rad = geo.lat.to_radians();
        let y = EARTH_RADIUS * ((PI / 4.0 + lat_rad / 2.0).tan()).ln();
        WorldCoord::new(x, y, geo.alt)
    }

    fn unproject(&self, world: &WorldCoord) -> GeoCoord {
        // EPSG:3857 inverse equations:
        // lon = x / R
        // lat = 2 * atan(exp(y / R)) - pi/2
        let mut lon = (world.position.x / EARTH_RADIUS).to_degrees();
        lon = ((lon + 180.0).rem_euclid(360.0)) - 180.0;
        let lat = (2.0 * (world.position.y / EARTH_RADIUS).exp().atan() - PI / 2.0).to_degrees();
        GeoCoord::new(lat, lon, world.position.z)
    }

    /// Scale factor for Web Mercator: `sec(lat) = 1 / cos(lat)`.
    ///
    /// Approaches infinity near the poles, which is why the projection
    /// is limited to approximately 85.06 degrees latitude.
    fn scale_factor(&self, geo: &GeoCoord) -> f64 {
        1.0 / geo.lat.to_radians().cos()
    }

    /// Web Mercator is valid within approximately 85.06 degrees latitude
    /// and 180 degrees longitude.
    fn projection_bounds(&self) -> GeoBounds {
        GeoBounds::new(
            GeoCoord::from_lat_lon(-MAX_MERCATOR_LAT, -180.0),
            GeoCoord::from_lat_lon(MAX_MERCATOR_LAT, 180.0),
        )
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn roundtrip_origin() {
        let geo = GeoCoord::from_lat_lon(0.0, 0.0);
        let world = WebMercator::project(&geo);
        let back = WebMercator::unproject(&world);
        assert!((back.lat - geo.lat).abs() < 1e-10);
        assert!((back.lon - geo.lon).abs() < 1e-10);
    }

    #[test]
    fn roundtrip_nonzero() {
        let geo = GeoCoord::from_lat_lon(51.09916, 17.03664);
        let world = WebMercator::project(&geo);
        let back = WebMercator::unproject(&world);
        assert!((back.lat - geo.lat).abs() < 1e-8);
        assert!((back.lon - geo.lon).abs() < 1e-8);
    }

    #[test]
    fn project_checked_rejects_invalid_lat() {
        let geo = GeoCoord::from_lat_lon(89.0, 0.0);
        assert!(WebMercator::project_checked(&geo).is_none());
    }

    #[test]
    fn project_clamped_accepts_invalid_lat() {
        let geo = GeoCoord::from_lat_lon(89.0, 0.0);
        let world = WebMercator::project_clamped(&geo);
        assert!(world.position.y.is_finite());
    }

    #[test]
    fn world_size_is_double_extent() {
        assert!((WebMercator::world_size() - 2.0 * WebMercator::max_extent()).abs() < 1e-10);
    }

    #[test]
    fn scale_factor_equator() {
        let sf = WebMercator.scale_factor(&GeoCoord::from_lat_lon(0.0, 0.0));
        assert!((sf - 1.0).abs() < 1e-10);
    }

    #[test]
    fn scale_factor_60_degrees() {
        // sec(60 deg) = 2.0
        let sf = WebMercator.scale_factor(&GeoCoord::from_lat_lon(60.0, 0.0));
        assert!((sf - 2.0).abs() < 1e-10);
    }

    #[test]
    fn projection_bounds_mercator() {
        let bounds = WebMercator.projection_bounds();
        assert!((bounds.sw().lat - (-MAX_MERCATOR_LAT)).abs() < 1e-10);
        assert!((bounds.ne().lat - MAX_MERCATOR_LAT).abs() < 1e-10);
    }

    #[test]
    fn project_clamped_wraps_longitude() {
        let a = GeoCoord {
            lat: 0.0,
            lon: 190.0,
            alt: 0.0,
        };
        let b = GeoCoord::from_lat_lon(0.0, -170.0);
        let wa = WebMercator::project_clamped(&a);
        let wb = WebMercator::project_clamped(&b);
        assert!((wa.position.x - wb.position.x).abs() < 1e-10);
    }

    #[test]
    fn altitude_passthrough() {
        let geo = GeoCoord::new(45.0, 12.0, 1234.5);
        let world = WebMercator::project(&geo);
        let back = WebMercator::unproject(&world);
        assert!((world.position.z - 1234.5).abs() < 1e-12);
        assert!((back.alt - 1234.5).abs() < 1e-12);
    }

    #[test]
    fn unproject_wraps_longitude_back_into_valid_range() {
        let extent = WebMercator::max_extent();
        let world = WorldCoord::new(extent + 1000.0, 0.0, 0.0);
        let geo = WebMercator::unproject(&world);
        assert!(geo.lon >= -180.0 && geo.lon <= 180.0);
    }
}