webmerc/
lib.rs

1use std::f64::consts::{PI, TAU, FRAC_PI_4};
2
3
4#[cfg(test)]
5#[macro_use]
6extern crate approx;
7
8pub const MAX_ZOOM: u64 = 29;
9
10type D = f64;
11type M = f64;
12type P = u64;
13type Z = u64;
14type T = u64;
15
16
17#[derive(Debug)]
18pub struct GlobalMercator {
19    tile_size: u64,
20    a: f64,
21    initial_resolution: f64,
22}
23
24impl GlobalMercator {
25    /// Create a new instance of the coordinate transformer
26    ///
27    /// Example:
28    /// ```
29    ///  let transformer = webmerc::GlobalMercator::new(256);
30    ///  assert_eq!(transformer.tile_size(), 256);
31    /// ```
32    pub fn new(tile_size: u64) -> Self {
33        let a = 6378137.0;
34        let initial_resolution = TAU * a / (tile_size as f64);
35
36        Self {
37            tile_size,
38            initial_resolution,
39            a,
40        }
41    }
42
43    /// Get the tile_size used to instantiate the transformer
44    ///
45    /// Example:
46    /// ```
47    /// let transformer = webmerc::GlobalMercator::new(256);
48    /// assert_eq!(transformer.tile_size(), 256);
49    /// ```
50    pub fn tile_size(&self) -> u64 {
51        self.tile_size
52    }
53
54    /// Converts given lat/lon in WGS84 Datum to XY in Spherical Mercator EPSG:3857
55    ///
56    /// Example:
57    /// ```
58    /// # #[macro_use] extern crate approx;
59    /// let transformer = webmerc::GlobalMercator::new(256);
60    /// let (x, y) = transformer.lat_lon_to_meters(10.0, 10.0);
61    /// assert_relative_eq!(x, 1113194.91, epsilon = 1e-2);
62    /// assert_relative_eq!(y, 1118889.97, epsilon = 1e-2);
63    /// ```
64    pub fn lat_lon_to_meters(&self, lon: D, lat: D) -> (M, M) {
65        let mx = self.a * lon.to_radians();
66        let my = self.a * f64::ln(f64::tan(FRAC_PI_4 + lat.to_radians() / 2.0));
67        (mx, my)
68    }
69
70    /// Converts XY point from Spherical Mercator EPSG:3857 to lat/lon in WGS84 Datum
71    ///
72    /// ```
73    /// # #[macro_use] extern crate approx;
74    /// let transformer = webmerc::GlobalMercator::new(256);
75    /// let (lon, lat) = transformer.meters_2_lat_lon(-20037508.0, -20037508.0);
76    /// assert_relative_eq!(lon, -179.999997, epsilon = 1e-6);
77    /// assert_relative_eq!(lat, -85.051129, epsilon = 1e-6);
78    /// ```
79    pub fn meters_2_lat_lon(&self, mx: M, my: M) -> (D, D) {
80        let lon = (mx / self.a).to_degrees();
81        let lat = (f64::atan(f64::sinh(my / self.a))).to_degrees();
82        (lon, lat)
83    }
84
85    /// Converts pixel coordinates in given Zoom level of pyramid to EPSG:3857
86    pub fn pixels_to_meters(&self, px: P, py: P, zoom: Z) -> (M, M) {
87        let res = self.resolution(zoom);
88        let mx = (px as f64) * res - PI * self.a;
89        let my = (py as f64) * res - PI * self.a;
90        (mx, my)
91    }
92
93    /// Converts EPSG:3857 to pyramid pixel coordinates in given Zoom level
94    pub fn meters_to_pixels(&self, mx: M, my: M, zoom: Z) -> (P, P) {
95        let res = self.resolution(zoom);
96        let px = M::floor((mx + PI * self.a) / res) as P;
97        let py = M::floor((my + PI * self.a) / res) as P;
98        (px, py)
99    }
100
101    /// Returns a Tile covering region in given pixel coordinates
102    pub fn pixels_to_tile(&self, px: P, py: P) -> (T, T) {
103        let tx = f64::ceil(px as f64 / self.tile_size as f64) as u64 - 1;
104        let ty = f64::ceil(py as f64 / self.tile_size as f64) as u64 - 1;
105        (tx, ty)
106    }
107
108    /// Move the origin of pixel coordinates to top-left corner
109    pub fn pixels_to_raster(&self, px: P, py: P, zoom: Z) -> (P, P) {
110        let map_size = (self.tile_size) << zoom;
111        (px, map_size - py)
112    }
113
114    /// Returns tile for given mercator coordinates
115    pub fn meters_to_tile(&self, mx: M, my: M, zoom: Z) -> (T, T) {
116        let (px, py) = self.meters_to_pixels(mx, my, zoom);
117        self.pixels_to_tile(px, py)
118    }
119
120    /// Returns bounds of the given tile in EPSG:3857 coordinates
121    pub fn tile_bounds(&self, tx: T, ty: T, zoom: Z) -> (M, M, M, M) {
122        let (minx, miny) = self.pixels_to_meters(tx * self.tile_size, ty * self.tile_size, zoom);
123        let (maxx, maxy) = self.pixels_to_meters(
124            (tx + 1) * self.tile_size,
125            (ty + 1) * self.tile_size,
126            zoom,
127        );
128        (minx, miny, maxx, maxy)
129    }
130
131    /// Returns bounds of the given tile in latitude/longitude using WGS84 datum
132    pub fn tile_lat_lon_bounds(&self, tx: T, ty: T, zoom: Z) -> (D, D, D, D) {
133        let (min_x, min_y, max_x, max_y) = self.tile_bounds(tx, ty, zoom);
134        let (min_lon, min_lat) = self.meters_2_lat_lon(min_x, min_y);
135        let (max_lon, max_lat) = self.meters_2_lat_lon(max_x, max_y);
136        (min_lon, min_lat, max_lon, max_lat)
137    }
138
139    /// Resolution (meters/pixel) for given zoom level (measured at Equator)
140    pub fn resolution(&self, zoom: Z) -> f64 {
141        self.initial_resolution / (2.0f64.powi(zoom as i32))
142    }
143
144    /// Maximal scaledown zoom of the pyramid closest to the pixelSize.
145    pub fn zoom_for_pixel_size(&self, pixel_size: f64) -> Z {
146        for i in 0..=MAX_ZOOM {
147            if pixel_size > self.resolution(i) {
148                return std::cmp::max(0, i - 1);
149            }
150        }
151        panic!("Invalid pixel_size: {}", pixel_size);
152    }
153
154    /// Converts TMS tile coordinates to Google Tile coordinates and vice versa
155    pub fn google_tile(&self, tx: T, ty: T, zoom: Z) -> (T, T) {
156        (tx, (2u64.pow(zoom as u32) - 1) - ty)
157    }
158
159    /// Converts TMS tile coordinates to Microsoft quad_tree
160    pub fn quad_tree(&self, tx: T, ty: T, zoom: Z) -> String {
161        let mut quad_key = String::new();
162        let ty = (f64::powi(2.0, zoom as i32) - 1.0) as T - ty;
163        for i in (1..(zoom + 1) as i32).rev() {
164            let mut digit = 0;
165            let mask = 1 << (i - 1);
166            if (tx & mask) != 0 {
167                digit += 1;
168            }
169            if (ty & mask) != 0 {
170                digit += 2;
171            }
172            quad_key.push_str(format!("{}", digit).as_str());
173        }
174        quad_key
175    }
176}
177
178// Add tests
179#[cfg(test)]
180mod tests {
181    use super::*;
182
183    #[test]
184    fn test_new() {
185        let transformer = GlobalMercator::new(256);
186        assert_eq!(transformer.tile_size, 256);
187    }
188
189    #[test]
190    fn lat_lon_to_meters() {
191        let transformer = GlobalMercator::new(256);
192        let (x, y) = transformer.lat_lon_to_meters(0.0, 0.0);
193        assert_relative_eq!(x, 0.0, epsilon = 1e-9);
194        assert_relative_eq!(y, 0.0, epsilon = 1e-9);
195
196        let (x, y) = transformer.lat_lon_to_meters(10.0, 10.0);
197        assert_relative_eq!(x, 1113194.91, epsilon = 1e-2);
198        assert_relative_eq!(y, 1118889.97, epsilon = 1e-2);
199
200        // These are the max lat lon for Google maps. More North or South and numerical errors occur
201        let (x, y) = transformer.lat_lon_to_meters(179.999_996_920_671_83, 85.051_128_514_163);
202        assert_relative_eq!(x, 20037508.0, epsilon = 1e-2);
203        assert_relative_eq!(y, 20037508.0, epsilon = 1e-2);
204
205        let (x, y) = transformer.lat_lon_to_meters(180.0, 85.0);
206        assert_relative_eq!(x, 20037508.34, epsilon = 5e-3);
207        assert_relative_eq!(y, 19971868.88, epsilon = 5e-3);
208        //
209        let (x, y) = transformer.lat_lon_to_meters(-180.0, -85.0);
210        assert_relative_eq!(x, -20037508.34, epsilon = 5e-3);
211        assert_relative_eq!(y, -19971868.88, epsilon = 5e-3);
212    }
213
214    #[test]
215    fn meters_to_lat_lon() {
216        let transformer = GlobalMercator::new(256);
217        let (lon, lat) = transformer.meters_2_lat_lon(0.0, 0.0);
218        assert_relative_eq!(lon, 0.0, epsilon = f64::EPSILON);
219        assert_relative_eq!(lat, 0.0, epsilon = f64::EPSILON);
220
221        let (lon, lat) = transformer.meters_2_lat_lon(-20037508.0, -20037508.0);
222        assert_relative_eq!(lon, -179.999997, epsilon = 1e-6);
223        assert_relative_eq!(lat, -85.051129, epsilon = 1e-6);
224
225        let (lon, lat) = transformer.meters_2_lat_lon(20037508.0, 20037508.0);
226        assert_relative_eq!(lon, 179.999997, epsilon = 1e-6);
227        assert_relative_eq!(lat, 85.051129, epsilon = 1e-6);
228    }
229
230    #[test]
231    fn tile_bounds() {
232        let transformer = GlobalMercator::new(256);
233
234        let c = 20_037_508.34;
235
236        for z in 0..MAX_ZOOM {
237            let num_tiles = 2u64.pow(z as u32);
238            let m = c * 2.0 / (num_tiles as f64) - c;
239
240            let bounds = transformer.tile_bounds(0, 0, z);
241            assert_relative_eq!(bounds.0, -c, epsilon = 5e-3);
242            assert_relative_eq!(bounds.1, -c, epsilon = 5e-3);
243            assert_relative_eq!(bounds.2, m, epsilon = 5e-3);
244            assert_relative_eq!(bounds.3, m, epsilon = 5e-3);
245
246            let bounds = transformer.tile_bounds(num_tiles - 1, 0, z);
247            assert_relative_eq!(bounds.0, -m, epsilon = 5e-3);
248            assert_relative_eq!(bounds.1, -c, epsilon = 5e-3);
249            assert_relative_eq!(bounds.2, c, epsilon = 5e-3);
250            assert_relative_eq!(bounds.3, m, epsilon = 5e-3);
251        }
252    }
253}
webmerc/lib.rs

webmerc/
lib.rs
