Fix bbox_to_xyz (#1070)

Figure out and fix the bug with `tile-grid` crate or remove the crate and implement it ourself. --------- Co-authored-by: Yuri Astrakhan <[email protected]>
maplibre · Dec 18, 2023 · f21119d · f21119d
1 parent c650d4b
commit f21119d
Show file tree

Hide file tree

Showing 6 changed files with 133 additions and 103 deletions.
diff --git a/Cargo.lock b/Cargo.lock
diff --git a/martin-tile-utils/Cargo.toml b/martin-tile-utils/Cargo.toml
@@ -17,7 +17,7 @@ repository.workspace = true
 rust-version.workspace = true
 
 [dependencies]
-tile-grid.workspace = true
 
 [dev-dependencies]
 approx.workspace = true
+insta.workspace = true
diff --git a/martin-tile-utils/src/lib.rs b/martin-tile-utils/src/lib.rs
@@ -6,18 +6,10 @@
 use std::f64::consts::PI;
 use std::fmt::Display;
 
-use tile_grid::{tms, Tms, Xyz};
-
 pub const EARTH_CIRCUMFERENCE: f64 = 40_075_016.685_578_5;
 pub const EARTH_RADIUS: f64 = EARTH_CIRCUMFERENCE / 2.0 / PI;
 
-pub const MAX_ZOOM: u8 = 24;
-use std::sync::OnceLock;
-
-fn web_merc() -> &'static Tms {
-    static TMS: OnceLock<Tms> = OnceLock::new();
-    TMS.get_or_init(|| tms().lookup("WebMercatorQuad").unwrap())
-}
+pub const MAX_ZOOM: u8 = 30;
 
 #[derive(Clone, Copy, Debug, PartialEq, Eq)]
 pub enum Format {
@@ -196,37 +188,46 @@ impl Display for TileInfo {
 
 /// Convert longitude and latitude to a tile (x,y) coordinates for a given zoom
 #[must_use]
-#[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
-pub fn tile_index(lon: f64, lat: f64, zoom: u8) -> (u32, u32) {
-    assert!(zoom <= MAX_ZOOM, "zoom {zoom} must be <= {MAX_ZOOM}");
-    let tile = web_merc().tile(lon, lat, zoom).unwrap();
-    let max_value = (1_u64 << zoom) - 1;
-    (tile.x.min(max_value) as u32, tile.y.min(max_value) as u32)
+#[allow(clippy::cast_possible_truncation)]
+#[allow(clippy::cast_sign_loss)]
+pub fn tile_index(lng: f64, lat: f64, zoom: u8) -> (u32, u32) {
+    let tile_size = EARTH_CIRCUMFERENCE / f64::from(1_u32 << zoom);
+    let (x, y) = wgs84_to_webmercator(lng, lat);
+    let col = (((x - (EARTH_CIRCUMFERENCE * -0.5)).abs() / tile_size) as u32).min((1 << zoom) - 1);
+    let row = ((((EARTH_CIRCUMFERENCE * 0.5) - y).abs() / tile_size) as u32).min((1 << zoom) - 1);
+    (col, row)
 }
 
 /// Convert min/max XYZ tile coordinates to a bounding box values.
 /// The result is `[min_lng, min_lat, max_lng, max_lat]`
 #[must_use]
 pub fn xyz_to_bbox(zoom: u8, min_x: u32, min_y: u32, max_x: u32, max_y: u32) -> [f64; 4] {
     assert!(zoom <= MAX_ZOOM, "zoom {zoom} must be <= {MAX_ZOOM}");
-    assert!(min_x <= max_x, "min_x {min_x} must be <= max_x {max_x}");
-    assert!(min_y <= max_y, "min_y {min_y} must be <= max_y {max_y}");
-    let left_top_bounds = web_merc().xy_bounds(&Xyz::new(u64::from(min_x), u64::from(min_y), zoom));
-    let right_bottom_bounds =
-        web_merc().xy_bounds(&Xyz::new(u64::from(max_x), u64::from(max_y), zoom));
-    let (min_lng, min_lat) = webmercator_to_wgs84(left_top_bounds.left, right_bottom_bounds.bottom);
-    let (max_lng, max_lat) = webmercator_to_wgs84(right_bottom_bounds.right, left_top_bounds.top);
 
+    let tile_length = EARTH_CIRCUMFERENCE / f64::from(1_u32 << zoom);
+
+    let left_down_bbox = tile_bbox(min_x, max_y, tile_length);
+    let right_top_bbox = tile_bbox(max_x, min_y, tile_length);
+
+    let (min_lng, min_lat) = webmercator_to_wgs84(left_down_bbox[0], left_down_bbox[1]);
+    let (max_lng, max_lat) = webmercator_to_wgs84(right_top_bbox[2], right_top_bbox[3]);
     [min_lng, min_lat, max_lng, max_lat]
 }
 
+#[allow(clippy::cast_lossless)]
+fn tile_bbox(x: u32, y: u32, tile_length: f64) -> [f64; 4] {
+    let min_x = EARTH_CIRCUMFERENCE * -0.5 + x as f64 * tile_length;
+    let max_y = EARTH_CIRCUMFERENCE * 0.5 - y as f64 * tile_length;
+
+    [min_x, max_y - tile_length, min_x + tile_length, max_y]
+}
+
 /// Convert bounding box to a tile box `(min_x, min_y, max_x, max_y)` for a given zoom
 #[must_use]
-#[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
 pub fn bbox_to_xyz(left: f64, bottom: f64, right: f64, top: f64, zoom: u8) -> (u32, u32, u32, u32) {
-    let (min_x, min_y) = tile_index(left, top, zoom);
-    let (max_x, max_y) = tile_index(right, bottom, zoom);
-    (min_x, min_y, max_x, max_y)
+    let (min_col, min_row) = tile_index(left, top, zoom);
+    let (max_col, max_row) = tile_index(right, bottom, zoom);
+    (min_col, min_row, max_col, max_row)
 }
 
 /// Compute precision of a zoom level, i.e. how many decimal digits of the longitude and latitude are relevant
@@ -250,13 +251,27 @@ pub fn webmercator_to_wgs84(x: f64, y: f64) -> (f64, f64) {
     (lng, lat)
 }
 
+/// transform WGS84 to `WebMercator`
+// from https://github.com/Esri/arcgis-osm-editor/blob/e4b9905c264aa22f8eeb657efd52b12cdebea69a/src/OSMWeb10_1/Utils/WebMercator.cs
+#[must_use]
+pub fn wgs84_to_webmercator(lon: f64, lat: f64) -> (f64, f64) {
+    let x = lon * PI / 180.0 * EARTH_RADIUS;
+
+    let rad = lat * PI / 180.0;
+    let sin = rad.sin();
+    let y = EARTH_RADIUS / 2.0 * ((1.0 + sin) / (1.0 - sin)).ln();
+
+    (x, y)
+}
+
 #[cfg(test)]
 mod tests {
     #![allow(clippy::unreadable_literal)]
 
     use std::fs::read;
 
     use approx::assert_relative_eq;
+    use insta::assert_snapshot;
     use Encoding::{Internal, Uncompressed};
     use Format::{Jpeg, Json, Png, Webp};
 
@@ -296,61 +311,75 @@ mod tests {
     }
 
     #[test]
-    fn test_tile_index() {
+    fn test_tile_colrow() {
         assert_eq!((0, 0), tile_index(-180.0, 85.0511, 0));
     }
 
     #[test]
     fn test_xyz_to_bbox() {
         // you could easily get test cases from maptiler: https://www.maptiler.com/google-maps-coordinates-tile-bounds-projection/#4/-118.82/71.02
         let bbox = xyz_to_bbox(0, 0, 0, 0, 0);
-        assert_relative_eq!(bbox[0], -179.99999999999955, epsilon = f64::EPSILON * 2.0);
+        assert_relative_eq!(bbox[0], -180.0, epsilon = f64::EPSILON * 2.0);
         assert_relative_eq!(bbox[1], -85.0511287798066, epsilon = f64::EPSILON * 2.0);
-        assert_relative_eq!(bbox[2], 179.99999999999986, epsilon = f64::EPSILON * 2.0);
-        assert_relative_eq!(bbox[3], 85.05112877980655, epsilon = f64::EPSILON * 2.0);
-
-        let xyz = bbox_to_xyz(bbox[0], bbox[1], bbox[2], bbox[3], 0);
-        assert_eq!(xyz, (0, 0, 0, 0));
+        assert_relative_eq!(bbox[2], 180.0, epsilon = f64::EPSILON * 2.0);
+        assert_relative_eq!(bbox[3], 85.0511287798066, epsilon = f64::EPSILON * 2.0);
 
         let bbox = xyz_to_bbox(1, 0, 0, 0, 0);
-        assert_relative_eq!(bbox[0], -179.99999999999955, epsilon = f64::EPSILON * 2.0);
-        assert_relative_eq!(bbox[1], 2.007891127734306e-13, epsilon = f64::EPSILON * 2.0);
-        assert_relative_eq!(
-            bbox[2],
-            -2.007891127734306e-13,
-            epsilon = f64::EPSILON * 2.0
-        );
-        assert_relative_eq!(bbox[3], 85.05112877980655, epsilon = f64::EPSILON * 2.0);
-
-        let xyz = bbox_to_xyz(bbox[0], bbox[1], bbox[2], bbox[3], 1);
-        assert!(xyz.0 == 0 || xyz.0 == 1);
-        assert!(xyz.1 == 0);
-        assert!(xyz.2 == 0 || xyz.2 == 1);
-        assert!(xyz.3 == 0 || xyz.3 == 1);
+        assert_relative_eq!(bbox[0], -180.0, epsilon = f64::EPSILON * 2.0);
+        assert_relative_eq!(bbox[1], 0.0, epsilon = f64::EPSILON * 2.0);
+        assert_relative_eq!(bbox[2], 0.0, epsilon = f64::EPSILON * 2.0);
+        assert_relative_eq!(bbox[3], 85.0511287798066, epsilon = f64::EPSILON * 2.0);
 
         let bbox = xyz_to_bbox(5, 1, 1, 2, 2);
-        assert_relative_eq!(bbox[0], -168.74999999999955, epsilon = f64::EPSILON * 2.0);
-        assert_relative_eq!(bbox[1], 81.09321385260832, epsilon = f64::EPSILON * 2.0);
-        assert_relative_eq!(bbox[2], -146.2499999999996, epsilon = f64::EPSILON * 2.0);
-        assert_relative_eq!(bbox[3], 83.979259498862, epsilon = f64::EPSILON * 2.0);
-
-        let xyz = bbox_to_xyz(bbox[0], bbox[1], bbox[2], bbox[3], 5);
-        assert!(xyz.0 == 1 || xyz.0 == 2);
-        assert!(xyz.1 == 0 || xyz.1 == 1);
-        assert!(xyz.2 == 2 || xyz.2 == 3);
-        assert!(xyz.3 == 2 || xyz.3 == 3);
+        assert_relative_eq!(bbox[0], -168.75, epsilon = f64::EPSILON * 2.0);
+        assert_relative_eq!(bbox[1], 81.09321385260837, epsilon = f64::EPSILON * 2.0);
+        assert_relative_eq!(bbox[2], -146.25, epsilon = f64::EPSILON * 2.0);
+        assert_relative_eq!(bbox[3], 83.97925949886205, epsilon = f64::EPSILON * 2.0);
 
         let bbox = xyz_to_bbox(5, 1, 3, 2, 5);
-        assert_relative_eq!(bbox[0], -168.74999999999955, epsilon = f64::EPSILON * 2.0);
-        assert_relative_eq!(bbox[1], 74.01954331150218, epsilon = f64::EPSILON * 2.0);
-        assert_relative_eq!(bbox[2], -146.2499999999996, epsilon = f64::EPSILON * 2.0);
-        assert_relative_eq!(bbox[3], 81.09321385260832, epsilon = f64::EPSILON * 2.0);
-
-        let xyz = bbox_to_xyz(bbox[0], bbox[1], bbox[2], bbox[3], 5);
-        assert!(xyz.0 == 1 || xyz.0 == 2);
-        assert!(xyz.1 == 2 || xyz.1 == 3);
-        assert!(xyz.2 == 2 || xyz.2 == 3);
-        assert!(xyz.3 == 5 || xyz.3 == 6);
+        assert_relative_eq!(bbox[0], -168.75, epsilon = f64::EPSILON * 2.0);
+        assert_relative_eq!(bbox[1], 74.01954331150226, epsilon = f64::EPSILON * 2.0);
+        assert_relative_eq!(bbox[2], -146.25, epsilon = f64::EPSILON * 2.0);
+        assert_relative_eq!(bbox[3], 81.09321385260837, epsilon = f64::EPSILON * 2.0);
+    }
+
+    #[test]
+    fn test_box_to_xyz() {
+        fn tst(left: f64, bottom: f64, right: f64, top: f64, zoom: u8) -> String {
+            let (x0, y0, x1, y1) = bbox_to_xyz(left, bottom, right, top, zoom);
+            format!("({x0}, {y0}, {x1}, {y1})")
+        }
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 0), @"(0, 0, 0, 0)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 1), @"(0, 1, 0, 1)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 2), @"(0, 3, 0, 3)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 3), @"(0, 7, 0, 7)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 4), @"(0, 14, 1, 15)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 5), @"(0, 29, 2, 31)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 6), @"(0, 58, 5, 63)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 7), @"(0, 116, 11, 126)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 8), @"(0, 233, 23, 253)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 9), @"(0, 466, 47, 507)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 10), @"(1, 933, 94, 1014)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 11), @"(3, 1866, 188, 2029)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 12), @"(6, 3732, 377, 4059)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 13), @"(12, 7465, 755, 8119)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 14), @"(25, 14931, 1510, 16239)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 15), @"(51, 29863, 3020, 32479)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 16), @"(102, 59727, 6041, 64958)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 17), @"(204, 119455, 12083, 129917)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 18), @"(409, 238911, 24166, 259834)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 19), @"(819, 477823, 48332, 519669)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 20), @"(1638, 955647, 96665, 1039339)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 21), @"(3276, 1911295, 193331, 2078678)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 22), @"(6553, 3822590, 386662, 4157356)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 23), @"(13107, 7645181, 773324, 8314713)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 24), @"(26214, 15290363, 1546649, 16629427)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 25), @"(52428, 30580726, 3093299, 33258855)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 26), @"(104857, 61161453, 6186598, 66517711)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 27), @"(209715, 122322907, 12373196, 133035423)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 28), @"(419430, 244645814, 24746393, 266070846)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 29), @"(838860, 489291628, 49492787, 532141692)");
+        assert_snapshot!(tst(-179.43749999999955,-84.76987877980656,-146.8124999999996,-81.37446385260833, 30), @"(1677721, 978583256, 98985574, 1064283385)");
     }
 
     #[test]

diff --git a/mbtiles/src/summary.rs b/mbtiles/src/summary.rs
@@ -239,10 +239,10 @@ mod tests {
         max_tile_size: 1107
         avg_tile_size: 96.2295918367347
         bbox:
-          - -179.99999999999955
-          - -85.05112877980659
-          - 180.00000000000028
-          - 85.05112877980655
+          - -180
+          - -85.0511287798066
+          - 180.00000000000003
+          - 85.0511287798066
         min_zoom: 0
         max_zoom: 6
         zoom_info:
@@ -252,70 +252,70 @@ mod tests {
             max_tile_size: 1107
             avg_tile_size: 1107
             bbox:
-              - -179.99999999999955
-              - -85.05112877980659
-              - 179.99999999999986
-              - 85.05112877980655
+              - -180
+              - -85.0511287798066
+              - 180
+              - 85.0511287798066
           - zoom: 1
             tile_count: 4
             min_tile_size: 160
             max_tile_size: 650
             avg_tile_size: 366.5
             bbox:
-              - -179.99999999999955
-              - -85.05112877980652
-              - 179.99999999999915
-              - 85.05112877980655
+              - -180
+              - -85.0511287798066
+              - 180
+              - 85.0511287798066
           - zoom: 2
             tile_count: 7
             min_tile_size: 137
             max_tile_size: 495
             avg_tile_size: 239.57142857142858
             bbox:
-              - -179.99999999999955
-              - -66.51326044311165
-              - 179.99999999999915
-              - 66.51326044311182
+              - -180
+              - -66.51326044311186
+              - 180.00000000000003
+              - 66.51326044311186
           - zoom: 3
             tile_count: 17
             min_tile_size: 67
             max_tile_size: 246
             avg_tile_size: 134
             bbox:
-              - -134.99999999999957
-              - -40.979898069620376
-              - 180.00000000000028
-              - 66.51326044311169
+              - -135
+              - -40.97989806962013
+              - 179.99999999999997
+              - 66.51326044311186
           - zoom: 4
             tile_count: 38
             min_tile_size: 64
             max_tile_size: 175
             avg_tile_size: 86
             bbox:
-              - -134.99999999999963
-              - -40.979898069620106
-              - 179.99999999999966
-              - 66.51326044311175
+              - -135
+              - -40.97989806962014
+              - 180.00000000000003
+              - 66.51326044311186
           - zoom: 5
             tile_count: 57
             min_tile_size: 64
             max_tile_size: 107
             avg_tile_size: 72.7719298245614
             bbox:
-              - -123.74999999999966
-              - -40.979898069620106
-              - 179.99999999999966
-              - 61.606396371386154
+              - -123.75000000000001
+              - -40.97989806962013
+              - 179.99999999999997
+              - 61.60639637138628
           - zoom: 6
             tile_count: 72
             min_tile_size: 64
             max_tile_size: 97
             avg_tile_size: 68.29166666666667
             bbox:
-              - -123.74999999999957
-              - -40.979898069620305
-              - 180.00000000000009
-              - 61.606396371386104
+              - -123.75000000000001
+              - -40.97989806962015
+              - 180.00000000000003
+              - 61.60639637138628
         "###);
 
         Ok(())

diff --git a/tests/expected/martin-cp/flat-with-hash_summary.txt b/tests/expected/martin-cp/flat-with-hash_summary.txt
@@ -7,7 +7,7 @@ Page size: 512B
     1 |         4 |      474B |      983B |      609B | -180,-85,180,85
     2 |         5 |      150B |      865B |      451B | -90,-67,180,67
     3 |         8 |       57B |      839B |      264B | -45,-41,180,67
-    4 |        13 |       57B |      751B |      216B | -22,-22,157,56
+    4 |        13 |       57B |      751B |      216B | -23,-22,158,56
     5 |        27 |       57B |      666B |      167B | -11,-11,146,49
     6 |        69 |       57B |      636B |      127B | -6,-6,146,45
   all |       127 |       57B |      983B |      187B | -180,-85,180,85

diff --git a/tests/expected/martin-cp/flat_summary.txt b/tests/expected/martin-cp/flat_summary.txt
@@ -4,10 +4,10 @@ Page size: 512B
 
  Zoom |   Count   | Smallest  |  Largest  |  Average  | Bounding Box
     0 |         1 |      643B |      643B |      643B | -180,-85,180,85
-    1 |         2 |      150B |      172B |      161B | -180,-85,-0,85
+    1 |         2 |      150B |      172B |      161B | -180,-85,0,85
     2 |         4 |      291B |      690B |      414B | -90,-67,90,67
     3 |         7 |       75B |      727B |      263B | -45,-41,90,67
-    4 |        13 |       75B |      684B |      225B | -22,-22,157,56
+    4 |        13 |       75B |      684B |      225B | -23,-22,158,56
     5 |        27 |       75B |      659B |      195B | -11,-11,146,49
     6 |        69 |       75B |      633B |      155B | -6,-6,146,45
   all |       123 |       75B |      727B |      190B | -180,-85,180,85