Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix bbox_to_xyz #1070

Merged
merged 20 commits into from
Dec 18, 2023
Merged

Fix bbox_to_xyz #1070

Show file tree
Hide file tree
Changes from 17 commits
Commits
Show all changes
20 commits
Select commit Hold shift + click to select a range
ebe8303
Add wgs84_to_webmercator to tile util
sharkAndshark Dec 15, 2023
d140f4e
Add tile_colrow to tile-util
sharkAndshark Dec 15, 2023
2581bff
Add bbox_to_colrow to tile-util
sharkAndshark Dec 15, 2023
e33c787
Use bbox_to_colrow
sharkAndshark Dec 15, 2023
1ba2d92
Update test
sharkAndshark Dec 15, 2023
c3b4c42
Rename bbox_to_colrow to bbox_to_xyz and use u32
sharkAndshark Dec 15, 2023
ae8b98d
fmt and clippy
sharkAndshark Dec 15, 2023
56846a8
clippy
sharkAndshark Dec 15, 2023
7945f2c
clippy
sharkAndshark Dec 15, 2023
2dfff7d
Update martin-tile-utils/src/lib.rs
nyurik Dec 16, 2023
1bb848e
Update martin-tile-utils/src/lib.rs
nyurik Dec 16, 2023
2d69b1a
a few minor fixes
nyurik Dec 16, 2023
2719626
Remove tile_index
sharkAndshark Dec 17, 2023
0fd5e16
Minior imrovement to tile_colrow function
sharkAndshark Dec 17, 2023
4bc205b
remove tile-grid
sharkAndshark Dec 17, 2023
8bde903
Remove tile-grid, udpate mercator and udpate test
sharkAndshark Dec 17, 2023
8f9c619
Add comment
sharkAndshark Dec 17, 2023
93d9f27
Rename arguments and remove tiny functions
sharkAndshark Dec 18, 2023
a83ccd7
Update test results
sharkAndshark Dec 18, 2023
eea4a86
Merge remote-tracking branch 'maplibre/main' into tilecolrow
sharkAndshark Dec 18, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions Cargo.lock
View file
Open in desktop

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

2 changes: 1 addition & 1 deletion martin-tile-utils/Cargo.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ repository.workspace = true
rust-version.workspace = true

[dependencies]
tile-grid.workspace = true

[dev-dependencies]
approx.workspace = true
insta.workspace = true
167 changes: 101 additions & 66 deletions martin-tile-utils/src/lib.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 = 30;
use std::sync::OnceLock;

fn web_merc() -> &'static Tms {
static TMS: OnceLock<Tms> = OnceLock::new();
TMS.get_or_init(|| tms().lookup("WebMercatorQuad").unwrap())
}

#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum Format {
Expand Down Expand Up @@ -196,35 +188,46 @@ 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_colrow(lng: f64, lat: f64, zoom: u8) -> (u32, u32) {
sharkAndshark marked this conversation as resolved.
Show resolved Hide resolved
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.
/// The result is `[min_lng, min_lat, max_lng, max_lat]`
#[must_use]
pub fn xyz_to_bbox(zoom: u8, min_x: u32, min_y: u32, max_x: u32, max_y: u32) -> [f64; 4] {
pub fn xyz_to_bbox(zoom: u8, min_col: u32, min_row: u32, max_col: u32, max_row: u32) -> [f64; 4] {
sharkAndshark marked this conversation as resolved.
Show resolved Hide resolved
assert!(zoom <= MAX_ZOOM, "zoom {zoom} must be <= {MAX_ZOOM}");
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_col, max_row, tile_length);
let right_top_bbox = tile_bbox(max_col, min_row, 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_colrow(left, top, zoom);
let (max_col, max_row) = tile_colrow(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
Expand All @@ -248,13 +251,31 @@ 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
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's from esri.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this code is actually from the OSM wiki that has all this math

#[must_use]
pub fn wgs84_to_webmercator(lon: f64, lat: f64) -> (f64, f64) {
(to_mercator_x(lon), to_mercator_y(lat))
}

fn to_mercator_y(lat: f64) -> f64 {
sharkAndshark marked this conversation as resolved.
Show resolved Hide resolved
let rad = lat * PI / 180.0;
let sin = rad.sin();
EARTH_RADIUS / 2.0 * ((1.0 + sin) / (1.0 - sin)).ln()
}

fn to_mercator_x(lon: f64) -> f64 {
lon * PI / 180.0 * EARTH_RADIUS
}

#[cfg(test)]
mod tests {
#![allow(clippy::unreadable_literal)]

use std::fs::read;

use approx::assert_relative_eq;
use insta::assert_snapshot;
use Encoding::{Internal, Uncompressed};
use Format::{Jpeg, Json, Png, Webp};

Expand Down Expand Up @@ -294,61 +315,75 @@ mod tests {
}

#[test]
fn test_tile_index() {
assert_eq!((0, 0), tile_index(-180.0, 85.0511, 0));
fn test_tile_colrow() {
assert_eq!((0, 0), tile_colrow(-180.0, 85.0511, 0));
}

#[test]
fn test_xyz_to_bbox() {
Copy link
Collaborator Author

@sharkAndshark sharkAndshark Dec 17, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Due to the inherent precision limitations of floating point numbers, performing back-and-forth transformations between XYZ and bounding box coordinates within a single unit test seems not a good idea?
I removed the transform back test code. And actually if you add them back, you can see some failed. Any idea? I can't see logic error in these methods. Is it our const not accurate. My mind get lost now.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Or maybe we should reserve tile-grid in this?

// 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]
Expand Down