gistools/geometry/wm/
coords.rs

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

gistools/geometry/wm/
coords.rs
