irox_carto/
epsg3857.rs

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
// SPDX-License-Identifier: MIT
// Copyright 2025 IROX Contributors
//

//!
//! Spherical (Web) Mercator Map Projection

use core::f64::consts::{PI, TAU};

use irox_units::units::angle::{self, Angle};

use crate::coordinate::{CartesianCoordinate, EllipticalCoordinate, Latitude, Longitude};
use crate::geo::standards::wgs84::WGS84_SHAPE;
use crate::geo::EllipticalShape;
use crate::proj::Projection;

pub const SPHERICAL_MERCATOR_SHAPE: EllipticalShape = EllipticalShape::EpsgDatum(3857);

pub struct SphericalMercatorProjection {
    zoom_level: u8,
}

impl SphericalMercatorProjection {
    #[must_use]
    pub fn new(zoom_level: u8) -> SphericalMercatorProjection {
        SphericalMercatorProjection { zoom_level }
    }

    #[must_use]
    pub fn tile_x_index(&self, coordinate: &EllipticalCoordinate) -> f64 {
        let lon_deg = coordinate.get_longitude().0.as_degrees().value();
        let offset = (lon_deg + 180.) / 360.;
        let max_tile = f64::from(1 << self.zoom_level);
        offset * max_tile
    }

    #[must_use]
    pub fn tile_y_index(&self, coordinate: &EllipticalCoordinate) -> f64 {
        let lat_rad = coordinate.get_latitude().0.as_radians().value();

        let y = lat_rad.tan().asinh();
        let y = (1. - (y / PI)) / 2.;
        let max_tile = f64::from(1 << self.zoom_level);
        max_tile * y
    }

    #[must_use]
    pub fn latitude(&self, tile_y: f64) -> Latitude {
        let offset = 1. - (2. * tile_y) / f64::from(1 << self.zoom_level);
        Latitude(Angle::new_radians((PI * offset).sinh().atan()))
    }

    #[must_use]
    pub fn longitude(&self, tile_x: f64) -> Longitude {
        let offset = tile_x / f64::from(1 << self.zoom_level);
        Longitude(Angle::new_radians(offset * TAU - PI))
    }

    #[must_use]
    pub fn max_tile_index(&self) -> u64 {
        (1 << self.zoom_level) - 1
    }
}

impl Projection for SphericalMercatorProjection {
    fn get_center_coords(&self) -> &EllipticalCoordinate {
        &CENTER_COORDS
    }

    fn project_to_cartesian(&self, coord: &EllipticalCoordinate) -> CartesianCoordinate {
        let x = self.tile_x_index(coord) * TILE_TO_PIXEL;
        let y = self.tile_y_index(coord) * TILE_TO_PIXEL;
        let z = f64::from(self.zoom_level);

        CartesianCoordinate::new_meters(x, y, z)
    }

    fn project_to_elliptical(&self, coord: &CartesianCoordinate) -> EllipticalCoordinate {
        let lat = self.latitude(coord.get_y().as_meters().value());
        let lon = self.longitude(coord.get_x().as_meters().value());

        EllipticalCoordinate::new(lat, lon, WGS84_SHAPE)
    }
}

pub const UPPER_LEFT_COORDINATE_X: f64 = -180.0;
pub const UPPER_LEFT_COORDINATE_Y: f64 = 85.051_128_779_806_59;

pub const LOWER_LEFT_COORDINATE_X: f64 = -180.0;
pub const LOWER_LEFT_COORDINATE_Y: f64 = -85.051_128_779_806_59;

pub const UPPER_RIGHT_COORDINATE_X: f64 = -UPPER_LEFT_COORDINATE_X;
pub const UPPER_RIGHT_COORDINATE_Y: f64 = UPPER_LEFT_COORDINATE_Y;

pub const LOWER_RIGHT_COORDINATE_X: f64 = -LOWER_LEFT_COORDINATE_X;
pub const LOWER_RIGHT_COORDINATE_Y: f64 = LOWER_LEFT_COORDINATE_Y;

pub static CENTER_COORDS: EllipticalCoordinate =
    EllipticalCoordinate::new(Latitude(angle::ZERO), Longitude(angle::ZERO), WGS84_SHAPE);

// pub const BOUNDS: Bounds<CartesianCoordinate> = Bounds::new()

const TILE_TO_PIXEL: f64 = 40.743_665_431_525_21;

#[cfg(test)]
mod test {
    use crate::coordinate::EllipticalCoordinate;

    use super::SphericalMercatorProjection;

    #[test]
    pub fn test1() {
        let sm1 = SphericalMercatorProjection::new(1);

        assert_eq!(0.0, sm1.latitude(1.0).0.as_degrees().value());
        assert_eq!(0.0, sm1.longitude(1.0).0.as_degrees().value());
    }

    #[test]
    pub fn test2() {
        let sm = SphericalMercatorProjection::new(10);

        let coord = EllipticalCoordinate::new_degrees_wgs84(24.846_562, -81.914);

        assert_eq!(439, sm.tile_y_index(&coord) as u64);
        assert_eq!(279, sm.tile_x_index(&coord) as u64);

        let max_tile = 1 << sm.zoom_level;
        assert_eq!(2_u64.pow(sm.zoom_level.into()), max_tile as u64);

        let invy = max_tile - 439 - 1;
        assert_eq!(invy, 584);
    }
}