Move math utils to martin-tile-utils

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 <[email protected]>", "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

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);
+    }
 }

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<TileRect> {
     let mut ranges = Vec::new();
     let mut zooms_vec = Vec::new();
@@ -179,8 +167,8 @@ fn compute_tile_ranges(args: &CopyArgs) -> Vec<TileRect> {
     };
     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;

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()
     };

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"] }

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<MbtTypeCli> = Some(MbtTypeCli::Flat);
     const FLAT_WITH_HASH: Option<MbtTypeCli> = Some(MbtTypeCli::FlatWithHash);
     const NORM_CLI: Option<MbtTypeCli> = Some(MbtTypeCli::Normalized);

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(())

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

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