utiles_core/
fns.rs

1//! Core util(e)ity functions for working with web mercator tiles, bounding boxes, et al
2#![deny(clippy::missing_const_for_fn)]
3
4use std::f64::consts::PI;
5use std::num::FpCategory;
6
7use crate::Point2d;
8use crate::bbox::{BBox, WebBBox};
9use crate::constants::{DEG2RAD, EARTH_CIRCUMFERENCE, EARTH_RADIUS, LL_EPSILON};
10use crate::errors::UtilesCoreResult;
11use crate::point2d;
12use crate::sibling_relationship::SiblingRelationship;
13use crate::tile_zbox::{TileZBox, TileZBoxes};
14use crate::utile;
15use crate::zoom::ZoomOrZooms;
16use crate::{LngLat, Tile, UtilesCoreError};
17
18/// Return upper left lnglat as tuple `(f64, f64)` from x,y,z
19#[must_use]
20pub fn ult(x: u32, y: u32, z: u8) -> (f64, f64) {
21    let z2 = f64::from(2_u32.pow(u32::from(z)));
22    let lon_deg = (f64::from(x) / z2).mul_add(360.0, -180.0);
23    let lat_rad = ((1.0 - 2.0 * f64::from(y) / z2) * PI).sinh().atan();
24    (lon_deg, lat_rad.to_degrees())
25}
26
27/// Return upper left lnglat as `LngLat` from x,y,z
28#[must_use]
29pub fn ul(x: u32, y: u32, z: u8) -> LngLat {
30    let (lon_deg, lat_deg) = ult(x, y, z);
31    LngLat::new(lon_deg, lat_deg)
32}
33
34/// Return lower left lnglat as `LngLat` from x,y,z
35#[must_use]
36pub fn ll(x: u32, y: u32, z: u8) -> LngLat {
37    ul(x, y + 1, z)
38}
39
40/// Return upper right lnglat as `LngLat` from x,y,z
41#[must_use]
42pub fn ur(x: u32, y: u32, z: u8) -> LngLat {
43    ul(x + 1, y, z)
44}
45
46/// Return lower right lnglat as `LngLat` from x,y,z
47#[must_use]
48pub fn lr(x: u32, y: u32, z: u8) -> LngLat {
49    ul(x + 1, y + 1, z)
50}
51
52/// Return tuple (min, max) x/y for a zoom-level
53#[must_use]
54pub fn minmax(zoom: u8) -> (u32, u32) {
55    (0, 2_u32.pow(u32::from(zoom)) - 1)
56}
57
58/// Return true if x, y, z is a valid tile coordinate
59#[must_use]
60pub fn valid(x: u32, y: u32, z: u8) -> bool {
61    let (minx, maxx) = minmax(z);
62    let (miny, maxy) = minmax(z);
63    x >= minx && x <= maxx && y >= miny && y <= maxy
64}
65
66/// Return the inverted/y-flipped y coordinate for a given y and z
67#[must_use]
68#[inline]
69pub fn flipy(y: u32, z: u8) -> u32 {
70    2_u32.pow(u32::from(z)) - 1 - y
71}
72
73/// Return the y-flipped y coordinate for a given y and z
74#[must_use]
75#[inline]
76pub fn yflip(y: u32, z: u8) -> u32 {
77    flipy(y, z)
78}
79
80/// Return cumulative base-tile-id for a tile and the zoom level of the tile
81///
82/// Base-tile-id is the sum of all tiles in all zoom levels below the zoom level
83/// of the tile.
84#[must_use]
85#[inline]
86pub const fn int_2_offset_zoom(i: u64) -> (u64, u8) {
87    if i == 0 {
88        return (0, 0);
89    }
90    let mut acc: u64 = 0;
91    let mut z: u8 = 0;
92    loop {
93        let num_tiles: u64 = (1 << z) * (1 << z);
94        if acc + num_tiles > i {
95            return (i - acc, z);
96        }
97        acc += num_tiles;
98        z += 1;
99    }
100}
101
102/// Calculate the row-major-id for a tile which is the
103/// index of the tile for the zoom level PLUS the total number of tiles
104/// in all zoom levels below it for the zoom level.
105///
106/// (x, y, z) => x + y * 2^z + 1 + 2^z * 2^z
107/// (0,0,0) is 0
108/// (0,0,1) is 1
109/// (0,1,1) is 2
110///
111/// # Examples
112/// ```
113/// use utiles_core::xyz2rmid;
114/// let zzz = xyz2rmid(0, 0, 0);
115/// assert_eq!(zzz, 0);
116/// ```
117///
118/// ```
119/// use utiles_core::xyz2rmid;
120/// let xyz_0_0_1 = xyz2rmid(0, 0, 1);
121/// assert_eq!(xyz_0_0_1, 1);
122/// ```
123///
124/// ```
125/// use utiles_core::xyz2rmid;
126/// let xyz_0_1_1 = xyz2rmid(1, 0, 1);
127/// assert_eq!(xyz_0_1_1, 2);
128/// ```
129///
130/// ```
131/// use utiles_core::xyz2rmid;
132/// let xyz_1_0_1 = xyz2rmid(0, 1, 1);
133/// assert_eq!(xyz_1_0_1, 3);
134/// ```
135///
136/// ```
137/// use utiles_core::xyz2rmid;
138/// let xyz_1_1_1 = xyz2rmid(1, 1, 1);
139/// assert_eq!(xyz_1_1_1, 4);
140/// ```
141///
142/// ```
143/// use utiles_core::xyz2rmid;
144/// let one_two_three = xyz2rmid(1, 2, 3);
145/// assert_eq!(one_two_three, 38);
146/// ```
147///
148/// ```
149/// use utiles_core::xyz2rmid;
150/// let last_tile_in_z12 = xyz2rmid(4095, 4095, 12);
151/// assert_eq!(last_tile_in_z12, 22369621 - 1); // total tiles thru z12 - 1
152/// ```
153#[must_use]
154pub fn xyz2rmid(x: u32, y: u32, z: u8) -> u64 {
155    if z == 0 {
156        return u64::from(x) + u64::from(y) * 2u64.pow(u32::from(z));
157    }
158    let base_id: u64 = (4u64.pow(u32::from(z)) - 1) / 3;
159    base_id + u64::from(x) + u64::from(y) * 2u64.pow(u32::from(z))
160}
161
162/// Calculate the xyz of the tile from a row-major-id
163///
164/// # Examples
165/// ```
166/// use utiles_core::rmid2xyz;
167/// let zzz = rmid2xyz(0);
168/// assert_eq!(zzz, (0, 0, 0));
169/// ```
170///
171/// ```
172/// use utiles_core::rmid2xyz;
173/// let xyz_0_0_1 = rmid2xyz(1);
174/// assert_eq!(xyz_0_0_1, (0, 0, 1));
175/// ```
176///
177/// ```
178/// use utiles_core::rmid2xyz;
179/// let xyz_0_1_1 = rmid2xyz(2);
180/// assert_eq!(xyz_0_1_1, (1, 0, 1));
181/// ```
182///
183/// ```
184/// use utiles_core::rmid2xyz;
185/// let xyz_1_0_1 = rmid2xyz(3);
186/// assert_eq!(xyz_1_0_1, (0, 1, 1));
187/// ```
188///
189/// ```
190/// use utiles_core::rmid2xyz;
191/// let xyz_1_1_1 = rmid2xyz(4);
192/// assert_eq!(xyz_1_1_1, (1, 1, 1));
193/// ```
194///
195/// ```
196/// use utiles_core::rmid2xyz;
197/// let one_two_three = rmid2xyz(38);
198/// assert_eq!(one_two_three, (1, 2, 3));
199/// ```
200///
201/// ```
202/// use utiles_core::rmid2xyz;
203/// let last_tile_in_z12 = rmid2xyz(22369621 - 1); // total tiles thru z12 - 1
204/// assert_eq!(last_tile_in_z12, (4095, 4095, 12));
205/// ```
206#[must_use]
207#[allow(clippy::cast_possible_truncation)]
208pub fn rmid2xyz(i: u64) -> (u32, u32, u8) {
209    if i == 0 {
210        return (0, 0, 0);
211    }
212    let (i_o, z) = int_2_offset_zoom(i);
213    let pow_z = 2u64.pow(u32::from(z));
214    let x = i_o % pow_z;
215    let y = i_o / pow_z;
216    (x as u32, y as u32, z)
217}
218
219/// Calculate the zoom level for the bounding-tile of a bbox
220#[must_use]
221pub fn bbox2zoom(bbox: (u32, u32, u32, u32)) -> u8 {
222    let max_zoom = 28;
223    let (west, south, east, north) = bbox;
224    for z in 0..max_zoom {
225        let mask = 1 << (32 - (z + 1));
226        if (west & mask) != (east & mask) || (south & mask) != (north & mask) {
227            return z;
228        }
229    }
230    max_zoom
231}
232
233/// Return the bbox tuple given x, y, z.
234#[must_use]
235pub fn bounds(x: u32, y: u32, z: u8) -> (f64, f64, f64, f64) {
236    let ul_corner = ul(x, y, z);
237    let lr_corner = ul(x + 1, y + 1, z);
238    (
239        ul_corner.lng(),
240        lr_corner.lat(),
241        lr_corner.lng(),
242        ul_corner.lat(),
243    )
244}
245
246/// Truncate a longitude to the valid range of -180 to 180.
247#[must_use]
248pub const fn truncate_lng(lng: f64) -> f64 {
249    lng.clamp(-180.0, 180.0)
250}
251
252/// Truncate a latitude to the valid range of -90 to 90.
253#[must_use]
254pub const fn truncate_lat(lat: f64) -> f64 {
255    lat.clamp(-90.0, 90.0)
256}
257
258/// Truncate a `LngLat` to valid range of longitude and latitude.
259#[must_use]
260pub fn truncate_lnglat(lnglat: &LngLat) -> LngLat {
261    LngLat::new(truncate_lng(lnglat.lng()), truncate_lat(lnglat.lat()))
262}
263
264/// Return the parent tile of a tile given x, y, z, and n (number of ancestors).
265#[must_use]
266pub fn parent(x: u32, y: u32, z: u8, n: Option<u8>) -> Option<Tile> {
267    let n = n.unwrap_or(0);
268    if n == 0 {
269        if z == 0 {
270            None
271        } else {
272            Some(utile!(x >> 1, y >> 1, z - 1))
273        }
274    } else {
275        parent(x >> 1, y >> 1, z - 1, Some(n - 1))
276    }
277}
278/// Return the 4 direct children of a tile
279#[must_use]
280pub fn children1_zorder(x: u32, y: u32, z: u8) -> [Tile; 4] {
281    [
282        utile!(x * 2, y * 2, z + 1),         // top-left
283        utile!(x * 2 + 1, y * 2, z + 1),     // top-right
284        utile!(x * 2, y * 2 + 1, z + 1),     // bottom-left
285        utile!(x * 2 + 1, y * 2 + 1, z + 1), // bottom-right
286    ]
287}
288/// Return the children of a tile given x, y, z, and zoom in z-order.
289///
290/// # Examples
291/// ```
292/// use utiles_core::{children_zorder, utile, Tile};
293/// let children = children_zorder(0, 0, 0, Some(1));
294/// assert_eq!(children.len(), 4);
295/// assert_eq!(children, vec![
296///     utile!(0, 0, 1),
297///     utile!(1, 0, 1),
298///     utile!(0, 1, 1),
299///     utile!(1, 1, 1),
300/// ]);
301/// ```
302///
303/// ```
304/// use utiles_core::{children_zorder, utile, Tile};
305/// let children = children_zorder(0, 0, 0, Some(2));
306/// assert_eq!(children.len(), 16);
307/// let expected = [
308///     utile!(0,0,2),
309///     utile!(1,0,2),
310///     utile!(0,1,2),
311///     utile!(1,1,2),
312///     utile!(2,0,2),
313///     utile!(3,0,2),
314///     utile!(2,1,2),
315///     utile!(3,1,2),
316///     utile!(0,2,2),
317///     utile!(1,2,2),
318///     utile!(0,3,2),
319///     utile!(1,3,2),
320///     utile!(2,2,2),
321///     utile!(3,2,2),
322///     utile!(2,3,2),
323///     utile!(3,3,2),
324/// ];
325/// assert_eq!(children, expected);
326/// ```
327///
328#[must_use]
329pub fn children_zorder(x: u32, y: u32, z: u8, zoom: Option<u8>) -> Vec<Tile> {
330    let zoom = zoom.unwrap_or(z + 1);
331    let tile = utile!(x, y, z);
332    let mut tiles = vec![tile];
333    while tiles[0].z < zoom {
334        let (xtile, ytile, ztile) = (tiles[0].x, tiles[0].y, tiles[0].z);
335        tiles.append(&mut vec![
336            utile!(xtile * 2, ytile * 2, ztile + 1), // top-left
337            utile!(xtile * 2 + 1, ytile * 2, ztile + 1), // top-right
338            utile!(xtile * 2, ytile * 2 + 1, ztile + 1), // bottom-left
339            utile!(xtile * 2 + 1, ytile * 2 + 1, ztile + 1), // bottom-right
340        ]);
341        tiles.remove(0);
342    }
343    tiles
344}
345
346/// Return the children of a tile given x, y, z, and zoom; returns children
347/// in stupid `a, b, d, c` orderl; but this is the mercantile way... and I
348/// am not gonna fix it right now
349///
350/// # Examples
351/// ```
352/// use utiles_core::{children, utile, Tile};
353/// let children = children(0, 0, 0, Some(1));
354/// assert_eq!(children.len(), 4);
355/// assert_eq!(children, vec![
356///     utile!(0, 0, 1),
357///     utile!(1, 0, 1),
358///     utile!(1, 1, 1),
359///     utile!(0, 1, 1),
360/// ]);
361/// ```
362#[must_use]
363pub fn children(x: u32, y: u32, z: u8, zoom: Option<u8>) -> Vec<Tile> {
364    let zoom = zoom.unwrap_or(z + 1);
365    let tile = utile!(x, y, z);
366    let mut tiles = vec![tile];
367    while tiles[0].z < zoom {
368        let (xtile, ytile, ztile) = (tiles[0].x, tiles[0].y, tiles[0].z);
369        tiles.append(&mut vec![
370            utile!(xtile * 2, ytile * 2, ztile + 1), // top-left
371            utile!(xtile * 2 + 1, ytile * 2, ztile + 1), // top-right
372            utile!(xtile * 2 + 1, ytile * 2 + 1, ztile + 1), // bottom-right
373            utile!(xtile * 2, ytile * 2 + 1, ztile + 1), // bottom-left
374        ]);
375        tiles.remove(0);
376    }
377    tiles
378}
379
380/// Return the siblings of a tile given x, y, z
381///
382/// Siblings are tiles that share the same parent, NOT neighbors.
383#[must_use]
384pub fn siblings(x: u32, y: u32, z: u8) -> Vec<Tile> {
385    let sibrel = SiblingRelationship::from((x, y));
386    match sibrel {
387        SiblingRelationship::UpperLeft => vec![
388            utile!(x + 1, y, z),
389            utile!(x, y + 1, z),
390            utile!(x + 1, y + 1, z),
391        ],
392        SiblingRelationship::UpperRight => vec![
393            utile!(x - 1, y, z),
394            utile!(x, y + 1, z),
395            utile!(x - 1, y + 1, z),
396        ],
397        SiblingRelationship::LowerLeft => vec![
398            utile!(x + 1, y, z),
399            utile!(x, y - 1, z),
400            utile!(x + 1, y - 1, z),
401        ],
402        SiblingRelationship::LowerRight => vec![
403            utile!(x - 1, y, z),
404            utile!(x, y - 1, z),
405            utile!(x - 1, y - 1, z),
406        ],
407    }
408}
409
410/// Truncate a bounding box to the valid range of longitude and latitude.
411#[must_use]
412pub fn bbox_truncate(
413    west: f64,
414    south: f64,
415    east: f64,
416    north: f64,
417    truncate: Option<bool>,
418) -> (f64, f64, f64, f64) {
419    let trunc = truncate.unwrap_or(false);
420    let mut west = west;
421    let mut east = east;
422    let mut south = south;
423    let mut north = north;
424    if trunc {
425        if west < -180.0 {
426            west = -180.0;
427        }
428        if east > 180.0 {
429            east = 180.0;
430        }
431        if south < -90.0 {
432            south = -90.0;
433        }
434        if north > 90.0 {
435            north = 90.0;
436        }
437    }
438    (west, south, east, north)
439}
440
441/// Convert lng lat to web mercator x and y
442///
443/// # Errors
444///
445/// Returns error if y can not be computed.
446pub fn _xy(lng: f64, lat: f64, truncate: Option<bool>) -> UtilesCoreResult<(f64, f64)> {
447    let (lng, lat) = if truncate.unwrap_or(false) {
448        (truncate_lng(lng), truncate_lat(lat))
449    } else {
450        (lng, lat)
451    };
452    let sinlat = lat.to_radians().sin();
453    let yish = (1.0 + sinlat) / (1.0 - sinlat);
454    match yish.classify() {
455        FpCategory::Infinite | FpCategory::Nan => {
456            Err(UtilesCoreError::LngLat2WebMercator(format!(
457                "Y can not be computed: lat={lat}"
458            )))
459        }
460        _ => {
461            let y = 0.5 - 0.25 * (yish.ln()) / PI;
462            let x = lng / 360.0 + 0.5;
463            Ok((x, y))
464        }
465    }
466}
467
468/// Convert lng lat to web mercator x and y
469#[must_use]
470pub fn lnglat2webmercator(lng: f64, lat: f64) -> (f64, f64) {
471    let x = EARTH_RADIUS * lng.to_radians();
472    let y = if (lat - 90.0).abs() < f64::EPSILON {
473        f64::INFINITY
474    } else if (lat + 90.0).abs() < f64::EPSILON {
475        f64::NEG_INFINITY
476    } else {
477        EARTH_RADIUS * PI.mul_add(0.25, 0.5 * lat.to_radians()).tan().ln()
478    };
479    (x, y)
480}
481
482/// Convert web mercator x and y to longitude and latitude.
483///
484/// # Examples
485/// ```
486/// use utiles_core::webmercator2lnglat;
487/// let (lng, lat) = webmercator2lnglat(0.5, 0.5);
488/// assert!((lng - 0.0).abs() < 0.0001, "lng: {}", lng);
489/// assert!((lat - 0.0).abs() < 0.0001, "lat: {}", lat);
490/// ```
491///
492#[must_use]
493#[inline]
494pub fn webmercator2lnglat(x: f64, y: f64) -> (f64, f64) {
495    let lng = (x / EARTH_RADIUS).to_degrees();
496    let lat = 2.0f64
497        .mul_add((y / EARTH_RADIUS).exp().atan(), -(PI * 0.5))
498        .to_degrees();
499    (lng, lat)
500}
501
502/// Convert longitude and latitude to web mercator x and y with optional truncation.
503///
504/// Name "xy" comes from mercantile python library.
505#[must_use]
506pub fn xy(lng: f64, lat: f64, truncate: Option<bool>) -> (f64, f64) {
507    let (lng, lat) = if truncate.unwrap_or(false) {
508        (truncate_lng(lng), truncate_lat(lat))
509    } else {
510        (lng, lat)
511    };
512    lnglat2webmercator(lng, lat)
513}
514
515/// Convert web mercator x and y to longitude and latitude with optional truncation.
516#[must_use]
517pub fn lnglat(x: f64, y: f64, truncate: Option<bool>) -> LngLat {
518    let (lng, lat) = webmercator2lnglat(x, y);
519    if truncate.is_some() {
520        truncate_lnglat(&LngLat::new(lng, lat))
521    } else {
522        LngLat::new(lng, lat)
523    }
524}
525
526enum TileEdgeInfo {
527    Bottom,
528    BottomLeft,
529    BottomRight,
530    Left,
531    Right,
532    Top,
533    TopLeft,
534    TopRight,
535    Middle,
536}
537
538fn tile_edge_info(x: u32, y: u32, z: u8) -> TileEdgeInfo {
539    if x == 0 && y == 0 {
540        return TileEdgeInfo::TopLeft;
541    }
542    let max_xy = 2u32.pow(u32::from(z));
543    if x == max_xy && y == max_xy {
544        return TileEdgeInfo::BottomRight;
545    }
546    match (x, y) {
547        (max, 0) if max == max_xy => TileEdgeInfo::TopRight,
548        (0, max) if max == max_xy => TileEdgeInfo::BottomLeft,
549        (0, _) => TileEdgeInfo::Left,
550        (max, _) if max == max_xy => TileEdgeInfo::Right,
551        (_, 0) => TileEdgeInfo::Top,
552        (_, max) if max == max_xy => TileEdgeInfo::Bottom,
553        _ => TileEdgeInfo::Middle,
554    }
555}
556
557fn neighbors_middle_tile(x: u32, y: u32, z: u8) -> Vec<Tile> {
558    vec![
559        utile!(x + 1, y, z),
560        utile!(x, y + 1, z),
561        utile!(x + 1, y + 1, z),
562        utile!(x - 1, y, z),
563        utile!(x, y - 1, z),
564        utile!(x - 1, y - 1, z),
565        utile!(x + 1, y - 1, z),
566        utile!(x - 1, y + 1, z),
567    ]
568}
569
570/// Return neighbors of a tile (non-wrapping).
571#[must_use]
572pub fn neighbors(x: u32, y: u32, z: u8) -> Vec<Tile> {
573    if z == 0 {
574        Vec::new()
575    } else if z == 1 {
576        siblings(x, y, z)
577    } else {
578        let edge_info = tile_edge_info(x, y, z);
579        match edge_info {
580            TileEdgeInfo::Middle => neighbors_middle_tile(x, y, z),
581            TileEdgeInfo::TopLeft => vec![
582                utile!(x + 1, y, z),
583                utile!(x, y + 1, z),
584                utile!(x + 1, y + 1, z),
585            ],
586            TileEdgeInfo::TopRight => vec![
587                utile!(x - 1, y, z),
588                utile!(x, y + 1, z),
589                utile!(x - 1, y + 1, z),
590            ],
591            TileEdgeInfo::BottomLeft => vec![
592                utile!(x + 1, y, z),
593                utile!(x, y - 1, z),
594                utile!(x + 1, y - 1, z),
595            ],
596            TileEdgeInfo::BottomRight => vec![
597                utile!(x - 1, y, z),
598                utile!(x, y - 1, z),
599                utile!(x - 1, y - 1, z),
600            ],
601            TileEdgeInfo::Left => vec![
602                utile!(x + 1, y, z),
603                utile!(x, y + 1, z),
604                utile!(x + 1, y + 1, z),
605                utile!(x, y - 1, z),
606                utile!(x + 1, y - 1, z),
607            ],
608            TileEdgeInfo::Right => vec![
609                utile!(x - 1, y, z),
610                utile!(x, y + 1, z),
611                utile!(x - 1, y + 1, z),
612                utile!(x, y - 1, z),
613                utile!(x - 1, y - 1, z),
614            ],
615            TileEdgeInfo::Top => vec![
616                utile!(x + 1, y, z),
617                utile!(x, y + 1, z),
618                utile!(x + 1, y + 1, z),
619                utile!(x - 1, y, z),
620                utile!(x - 1, y + 1, z),
621            ],
622            TileEdgeInfo::Bottom => vec![
623                utile!(x + 1, y, z),
624                utile!(x, y - 1, z),
625                utile!(x + 1, y - 1, z),
626                utile!(x - 1, y, z),
627                utile!(x - 1, y - 1, z),
628            ],
629        }
630    }
631}
632
633static NEIGHBOR_IDXS: &[(i64, i64)] = &[
634    (-1, -1),
635    (-1, 0),
636    (-1, 1),
637    (0, -1),
638    (0, 1),
639    (1, -1),
640    (1, 0),
641    (1, 1),
642];
643
644#[must_use]
645#[expect(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
646pub fn neighbors_wrap_x(x: u32, y: u32, z: u8) -> Vec<Tile> {
647    if z == 0 {
648        return Vec::new();
649    }
650    if z == 1 {
651        return siblings(x, y, z);
652    }
653
654    let max_xy = 1u32 << z; // 2^z
655    let max_xy_i = i64::from(max_xy);
656
657    NEIGHBOR_IDXS
658        .iter()
659        .filter_map(|&(dx, dy)| {
660            let nx_i = (i64::from(x) + dx).rem_euclid(max_xy_i);
661            let ny_i = i64::from(y) + dy;
662
663            if ny_i < 0 || ny_i >= max_xy_i {
664                return None; // off the top/bottom edge
665            }
666
667            Some(utile!(nx_i as u32, ny_i as u32, z))
668        })
669        .collect()
670}
671
672/// Return Tile struct from longitude, latitude, and zoom.
673///
674/// # Errors
675///
676/// Returns error if the lnglat can not be converted to web mercator.
677pub fn tile(
678    lng: f64,
679    lat: f64,
680    zoom: u8,
681    truncate: Option<bool>,
682) -> Result<Tile, UtilesCoreError> {
683    Tile::from_lnglat_zoom(lng, lat, zoom, truncate)
684}
685
686/// Converts longitude, latitude, and zoom level to fractional tile coordinates.
687///
688/// # Examples
689/// ```
690/// use utiles_core::lnglat2tile_frac;
691/// let (xf, yf, z) =lnglat2tile_frac(-95.939_655_303_955_08, 41.260_001_085_686_97, 9);
692/// assert!((xf - 119.552_490_234_375).abs() < 0.0001, "xf: {}", xf);
693/// assert!((yf - 191.471_191_406_25).abs() < 0.0001, "yf: {}", yf);
694/// assert!(z == 9);
695/// ```
696#[must_use]
697pub fn lnglat2tile_frac(lng: f64, lat: f64, z: u8) -> (f64, f64, u8) {
698    let sin = (lat * DEG2RAD).sin();
699    let z2 = 2f64.powi(i32::from(z));
700    let mut x = z2 * (lng / 360.0 + 0.5);
701    let y = z2 * (0.5 - (0.25 * ((1.0 + sin) / (1.0 - sin)).ln()) / PI);
702
703    // Wrap Tile X using rem_euclid
704    x = x.rem_euclid(z2);
705
706    (x, y, z)
707}
708
709/// Return the bounding tile for a bounding box.
710///
711/// # Errors
712///
713/// Returns error if the bounding tile can not be calculated for points on the bbox
714pub fn bounding_tile(
715    bbox: BBox,
716    truncate: Option<bool>,
717) -> Result<Tile, UtilesCoreError> {
718    let (west, south, east, north) =
719        bbox_truncate(bbox.west, bbox.south, bbox.east, bbox.north, truncate);
720    let tmin = tile(west, north, 32, truncate)?;
721    let tmax = tile(east - LL_EPSILON, south + LL_EPSILON, 32, truncate)?;
722
723    let cell = (tmin.x, tmin.y, tmax.x, tmax.y);
724    let z = bbox2zoom(cell);
725    if z == 0 {
726        Ok(utile!(0, 0, 0))
727    } else {
728        let x = cell.0 >> (32 - z);
729        let y = cell.1 >> (32 - z);
730        Ok(utile!(x, y, z))
731    }
732}
733
734/// Return web-mercator bbox from x, y, z.
735#[must_use]
736pub fn xyz2bbox(x: u32, y: u32, z: u8) -> WebBBox {
737    let tile_size = EARTH_CIRCUMFERENCE / 2.0_f64.powi(i32::from(z));
738    let left = f64::from(x).mul_add(tile_size, -(EARTH_CIRCUMFERENCE / 2.0));
739    let top = f64::from(y).mul_add(-tile_size, EARTH_CIRCUMFERENCE / 2.0);
740
741    WebBBox::new(left, top - tile_size, left + tile_size, top)
742}
743
744/// Return zooms-vec from a `ZoomOrZooms` enum
745#[must_use]
746pub fn as_zooms(zoom_or_zooms: ZoomOrZooms) -> Vec<u8> {
747    match zoom_or_zooms {
748        ZoomOrZooms::Zoom(zoom) => {
749            vec![zoom]
750        }
751        ZoomOrZooms::Zooms(zooms) => zooms,
752    }
753}
754
755fn tiles_range_zoom(
756    minx: u32,
757    maxx: u32,
758    miny: u32,
759    maxy: u32,
760    zoom: u8,
761) -> impl Iterator<Item = (u32, u32, u8)> {
762    (minx..=maxx).flat_map(move |i| (miny..=maxy).map(move |j| (i, j, zoom)))
763}
764
765/// Return `TileRanges` from a bounding box and zoom(s).
766///
767/// # Errors
768///
769/// Returns error tiles cannot be calculated due to invalid bbox/zbox
770pub fn tile_ranges(
771    bounds: (f64, f64, f64, f64),
772    zooms: ZoomOrZooms,
773) -> Result<TileZBoxes, UtilesCoreError> {
774    let zooms = as_zooms(zooms);
775    let bboxes = BBox::from(bounds).bboxes().into_iter().map(|bbox| {
776        // clip to web mercator extent
777        BBox {
778            north: bbox.north.min(85.051_129),
779            south: bbox.south.max(-85.051_129),
780            east: bbox.east.min(180.0),
781            west: bbox.west.max(-180.0),
782        }
783    });
784    let ranges: Vec<TileZBox> = bboxes
785        .into_iter()
786        .flat_map(move |bbox| {
787            let zooms = zooms.clone();
788            zooms.into_iter().map(move |zoom| {
789                let upper_left_lnglat = LngLat {
790                    xy: point2d! { x: bbox.west, y: bbox.north },
791                };
792                let lower_right_lnglat = LngLat {
793                    xy: point2d! { x: bbox.east, y: bbox.south },
794                };
795                let top_left_tile = Tile::from_lnglat_zoom(
796                    upper_left_lnglat.lng(),
797                    upper_left_lnglat.lat(),
798                    zoom,
799                    Some(false),
800                )?;
801                let bottom_right_tile = Tile::from_lnglat_zoom(
802                    lower_right_lnglat.lng() - LL_EPSILON,
803                    lower_right_lnglat.lat() + LL_EPSILON,
804                    zoom,
805                    Some(false),
806                )?;
807                Ok(TileZBox::new(
808                    top_left_tile.x,
809                    bottom_right_tile.x,
810                    top_left_tile.y,
811                    bottom_right_tile.y,
812                    zoom,
813                ))
814            })
815        })
816        .collect::<Result<Vec<TileZBox>, UtilesCoreError>>()?;
817    Ok(TileZBoxes::from(ranges))
818}
819
820/// Return the number of tiles for a bounding box and zoom(s).
821///
822/// # Errors
823///
824/// Returns error if the number of tiles cannot be calculated due to invalid bbox/zbox
825pub fn tiles_count(
826    bounds: (f64, f64, f64, f64),
827    zooms: ZoomOrZooms,
828) -> Result<u64, UtilesCoreError> {
829    let ranges = tile_ranges(bounds, zooms)?;
830    Ok(ranges.length())
831}
832
833/// Return an iterator of tiles for a bounding box and zoom(s).
834pub fn tiles(
835    bounds: (f64, f64, f64, f64),
836    zooms: ZoomOrZooms,
837) -> impl Iterator<Item = Tile> {
838    let zooms = as_zooms(zooms);
839    let bboxthing = BBox {
840        north: bounds.3,
841        south: bounds.1,
842        east: bounds.2,
843        west: bounds.0,
844    };
845    let bboxes = bboxthing.bboxes().into_iter().map(|bbox| {
846        // clip to web mercator extent
847        BBox {
848            north: bbox.north.min(85.051_129),
849            south: bbox.south.max(-85.051_129),
850            east: bbox.east.min(180.0),
851            west: bbox.west.max(-180.0),
852        }
853    });
854    bboxes.into_iter().flat_map(move |bbox| {
855        let zooms = zooms.clone();
856        zooms.into_iter().flat_map(move |zoom| {
857            let upper_left_lnglat = LngLat {
858                xy: point2d! { x: bbox.west, y: bbox.north },
859            };
860            let lower_right_lnglat = LngLat {
861                xy: point2d! { x: bbox.east, y: bbox.south },
862            };
863            let top_left_tile = Tile::from_lnglat_zoom(
864                upper_left_lnglat.lng(),
865                upper_left_lnglat.lat(),
866                zoom,
867                Some(false),
868            );
869            let bottom_right_tile = Tile::from_lnglat_zoom(
870                lower_right_lnglat.lng() - LL_EPSILON,
871                lower_right_lnglat.lat() + LL_EPSILON,
872                zoom,
873                Some(false),
874            );
875
876            match (top_left_tile, bottom_right_tile) {
877                (Ok(top_left), Ok(bottom_right)) => tiles_range_zoom(
878                    top_left.x,
879                    bottom_right.x,
880                    top_left.y,
881                    bottom_right.y,
882                    zoom,
883                )
884                .map(move |(x, y, z)| Tile { x, y, z })
885                .collect::<Vec<_>>()
886                .into_iter(),
887                _ => Vec::new().into_iter(),
888            }
889        })
890    })
891}
892
893/// Convert tile xyz to u64 tile id (based on mapbox coverage implementation)
894#[must_use]
895pub fn to_id(x: u32, y: u32, z: u8) -> u64 {
896    let dim = 2u64 * (1u64 << z);
897    ((dim * u64::from(y) + u64::from(x)) * 32u64) + u64::from(z)
898}
899
900/// Convert tile u64 id to tile xyz
901///
902/// # Panics
903///
904/// Errors on integer conversion error (should not happen) should not happen
905#[allow(clippy::cast_possible_truncation)]
906#[must_use]
907pub fn from_id(id: u64) -> Tile {
908    let z = (id % 32) as u8;
909    let dim = 2u64 * (1u64 << z);
910    let xy = (id - u64::from(z)) / 32u64;
911    let x = u32::try_from(xy % dim).expect("should never happen");
912    let y = ((xy - u64::from(x)) / dim) as u32;
913    utile!(x, y, z)
914}