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.
2024-12-19 12:51:37 +03:00 · 2023-12-08 23:38:15 -05:00 · 2023-12-08 23:38:15 -05:00 · 280dbe162f
commit 280dbe162f
parent 55fbf8f88a
10 changed files with 106 additions and 85 deletions
--- 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",
--- 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 <YuriAstrakhan@gmail.com>", "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
--- 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);
+    }
 }
--- 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<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;
--- 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()
    };
--- 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"] }
--- 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<MbtTypeCli> = Some(MbtTypeCli::Flat);
    const FLAT_WITH_HASH: Option<MbtTypeCli> = Some(MbtTypeCli::FlatWithHash);
    const NORM_CLI: Option<MbtTypeCli> = Some(MbtTypeCli::Normalized);
--- 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(())
--- 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
--- 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