From f21119d54895d77f0503dc3fc2dbf06740095686 Mon Sep 17 00:00:00 2001 From: Lucas Date: Mon, 18 Dec 2023 13:02:49 +0800 Subject: [PATCH] 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 --- Cargo.lock | 1 + martin-tile-utils/Cargo.toml | 2 +- martin-tile-utils/src/lib.rs | 163 +++++++++++------- mbtiles/src/summary.rs | 64 +++---- .../martin-cp/flat-with-hash_summary.txt | 2 +- tests/expected/martin-cp/flat_summary.txt | 4 +- 6 files changed, 133 insertions(+), 103 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index 7c9253342..0c6ed0352 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1904,6 +1904,7 @@ name = "martin-tile-utils" version = "0.2.0" dependencies = [ "approx", + "insta", "tile-grid", ] diff --git a/martin-tile-utils/Cargo.toml b/martin-tile-utils/Cargo.toml index 6ed71e7d4..63eca37c3 100644 --- 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 index 2b187e7e2..5e7aa4f0f 100644 --- 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 = 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,12 +188,14 @@ 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. @@ -209,24 +203,31 @@ pub fn tile_index(lon: f64, lat: f64, zoom: u8) -> (u32, u32) { #[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,6 +251,19 @@ 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)] @@ -257,6 +271,7 @@ mod tests { use std::fs::read; use approx::assert_relative_eq; + use insta::assert_snapshot; use Encoding::{Internal, Uncompressed}; use Format::{Jpeg, Json, Png, Webp}; @@ -296,7 +311,7 @@ mod tests { } #[test] - fn test_tile_index() { + fn test_tile_colrow() { assert_eq!((0, 0), tile_index(-180.0, 85.0511, 0)); } @@ -304,53 +319,67 @@ mod tests { 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 index b7eb474da..095a194a7 100644 --- 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 index 25b81ef4d..fedde3215 100644 --- 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 index 538d464ef..992a63bbb 100644 --- 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