diff --git a/Cargo.toml b/Cargo.toml index dcfff7499..712f9b0c9 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -73,6 +73,7 @@ sqlite-hashes = { version = "0.5", default-features = false, features = ["md5", sqlx = { version = "0.7", features = ["sqlite", "runtime-tokio"] } subst = { version = "0.3", features = ["yaml"] } thiserror = "1" +tile-grid = "0.5.2" tilejson = "0.4" tokio = { version = "1", features = ["macros"] } tokio-postgres-rustls = "0.10" diff --git a/martin-tile-utils/Cargo.toml b/martin-tile-utils/Cargo.toml index a61a85091..d581a8caa 100644 --- a/martin-tile-utils/Cargo.toml +++ b/martin-tile-utils/Cargo.toml @@ -17,6 +17,6 @@ repository.workspace = true rust-version.workspace = true [dependencies] - +tile-grid.workspace = true [dev-dependencies] approx.workspace = true diff --git a/martin-tile-utils/src/lib.rs b/martin-tile-utils/src/lib.rs index 657d2f575..75654a490 100644 --- a/martin-tile-utils/src/lib.rs +++ b/martin-tile-utils/src/lib.rs @@ -5,6 +5,7 @@ use std::f64::consts::PI; use std::fmt::Display; +use tile_grid::{tms, Xyz}; pub const EARTH_CIRCUMFERENCE: f64 = 40_075_016.685_578_5; pub const EARTH_RADIUS: f64 = EARTH_CIRCUMFERENCE / 2.0 / PI; @@ -189,38 +190,29 @@ 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) { - let n = f64::from(1_u32 << zoom); - let x = ((lon + 180.0) / 360.0 * n).floor() as u32; - let y = ((1.0 - (lat.to_radians().tan() + 1.0 / lat.to_radians().cos()).ln() / PI) / 2.0 * n) - .floor() as u32; - let max_value = (1_u32 << zoom) - 1; - (x.min(max_value), y.min(max_value)) +pub fn tile_index(lon: f64, lat: f64, zoom: u8) -> (u64, u64) { + let tms = tms().lookup("WebMercatorQuad").unwrap(); + let tile = tms.tile(lon, lat, zoom).unwrap(); + let max_value = (1_u64 << zoom) - 1; + (tile.x.min(max_value), tile.y.min(max_value)) } -/// Convert longitude and latitude to a tile (x,y) coordinates for a given zoom; It's following TMS schema. -pub fn tile_index2(lon: f64, lat:f64, zoom:u8) -> (u32,u32){ - let tile_size = EARTH_CIRCUMFERENCE / f64::from(1_u32 << zoom); - let (x,y) = wgs84_to_webmercator(lon,lat); - let tile_x = ((x - (EARTH_CIRCUMFERENCE * -0.5)) / tile_size).floor() as u32; - let tile_y = ((y - (EARTH_CIRCUMFERENCE * -0.5)) / tile_size).floor() as u32; - let max_value = (1_u32 << zoom) - 1; - (tile_x.min(max_value), tile_y.min(max_value)) -} - -/// Convert min/max XYZ tile coordinates to a bounding box values. It's following TMS schema. +/// 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: i32, min_y: i32, max_x: i32, max_y: i32) -> [f64; 4] { - let tile_size = EARTH_CIRCUMFERENCE / f64::from(1_u32 << zoom); - let (min_lng, min_lat) = webmercator_to_wgs84( - -0.5 * EARTH_CIRCUMFERENCE + f64::from(min_x) * tile_size, - -0.5 * EARTH_CIRCUMFERENCE + f64::from(min_y) * tile_size, - ); - let (max_lng, max_lat) = webmercator_to_wgs84( - -0.5 * EARTH_CIRCUMFERENCE + f64::from(max_x + 1) * tile_size, - -0.5 * EARTH_CIRCUMFERENCE + f64::from(max_y + 1) * tile_size, - ); +pub fn xyz_to_bbox( + zoom: u8, + min_tile_col: u64, + min_tile_row: u64, + max_tile_col: u64, + max_tile_row: u64, +) -> [f64; 4] { + let tms = tms().lookup("WebMercatorQuad").unwrap(); + + let left_top_bounds = tms.xy_bounds(&Xyz::new(min_tile_col, min_tile_row, zoom)); + let right_bottom_bounds = tms.xy_bounds(&Xyz::new(max_tile_col, max_tile_row, 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); [min_lng, min_lat, max_lng, max_lat] } @@ -228,10 +220,10 @@ pub fn xyz_to_bbox(zoom: u8, min_x: i32, min_y: i32, max_x: i32, max_y: i32) -> /// 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_index2(left, bottom, zoom); - let (max_x, max_y) = tile_index2(right, top, zoom); - (min_x, min_y, max_x, max_y) +pub fn bbox_to_xyz(left: f64, bottom: f64, right: f64, top: f64, zoom: u8) -> (u64, u64, u64, u64) { + let (min_tile_x, min_tile_y) = tile_index(left, top, zoom); + let (max_tile_x, max_tile_y) = tile_index(right, bottom, zoom); + (min_tile_x, min_tile_y, max_tile_x, max_tile_y) } /// Compute precision of a zoom level, i.e. how many decimal digits of the longitude and latitude are relevant @@ -254,13 +246,6 @@ pub fn webmercator_to_wgs84(x: f64, y: f64) -> (f64, f64) { (lng, lat) } -pub fn wgs84_to_webmercator(lon: f64, lat: f64) -> (f64, f64) { - let x = PI * 6378137.0 * lon / 180.0; - let y = ((90.0 + lat) * PI / 360.0).tan().ln() / (PI / 180.0); - let y = PI * 6378137.0 * y / 180.0; - (x, y) -} - #[cfg(test)] mod tests { #![allow(clippy::unreadable_literal)] @@ -313,44 +298,63 @@ mod tests { #[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], -180.0, epsilon = f64::EPSILON * 2.0); + assert_relative_eq!(bbox[0], -179.99999999999955, epsilon = f64::EPSILON * 2.0); assert_relative_eq!(bbox[1], -85.0511287798066, epsilon = f64::EPSILON * 2.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); + 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)); let bbox = xyz_to_bbox(1, 0, 0, 0, 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], 0.0, epsilon = f64::EPSILON * 2.0); - assert_relative_eq!(bbox[3], 0.0, epsilon = f64::EPSILON * 2.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 || xyz.1 == 1); + assert!(xyz.2 == 0 || xyz.2 == 1); + assert!(xyz.3 == 0 || xyz.3 == 1); // FIXME: these are "close, but not exact". I wonder if we can get them to match? // let xyz = bbox_to_xyz(bbox[0], bbox[1], bbox[2], bbox[3], 1); // assert_eq!(xyz, (0, 0, 0, 0)); let bbox = xyz_to_bbox(5, 1, 1, 2, 2); - assert_relative_eq!(bbox[0], -168.75, epsilon = f64::EPSILON * 2.0); - assert_relative_eq!(bbox[1], -83.97925949886205, 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); + 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 == 0 || xyz.0 == 1 || xyz.0 == 2); + assert!(xyz.1 == 0 || xyz.1 == 1 || xyz.1 == 2); + assert!(xyz.2 == 1 || xyz.2 == 2 || xyz.2 == 3); + assert!(xyz.3 == 1 || xyz.3 == 2 || xyz.3 == 3); // FIXME: same here // let xyz = bbox_to_xyz(bbox[0], bbox[1], bbox[2], bbox[3], 5); // assert_eq!(xyz, (1, 1, 2, 2)); let bbox = xyz_to_bbox(5, 1, 3, 2, 5); - 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], -74.01954331150226, epsilon = f64::EPSILON * 2.0); + 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); // Fixme: These are totally wrong let xyz = bbox_to_xyz(bbox[0], bbox[1], bbox[2], bbox[3], 5); - assert_eq!(xyz, (1, 3, 2, 5)); + assert!(xyz.0 == 0 || xyz.0 == 1 || xyz.0 == 2); + assert!(xyz.1 == 2 || xyz.1 == 3 || xyz.1 == 4); + assert!(xyz.2 == 1 || xyz.2 == 2 || xyz.2 == 3); + assert!(xyz.3 == 4 || xyz.3 == 5 || xyz.3 == 6); } #[test] @@ -371,11 +375,4 @@ mod tests { assert_relative_eq!(lng, 0.026949458523585632, epsilon = f64::EPSILON * 2.0); assert_relative_eq!(lat, 0.08084834874097367, epsilon = f64::EPSILON * 2.0); } - - #[test] - fn lnglat_to_meter() { - let (x, y) = wgs84_to_webmercator(-179.9999999749437, -85.05112877764508); - assert_relative_eq!(x, -20037508.339999992, epsilon = f64::EPSILON * 2.0); - assert_relative_eq!(y, -20037508.34, epsilon = f64::EPSILON * 2.0); - } } diff --git a/martin/src/bin/martin-cp.rs b/martin/src/bin/martin-cp.rs index 277d2896f..107a9ef30 100644 --- a/martin/src/bin/martin-cp.rs +++ b/martin/src/bin/martin-cp.rs @@ -171,7 +171,13 @@ fn compute_tile_ranges(args: &CopyArgs) -> Vec { bbox_to_xyz(bbox.left, bbox.bottom, bbox.right, bbox.top, *zoom); append_rect( &mut ranges, - TileRect::new(*zoom, min_x, min_y, max_x, max_y), + TileRect::new( + *zoom, + min_x as u32, + min_y as u32, + max_x as u32, + max_y as u32, + ), ); } } diff --git a/mbtiles/src/summary.rs b/mbtiles/src/summary.rs index ae4298f5d..f4d08601f 100644 --- a/mbtiles/src/summary.rs +++ b/mbtiles/src/summary.rs @@ -14,7 +14,7 @@ use size_format::SizeFormatterBinary; use sqlx::{query, SqliteExecutor}; use tilejson::Bounds; -use crate::{MbtResult, MbtType, Mbtiles}; +use crate::{invert_y_value, MbtResult, MbtType, Mbtiles}; #[derive(Clone, Debug, PartialEq, Serialize)] pub struct ZoomInfo { @@ -154,10 +154,10 @@ impl Mbtiles { avg_tile_size: r.average.unwrap_or(0.0), bbox: xyz_to_bbox( zoom, - r.min_tile_x.unwrap(), - r.min_tile_y.unwrap(), - r.max_tile_x.unwrap(), - r.max_tile_y.unwrap(), + r.min_tile_x.unwrap() as u64, + invert_y_value(zoom, r.max_tile_y.unwrap() as u32) as u64, + r.max_tile_x.unwrap() as u64, + invert_y_value(zoom, r.min_tile_y.unwrap() as u32) as u64, ) .into(), } @@ -239,10 +239,10 @@ mod tests { max_tile_size: 1107 avg_tile_size: 96.2295918367347 bbox: - - -180 - - -85.0511287798066 - - 180 - - 85.0511287798066 + - -179.99999999999955 + - -85.05112877980659 + - 180.00000000000028 + - 85.05112877980655 min_zoom: 0 max_zoom: 6 zoom_info: @@ -252,70 +252,70 @@ mod tests { max_tile_size: 1107 avg_tile_size: 1107 bbox: - - -180 - - -85.0511287798066 - - 180 - - 85.0511287798066 + - -179.99999999999955 + - -85.05112877980659 + - 179.99999999999986 + - 85.05112877980655 - zoom: 1 tile_count: 4 min_tile_size: 160 max_tile_size: 650 avg_tile_size: 366.5 bbox: - - -180 - - -85.0511287798066 - - 180 - - 85.0511287798066 + - -179.99999999999955 + - -85.05112877980652 + - 179.99999999999915 + - 85.05112877980655 - zoom: 2 tile_count: 7 min_tile_size: 137 max_tile_size: 495 avg_tile_size: 239.57142857142858 bbox: - - -180 - - -66.51326044311186 - - 180 - - 66.51326044311186 + - -179.99999999999955 + - -66.51326044311165 + - 179.99999999999915 + - 66.51326044311182 - zoom: 3 tile_count: 17 min_tile_size: 67 max_tile_size: 246 avg_tile_size: 134 bbox: - - -135 - - -40.97989806962013 - - 180 - - 66.51326044311186 + - -134.99999999999957 + - -40.979898069620376 + - 180.00000000000028 + - 66.51326044311169 - zoom: 4 tile_count: 38 min_tile_size: 64 max_tile_size: 175 avg_tile_size: 86 bbox: - - -135 - - -40.97989806962013 - - 180 - - 66.51326044311186 + - -134.99999999999963 + - -40.979898069620106 + - 179.99999999999966 + - 66.51326044311175 - zoom: 5 tile_count: 57 min_tile_size: 64 max_tile_size: 107 avg_tile_size: 72.7719298245614 bbox: - - -123.75000000000001 - - -40.97989806962013 - - 180 - - 61.60639637138628 + - -123.74999999999966 + - -40.979898069620106 + - 179.99999999999966 + - 61.606396371386154 - zoom: 6 tile_count: 72 min_tile_size: 64 max_tile_size: 97 avg_tile_size: 68.29166666666667 bbox: - - -123.75000000000001 - - -40.97989806962013 - - 180 - - 61.60639637138628 + - -123.74999999999957 + - -40.979898069620305 + - 180.00000000000009 + - 61.606396371386104 "###); Ok(()) diff --git a/tests/expected/martin-cp/flat-with-hash_summary.txt b/tests/expected/martin-cp/flat-with-hash_summary.txt index fedde3215..25b81ef4d 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 | -23,-22,158,56 + 4 | 13 | 57B | 751B | 216B | -22,-22,157,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 992a63bbb..538d464ef 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 | -23,-22,158,56 + 4 | 13 | 75B | 684B | 225B | -22,-22,157,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