gistools/geometry/wm/
coords.rs

1use crate::space::EARTH_CIRCUMFERENCE;
2use core::f64::consts::{PI, TAU};
3use libm::{atan, cos, exp, floor, fmax, fmin, log, pow, sin, tan};
4use s2json::{GetXY, NewXY, Point};
5
6/// 900913 (Web Mercator) constant
7pub const A: f64 = 6_378_137.0;
8/// 900913 (Web Mercator) max extent
9pub const MAXEXTENT: f64 = 20_037_508.342789244;
10/// 900913 (Web Mercator) maximum latitude
11pub const MAXLAT: f64 = 85.0511287798;
12
13/// The source of the coordinate inputs
14#[derive(Debug, Clone, PartialEq)]
15pub enum Source {
16    /// The WGS84 projection
17    WGS84,
18    /// The Google (900913) projection
19    Google,
20}
21
22/// Given a zoom and tilesize, build mercator positional attributes
23fn get_zoom_size(zoom: f64, tile_size: f64) -> (f64, f64, f64, f64) {
24    let size = tile_size * pow(2., zoom);
25    (size / 360., size / TAU, size / 2., size)
26}
27
28/// Convert Longitude and Latitude to a mercator pixel coordinate
29/// Return the mercator pixel coordinate
30pub fn ll_to_px<P: GetXY + NewXY>(
31    lonlat: &P,
32    zoom: f64,
33    anti_meridian: Option<bool>,
34    tile_size: Option<u64>,
35) -> P {
36    let anti_meridian = anti_meridian.unwrap_or(false);
37    let tile_size = tile_size.unwrap_or(512) as f64;
38
39    let (bc, cc, zc, ac) = get_zoom_size(zoom, tile_size);
40    let expansion = if anti_meridian { 2. } else { 1. };
41    let d = zc;
42    let f = sin((lonlat.y()).to_radians()).clamp(-0.999999999999, 0.999999999999);
43    let mut x = d + lonlat.x() * bc;
44    let mut y = d + 0.5 * log((1. + f) / (1. - f)) * -cc;
45    if x > ac * expansion {
46        x = ac * expansion;
47    }
48    if y > ac {
49        y = ac;
50    }
51
52    P::new_xy(x, y)
53}
54
55/// Convert mercator pixel coordinates to Longitude and Latitude
56/// Return the Longitude and Latitude
57pub fn px_to_ll<P: GetXY + NewXY>(xy: &P, zoom: f64, tile_size: Option<u64>) -> P {
58    let tile_size = tile_size.unwrap_or(512) as f64;
59    let (bc, cc, zc, _) = get_zoom_size(zoom, tile_size);
60    let g = (xy.y() - zc) / -cc;
61    let lon = (xy.x() - zc) / bc;
62    let lat = (2. * atan(exp(g)) - 0.5 * PI).to_degrees();
63
64    P::new_xy(lon, lat)
65}
66
67/// Convert Longitude and Latitude to a mercator x-y coordinates
68/// Return the mercator x-y coordinates
69pub fn ll_to_merc<P: GetXY + NewXY>(lonlan: &P) -> P {
70    let mut x = (A * lonlan.x()).to_radians();
71    let mut y = A * log(tan(PI * 0.25 + (0.5 * lonlan.y()).to_radians()));
72    // if xy value is beyond maxextent (e.g. poles), return maxextent.
73    x = x.clamp(-MAXEXTENT, MAXEXTENT);
74    y = y.clamp(-MAXEXTENT, MAXEXTENT);
75
76    P::new_xy(x, y)
77}
78
79/// Convert mercator x-y coordinates to Longitude and Latitude
80/// Return the Longitude and Latitude
81pub fn merc_to_ll<P: GetXY + NewXY>(merc: &P) -> P {
82    let x = (merc.x() / A).to_degrees();
83    let y = (0.5 * PI - 2. * atan(exp(-merc.y() / A))).to_degrees();
84    P::new_xy(x, y)
85}
86
87/// Convert a pixel coordinate to a tile x-y coordinate
88/// Return the tile x-y
89pub fn px_to_tile<P: GetXY>(px: &P, tile_size: Option<u64>) -> (i64, i64) {
90    let tile_size = tile_size.unwrap_or(512) as f64;
91    (floor(px.x() / tile_size) as i64, floor(px.y() / tile_size) as i64)
92}
93
94/// Convert a tile x-y-z to a bbox of the form `[w, s, e, n]`
95/// Return the bbox
96pub fn tile_to_bbox(tile: (u8, i64, i64), tile_size: Option<u64>) -> (i64, i64, i64, i64) {
97    let tile_size = tile_size.unwrap_or(512) as i64;
98    let (_zoom, x, y) = tile;
99    let min_x = x * tile_size;
100    let min_y = y * tile_size;
101    let max_x = min_x + tile_size;
102    let max_y = min_y + tile_size;
103
104    (min_x, min_y, max_x, max_y)
105}
106
107/// Convert a lat-lon and zoom to the tile's x-y coordinates
108///
109/// Note: You can just pass in a `&(f64, f64)` for `&P`
110///
111/// ## Returns
112/// The tile x-y
113pub fn ll_to_tile<P: GetXY + NewXY>(lonlat: &P, zoom: f64, tile_size: Option<u64>) -> (i64, i64) {
114    let px = ll_to_px(lonlat, zoom, Some(false), tile_size);
115    px_to_tile(&px, tile_size)
116}
117
118/// given a lon-lat and tile, find the offset in pixels
119/// return the tile xy pixel
120pub fn ll_to_tile_px<P: GetXY + NewXY>(
121    lonlat: &P,
122    tile: (u8, i64, i64),
123    tile_size: Option<u64>,
124) -> P {
125    let (zoom, x, y) = tile;
126    let tile_size = tile_size.unwrap_or(512);
127    let tile_size_f = tile_size as f64;
128    let px = ll_to_px(lonlat, zoom as f64, Some(false), Some(tile_size));
129    let tile_x_start = x as f64 * tile_size_f;
130    let tile_y_start = y as f64 * tile_size_f;
131
132    P::new_xy((px.x() - tile_x_start) / tile_size_f, (px.y() - tile_y_start) / tile_size_f)
133}
134
135/// Convert a bbox of the form `[w, s, e, n]` to a bbox of the form `[w, s, e, n]`
136/// The result can be in lon-lat (WGS84) or WebMercator (900913)
137pub fn convert_bbox(bbox: (f64, f64, f64, f64), source: Source) -> (f64, f64, f64, f64) {
138    let low: Point;
139    let high: Point;
140    match source {
141        Source::WGS84 => {
142            low = merc_to_ll(&Point(bbox.0, bbox.1));
143            high = merc_to_ll(&Point(bbox.2, bbox.3));
144        }
145        Source::Google => {
146            low = ll_to_merc(&Point(bbox.0, bbox.1));
147            high = ll_to_merc(&Point(bbox.2, bbox.3));
148        }
149    };
150    (low.0, low.1, high.0, high.1)
151}
152
153/// Convert a tile x-y-z to a bbox of the form `[w, s, e, n]`
154/// The result can be in lon-lat (WGS84) or WebMercator (900913)
155/// The default result is in WebMercator (900913)
156pub fn xyz_to_bbox(
157    x: i64,
158    y: i64,
159    zoom: f64,
160    tms_style: Option<bool>,
161    source: Option<Source>,
162    tile_size: Option<u64>,
163) -> (f64, f64, f64, f64) {
164    let x = x as f64;
165    let mut y = y as f64;
166    let tms_style = tms_style.unwrap_or(true);
167    let source = source.unwrap_or(Source::Google);
168    let tile_size = tile_size.unwrap_or(512);
169    let tile_size_f = tile_size as f64;
170    // Convert xyz into bbox with srs WGS84
171    // if tmsStyle, the y is inverted
172    if tms_style {
173        y = pow(2., zoom) - 1. - y;
174    }
175    // Use +y to make sure it's a number to avoid inadvertent concatenation.
176    let bl = Point(x * tile_size_f, (y + 1.) * tile_size_f);
177    // Use +x to make sure it's a number to avoid inadvertent concatenation.
178    let tr = Point((x + 1.) * tile_size_f, y * tile_size_f);
179    // to pixel-coordinates
180    let px_bl = px_to_ll(&bl, zoom, Some(tile_size));
181    let px_tr = px_to_ll(&tr, zoom, Some(tile_size));
182
183    match source {
184        Source::Google => {
185            let ll_bl = ll_to_merc(&px_bl);
186            let ll_tr = ll_to_merc(&px_tr);
187            (ll_bl.0, ll_bl.1, ll_tr.0, ll_tr.1)
188        }
189        _ => (px_bl.0, px_bl.1, px_tr.0, px_tr.1),
190    }
191}
192
193/// Convert a bbox of the form `[w, s, e, n]` to a tile's bounding box
194/// in the form of { minX, maxX, minY, maxY }
195/// The bbox can be in lon-lat (WGS84) or WebMercator (900913)
196/// The default expectation is in WebMercator (900913)
197/// returns the tile's bounding box
198pub fn bbox_to_xyz_bounds(
199    bbox: (f64, f64, f64, f64),
200    zoom: f64,
201    tms_style: Option<bool>,
202    source: Option<Source>,
203    tile_size: Option<u64>,
204) -> (i64, i64, i64, i64) {
205    let tms_style = tms_style.unwrap_or(true);
206    let source = source.unwrap_or(Source::WGS84);
207    let tile_size = tile_size.unwrap_or(512);
208    let tile_size_f: f64 = tile_size as f64;
209
210    let mut bl = Point(bbox.0, bbox.1); // bottom left
211    let mut tr = Point(bbox.2, bbox.3); // top right
212
213    if source == Source::Google {
214        bl = ll_to_merc(&bl);
215        tr = ll_to_merc(&tr);
216    }
217    let px_bl = ll_to_px(&bl, zoom, Some(false), Some(tile_size));
218    let px_tr = ll_to_px(&tr, zoom, Some(false), Some(tile_size));
219    let x = (floor(px_bl.0 / tile_size_f), floor((px_tr.0 - 1.0) / tile_size_f));
220    let y = (floor(px_tr.1 / tile_size_f), floor((px_bl.1 - 1.0) / tile_size_f));
221
222    let mut bounds = (fmin(x.0, x.1), fmin(y.0, y.1), fmax(x.0, x.1), fmax(y.0, y.1));
223
224    if tms_style {
225        let zoom_diff = pow(2., zoom) - 1.;
226        bounds.1 = zoom_diff - bounds.3;
227        bounds.3 = zoom_diff - bounds.1;
228    }
229
230    let min_x = fmax(bounds.0, 0.) as i64;
231    let min_y = fmax(bounds.1, 0.) as i64;
232    let max_x = fmin(bounds.2, pow(2., zoom) - 1.) as i64;
233    let max_y = fmin(bounds.3, pow(2., zoom) - 1.) as i64;
234
235    (min_x, min_y, max_x, max_y)
236}
237
238/// The circumference at a line of latitude in meters.
239pub fn circumference_at_latitude(latitude: f64, circumference: Option<f64>) -> f64 {
240    let circumference = circumference.unwrap_or(EARTH_CIRCUMFERENCE);
241    circumference * cos(latitude.to_radians())
242}
243
244/// Convert longitude to mercator projection X-Value
245/// returns the X-Value
246pub fn lng_to_mercator_x(lng: f64) -> f64 {
247    (180.0 + lng) / 360.0
248}
249
250/// Convert latitude to mercator projection Y-Value
251/// returns the Y-Value
252pub fn lat_to_mercator_y(lat: f64) -> f64 {
253    (180. - (180. / PI) * log(tan(PI / 4. + (lat * PI) / 360.))) / 360.
254}
255
256/// Convert altitude to mercator projection Z-Value
257/// returns the Z-Value
258pub fn altitude_to_mercator_z(altitude: f64, lat: f64, circumference: Option<f64>) -> f64 {
259    altitude / circumference_at_latitude(lat, circumference)
260}
261
262/// Convert mercator projection's X-Value to longitude
263/// returns the longitude
264pub fn lng_from_mercator_x(x: f64) -> f64 {
265    x * 360. - 180.
266}
267
268/// Convert mercator projection's Y-Value to latitude
269/// returns the latitude
270pub fn lat_from_mercator_y(y: f64) -> f64 {
271    let y2 = 180. - y * 360.;
272    (360. / PI) * atan(exp((y2 * PI) / 180.)) - 90.
273}
274
275/// Convert mercator projection's Z-Value to altitude
276/// returns the altitude
277pub fn altitude_from_mercator_z(z: f64, y: f64, circumference: Option<f64>) -> f64 {
278    z * circumference_at_latitude(lat_from_mercator_y(y), circumference)
279}
280
281/// Determine the Mercator scale factor for a given latitude, [See more](https://en.wikipedia.org/wiki/Mercator_projection#Scale_factor)
282///
283/// At the equator the scale factor will be 1, which increases at higher latitudes.
284/// returns the scale factor
285pub fn mercator_lat_scale(lat: f64) -> f64 {
286    1. / cos((lat * PI) / 180.)
287}
gistools/geometry/wm/coords.rs

gistools/geometry/wm/
coords.rs
