From 280dbe162f3301d7c8ba20ee8d27500dcb00c21d Mon Sep 17 00:00:00 2001 From: Yuri Astrakhan Date: Fri, 8 Dec 2023 23:38:15 -0500 Subject: [PATCH] Move math utils to martin-tile-utils (#1056) Mathematics should be consolidated in the utils crate. Eventually, I hope we will move it to [tile-grid](https://crates.io/crates/tile-grid) or [utiles](https://crates.io/crates/utiles) crates and consolidate our efforts there. Also, modify earth circumference by a few centimeters (to match various other implementations): from `.7` to `.685_578_5`. This modified a few unit tests. --- Cargo.lock | 6 +- martin-tile-utils/Cargo.toml | 5 +- martin-tile-utils/src/lib.rs | 86 ++++++++++++++++++- martin/src/bin/martin-cp.rs | 23 +---- martin/src/pg/table_source.rs | 3 +- mbtiles/Cargo.toml | 1 - mbtiles/src/copier.rs | 2 + mbtiles/src/summary.rs | 61 ++----------- .../martin-cp/flat-with-hash_summary.txt | 2 +- tests/expected/martin-cp/flat_summary.txt | 2 +- 10 files changed, 106 insertions(+), 85 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index c0b67f6f..d73132cd 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1894,7 +1894,10 @@ dependencies = [ [[package]] name = "martin-tile-utils" -version = "0.1.5" +version = "0.1.6" +dependencies = [ + "approx", +] [[package]] name = "mbtiles" @@ -1902,7 +1905,6 @@ version = "0.8.3" dependencies = [ "actix-rt", "anyhow", - "approx", "clap", "ctor", "enum-display", diff --git a/martin-tile-utils/Cargo.toml b/martin-tile-utils/Cargo.toml index deb38938..f8289357 100644 --- a/martin-tile-utils/Cargo.toml +++ b/martin-tile-utils/Cargo.toml @@ -2,7 +2,7 @@ lints.workspace = true [package] name = "martin-tile-utils" -version = "0.1.5" +version = "0.1.6" authors = ["Yuri Astrakhan ", "MapLibre contributors"] description = "Utilites to help with map tile processing, such as type and compression detection. Used by the MapLibre's Martin tile server." keywords = ["maps", "tiles", "mvt", "tileserver"] @@ -17,3 +17,6 @@ repository.workspace = true rust-version.workspace = true [dependencies] + +[dev-dependencies] +approx.workspace = true diff --git a/martin-tile-utils/src/lib.rs b/martin-tile-utils/src/lib.rs index da18e284..3c9b31a7 100644 --- a/martin-tile-utils/src/lib.rs +++ b/martin-tile-utils/src/lib.rs @@ -6,7 +6,7 @@ use std::f64::consts::PI; use std::fmt::Display; -pub const EARTH_CIRCUMFERENCE: f64 = 40_075_016.7; +pub const EARTH_CIRCUMFERENCE: f64 = 40_075_016.685_578_5; pub const EARTH_RADIUS: f64 = EARTH_CIRCUMFERENCE / 2.0 / PI; #[derive(Clone, Copy, Debug, PartialEq, Eq)] @@ -184,10 +184,69 @@ 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)) +} + +/// 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, + ); + + [min_lng, min_lat, max_lng, max_lat] +} + +/// 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) +} + +/// Compute precision of a zoom level, i.e. how many decimal digits of the longitude and latitude are relevant +#[must_use] +#[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)] +pub fn get_zoom_precision(zoom: u8) -> usize { + let lng_delta = webmercator_to_wgs84(EARTH_CIRCUMFERENCE / f64::from(1_u32 << zoom), 0.0).0; + let log = lng_delta.log10() - 0.5; + if log > 0.0 { + 0 + } else { + -log.ceil() as usize + } +} + +#[must_use] +pub fn webmercator_to_wgs84(x: f64, y: f64) -> (f64, f64) { + let lng = (x / EARTH_RADIUS).to_degrees(); + let lat = f64::atan(f64::sinh(y / EARTH_RADIUS)).to_degrees(); + (lng, lat) +} + #[cfg(test)] mod tests { use std::fs::read; + use approx::assert_relative_eq; use Encoding::{Internal, Uncompressed}; use Format::{Jpeg, Json, Png, Webp}; @@ -225,4 +284,29 @@ mod tests { info(Json, Uncompressed) ); } + + #[test] + fn test_tile_index() { + assert_eq!((0, 0), tile_index(-180.0, 85.0511, 0)); + } + + #[test] + #[allow(clippy::unreadable_literal)] + fn meter_to_lng_lat() { + let (lng, lat) = webmercator_to_wgs84(-20037508.34, -20037508.34); + assert_relative_eq!(lng, -179.9999999749437, epsilon = f64::EPSILON * 2.0); + assert_relative_eq!(lat, -85.05112877764508, epsilon = f64::EPSILON * 2.0); + + let (lng, lat) = webmercator_to_wgs84(20037508.34, 20037508.34); + assert_relative_eq!(lng, 179.9999999749437, epsilon = f64::EPSILON * 2.0); + assert_relative_eq!(lat, 85.05112877764508, epsilon = f64::EPSILON * 2.0); + + let (lng, lat) = webmercator_to_wgs84(0.0, 0.0); + assert_relative_eq!(lng, 0.0, epsilon = f64::EPSILON * 2.0); + assert_relative_eq!(lat, 0.0, epsilon = f64::EPSILON * 2.0); + + let (lng, lat) = webmercator_to_wgs84(3000.0, 9000.0); + assert_relative_eq!(lng, 0.026949458523585632, epsilon = f64::EPSILON * 2.0); + assert_relative_eq!(lat, 0.08084834874097367, epsilon = f64::EPSILON * 2.0); + } } diff --git a/martin/src/bin/martin-cp.rs b/martin/src/bin/martin-cp.rs index e6402320..5a29a8d1 100644 --- a/martin/src/bin/martin-cp.rs +++ b/martin/src/bin/martin-cp.rs @@ -1,4 +1,3 @@ -use std::f64::consts::PI; use std::fmt::{Debug, Display, Formatter}; use std::path::PathBuf; use std::sync::atomic::{AtomicU64, Ordering}; @@ -17,7 +16,7 @@ use martin::{ append_rect, read_config, Config, IdResolver, MartinError, MartinResult, ServerState, Source, TileCoord, TileData, TileRect, }; -use martin_tile_utils::TileInfo; +use martin_tile_utils::{bbox_to_xyz, TileInfo}; use mbtiles::sqlx::SqliteConnection; use mbtiles::{ init_mbtiles_schema, is_empty_database, CopyDuplicateMode, MbtType, MbtTypeCli, Mbtiles, @@ -151,17 +150,6 @@ async fn start(copy_args: CopierArgs) -> MartinCpResult<()> { run_tile_copy(copy_args.copy, sources).await } -/// Convert longitude and latitude to tile index -#[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)] -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)) -} - fn compute_tile_ranges(args: &CopyArgs) -> Vec { let mut ranges = Vec::new(); let mut zooms_vec = Vec::new(); @@ -179,8 +167,8 @@ fn compute_tile_ranges(args: &CopyArgs) -> Vec { }; for zoom in zooms { for bbox in &boxes { - let (min_x, min_y) = tile_index(bbox.left, bbox.top, *zoom); - let (max_x, max_y) = tile_index(bbox.right, bbox.bottom, *zoom); + let (min_x, min_y, max_x, max_y) = + 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), @@ -428,11 +416,6 @@ mod tests { use super::*; - #[test] - fn test_tile_index() { - assert_eq!((0, 0), tile_index(-180.0, 85.0511, 0)); - } - #[test] fn test_compute_tile_ranges() { let world = Bounds::MAX_TILED; diff --git a/martin/src/pg/table_source.rs b/martin/src/pg/table_source.rs index b0ed0ec1..16d307a2 100644 --- a/martin/src/pg/table_source.rs +++ b/martin/src/pg/table_source.rs @@ -159,8 +159,7 @@ pub async fn table_to_query( // TODO: we should use ST_Expand here, but it may require a bit more math work, // so might not be worth it as it is only used for PostGIS < v3.1. // v3.1 has been out for 2+ years (december 2020) - // let earth_circumference = 40075016.6855785; - // let val = earth_circumference * buffer as f64 / extent as f64; + // let val = EARTH_CIRCUMFERENCE * buffer as f64 / extent as f64; // format!("ST_Expand(ST_TileEnvelope($1::integer, $2::integer, $3::integer), {val}/2^$1::integer)") "ST_TileEnvelope($1::integer, $2::integer, $3::integer)".to_string() }; diff --git a/mbtiles/Cargo.toml b/mbtiles/Cargo.toml index a7fd1ae6..cb5c4f76 100644 --- a/mbtiles/Cargo.toml +++ b/mbtiles/Cargo.toml @@ -40,7 +40,6 @@ tokio = { workspace = true, features = ["rt-multi-thread"], optional = true } [dev-dependencies] # For testing, might as well use the same async framework as the Martin itself actix-rt.workspace = true -approx.workspace = true ctor.workspace = true env_logger.workspace = true insta = { workspace = true, features = ["toml", "yaml"] } diff --git a/mbtiles/src/copier.rs b/mbtiles/src/copier.rs index 6640cd19..bcfa4a18 100644 --- a/mbtiles/src/copier.rs +++ b/mbtiles/src/copier.rs @@ -613,6 +613,8 @@ mod tests { use super::*; + // TODO: Most of these tests are duplicating the tests from tests/mbtiles.rs, and should be cleaned up/removed. + const FLAT: Option = Some(MbtTypeCli::Flat); const FLAT_WITH_HASH: Option = Some(MbtTypeCli::FlatWithHash); const NORM_CLI: Option = Some(MbtTypeCli::Normalized); diff --git a/mbtiles/src/summary.rs b/mbtiles/src/summary.rs index 5ea82e3f..ae4298f5 100644 --- a/mbtiles/src/summary.rs +++ b/mbtiles/src/summary.rs @@ -8,7 +8,7 @@ use std::fmt::{Display, Formatter}; use std::path::PathBuf; use std::str::FromStr; -use martin_tile_utils::{EARTH_CIRCUMFERENCE, EARTH_RADIUS}; +use martin_tile_utils::{get_zoom_precision, xyz_to_bbox}; use serde::Serialize; use size_format::SizeFormatterBinary; use sqlx::{query, SqliteExecutor}; @@ -158,7 +158,8 @@ impl Mbtiles { r.min_tile_y.unwrap(), r.max_tile_x.unwrap(), r.max_tile_y.unwrap(), - ), + ) + .into(), } }) .collect(); @@ -186,66 +187,14 @@ impl Mbtiles { } } -/// Convert min/max XYZ tile coordinates to a bounding box -fn xyz_to_bbox(zoom: u8, min_x: i32, min_y: i32, max_x: i32, max_y: i32) -> Bounds { - 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, - ); - - Bounds::new(min_lng, min_lat, max_lng, max_lat) -} - -fn get_zoom_precision(zoom: u8) -> usize { - let lng_delta = webmercator_to_wgs84(EARTH_CIRCUMFERENCE / f64::from(1_u32 << zoom), 0.0).0; - let log = lng_delta.log10() - 0.5; - if log > 0.0 { - 0 - } else { - -log.ceil() as usize - } -} - -fn webmercator_to_wgs84(x: f64, y: f64) -> (f64, f64) { - let lng = (x / EARTH_RADIUS).to_degrees(); - let lat = (f64::atan(f64::sinh(y / EARTH_RADIUS))).to_degrees(); - (lng, lat) -} - #[cfg(test)] mod tests { #![allow(clippy::unreadable_literal)] - use approx::assert_relative_eq; use insta::assert_yaml_snapshot; - use crate::summary::webmercator_to_wgs84; use crate::{init_mbtiles_schema, MbtResult, MbtType, Mbtiles}; - #[actix_rt::test] - async fn meter_to_lng_lat() { - let (lng, lat) = webmercator_to_wgs84(-20037508.34, -20037508.34); - assert_relative_eq!(lng, -179.99999991016847, epsilon = f64::EPSILON); - assert_relative_eq!(lat, -85.05112877205713, epsilon = f64::EPSILON); - - let (lng, lat) = webmercator_to_wgs84(20037508.34, 20037508.34); - assert_relative_eq!(lng, 179.99999991016847, epsilon = f64::EPSILON); - assert_relative_eq!(lat, 85.05112877205713, epsilon = f64::EPSILON); - - let (lng, lat) = webmercator_to_wgs84(0.0, 0.0); - assert_relative_eq!(lng, 0.0, epsilon = f64::EPSILON); - assert_relative_eq!(lat, 0.0, epsilon = f64::EPSILON); - - let (lng, lat) = webmercator_to_wgs84(3000.0, 9000.0); - assert_relative_eq!(lng, 0.02694945851388753, epsilon = f64::EPSILON); - assert_relative_eq!(lat, 0.0808483487118794, epsilon = f64::EPSILON); - } - #[actix_rt::test] async fn summary_empty_file() -> MbtResult<()> { let mbt = Mbtiles::new("file:mbtiles_empty_summary?mode=memory&cache=shared")?; @@ -356,7 +305,7 @@ mod tests { - -123.75000000000001 - -40.97989806962013 - 180 - - 61.60639637138627 + - 61.60639637138628 - zoom: 6 tile_count: 72 min_tile_size: 64 @@ -366,7 +315,7 @@ mod tests { - -123.75000000000001 - -40.97989806962013 - 180 - - 61.60639637138627 + - 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 25b81ef4..fedde321 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 6ecb5f4e..992a63bb 100644 --- a/tests/expected/martin-cp/flat_summary.txt +++ b/tests/expected/martin-cp/flat_summary.txt @@ -7,7 +7,7 @@ Page size: 512B 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