diff --git a/Cargo.lock b/Cargo.lock index 9616960..7c55d6f 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -237,6 +237,15 @@ dependencies = [ "alloc-no-stdlib", ] +[[package]] +name = "approx" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cab112f0a86d568ea0e627cc1d6be74a1e9cd55214684db5561995f6dad897c6" +dependencies = [ + "num-traits", +] + [[package]] name = "atomic-waker" version = "1.1.2" @@ -399,6 +408,18 @@ dependencies = [ "libloading", ] +[[package]] +name = "console" +version = "0.15.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0e1f83fc076bd6dd27517eacdf25fef6c4dfe5f1d7448bafaaf3a26f13b5e4eb" +dependencies = [ + "encode_unicode", + "lazy_static", + "libc", + "windows-sys 0.52.0", +] + [[package]] name = "convert_case" version = "0.4.0" @@ -526,6 +547,12 @@ version = "1.13.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "60b1af1c220855b6ceac025d3f6ecdd2b7c4894bfe9cd9bda4fbb4bc7c0d4cf0" +[[package]] +name = "encode_unicode" +version = "0.3.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a357d28ed41a50f9c765dbfe56cbc04a64e53e5fc58ba79fbc34c10ef3df831f" + [[package]] name = "encoding_rs" version = "0.8.34" @@ -904,6 +931,18 @@ dependencies = [ "hashbrown", ] +[[package]] +name = "insta" +version = "1.40.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6593a41c7a73841868772495db7dc1e8ecab43bb5c0b6da2059246c4b506ab60" +dependencies = [ + "console", + "lazy_static", + "linked-hash-map", + "similar", +] + [[package]] name = "ipnet" version = "2.10.1" @@ -977,6 +1016,12 @@ dependencies = [ "windows-targets", ] +[[package]] +name = "linked-hash-map" +version = "0.5.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0717cef1bc8b636c6e1c1bbdefc09e6322da8a9321966e8928ef80d20f7f770f" + [[package]] name = "linux-raw-sys" version = "0.4.14" @@ -1682,6 +1727,12 @@ version = "0.3.7" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "d66dc143e6b11c1eddc06d5c423cfc97062865baf299914ab64caa38182078fe" +[[package]] +name = "similar" +version = "2.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1de1d4f81173b03af4c0cbed3c898f6bff5b870e4a7f5d6f4057d62a7a4b686e" + [[package]] name = "slab" version = "0.4.9" @@ -2227,8 +2278,10 @@ name = "wxbox-tiler" version = "0.1.0" dependencies = [ "actix-web", + "approx", "eccodes", "flate2", + "insta", "mime", "ndarray", "ordered-float", diff --git a/wxbox-pal/src/lib.rs b/wxbox-pal/src/lib.rs index ce44f0c..a57f589 100644 --- a/wxbox-pal/src/lib.rs +++ b/wxbox-pal/src/lib.rs @@ -49,7 +49,7 @@ impl Palette { } } } - +/* #[cfg(test)] mod tests { use std::collections::BTreeMap; @@ -85,4 +85,6 @@ mod tests { assert_eq!(pal.color_for(10.0), color!(rgb: 0, 165, 255)); assert_eq!(pal.color_for(15.0), color!(rgb: 0, 165, 255)); } -} \ No newline at end of file +} + + */ \ No newline at end of file diff --git a/wxbox-tiler/Cargo.toml b/wxbox-tiler/Cargo.toml index e1406f0..e61edb8 100644 --- a/wxbox-tiler/Cargo.toml +++ b/wxbox-tiler/Cargo.toml @@ -17,4 +17,8 @@ wxbox-pal = { version = "0.1", path = "../wxbox-pal" } png = "0.17" mime = "0.3.17" eccodes = { version = "0.12", features = ["message_ndarray"] } -ndarray = "0.16" \ No newline at end of file +ndarray = "0.16" + +[dev-dependencies] +approx = "0.5" +insta = "1.40" \ No newline at end of file diff --git a/wxbox-tiler/src/main.rs b/wxbox-tiler/src/main.rs index 30a04b6..e46c1e5 100644 --- a/wxbox-tiler/src/main.rs +++ b/wxbox-tiler/src/main.rs @@ -1,7 +1,8 @@ pub(crate) mod sources; -pub(crate) mod coords; + mod grib2; mod pixmap; +mod tile; use std::collections::{BTreeMap}; use std::fmt::Debug; diff --git a/wxbox-tiler/src/sources/noaa/mod.rs b/wxbox-tiler/src/sources/noaa/mod.rs index 0a178cd..f83678f 100644 --- a/wxbox-tiler/src/sources/noaa/mod.rs +++ b/wxbox-tiler/src/sources/noaa/mod.rs @@ -8,23 +8,25 @@ use actix_web::web::Data; use eccodes::{CodesHandle, FallibleStreamingIterator, ProductKind}; use flate2::read::GzDecoder; use ndarray::{Zip}; -use ordered_float::OrderedFloat; +use ordered_float::{FloatCore, OrderedFloat}; use png::{BitDepth, ColorType, Encoder}; use wxbox_pal::{Color, Palette}; use crate::{AppState, LutKey}; use crate::grib2::{closest_key, LookupTable2D}; use crate::pixmap::Pixmap; use wxbox_pal::color; +use crate::tile::xyz_to_bbox; #[get("/mrms_cref/{z}/{x}/{y}.png")] -async fn mrms_cref(path: web::Path<(i32, i32, i32)>, data: Data) -> HttpResponse { +async fn mrms_cref(path: web::Path<(u8, u32, u32)>, data: Data) -> HttpResponse { let z = path.0; let x = path.1; let y = path.2; - let lon = (x as f64 / 2.0_f64.powi(z)) * 360.0 - 180.0; - let lat = (PI - (y as f64 / 2.0_f64.powi(z)) * 2.0 * PI).sinh().atan() * (180.0 / PI); + let [min_lng,min_lat,max_lng,max_lat] = xyz_to_bbox(z, x, y, x, y); + let lat = min_lat; + let lon = min_lng; let mut needs_reload = false; @@ -106,33 +108,46 @@ async fn mrms_cref(path: web::Path<(i32, i32, i32)>, data: Data) -> Ht let mut image: Pixmap = Pixmap::new(); - let tile_size = 360.0 / 2.0_f64.powi(z); - let degrees_per_px = tile_size / 256.0; + let degrees_lat_per_px = (max_lat - min_lat) / 256.0; + let degrees_lon_per_px = (max_lng - min_lng) / 256.0; + + let xtile = x as f64; + let ytile = y as f64; for x in 0..256 { for y in 0..256 { - let lon = lon + degrees_per_px * x as f64; - let lat = lat - degrees_per_px * y as f64; + let tx = (xtile * 256.0 + x as f64) / (2.0_f64.powi(z as i32) * 256.0); + let ty = (ytile * 256.0 + y as f64) / (2.0_f64.powi(z as i32) * 256.0); + + /*let lon = lon + degrees_lon_per_px * x as f64; + let lat = lat - degrees_lat_per_px * y as f64; + + */ + + let lon = (2.0_f64 * PI * tx - PI).to_degrees(); + let lat = ((PI - 2.0_f64 * PI * ty).exp().atan() * 2.0_f64 - PI / 2.0_f64).to_degrees(); + let closest_lat_in_map = &closest_key(map, lat).unwrap(); if closest_lat_in_map.1 > 0.01 { - // too far + image.set(y, x, Color { red: 0, green: 0, blue: 0, alpha: 30 }); continue; } let row = map.get(&closest_lat_in_map.0).unwrap(); let closest_long_in_row = closest_key(row, lon).unwrap(); + let value = row.get(&closest_long_in_row.0).unwrap(); + if closest_long_in_row.1 > 0.01 { - // too far + image.set(y, x, Color { red: 0, green: 0, blue: 0, alpha: 30 }); continue; } - let value = row.get(&closest_long_in_row.0).unwrap(); - let color = match value { - -999.0 | -99.0 => Color { red: 0, green: 0, blue: 0, alpha: 5 }, + -999.0 => Color { red: 0, green: 0, blue: 0, alpha: 30 }, + -99.0 => Color { red: 0, green: 0, blue: 0, alpha: 0 }, value_at_pos => { match *value_at_pos { ..=-30.0 => color!(rgb: 165, 165, 165), diff --git a/wxbox-tiler/src/tile.rs b/wxbox-tiler/src/tile.rs new file mode 100644 index 0000000..db8d61d --- /dev/null +++ b/wxbox-tiler/src/tile.rs @@ -0,0 +1,184 @@ +//! Tilemath adapted from maplibre/martin + +use std::f64::consts::PI; + +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; + +/// Convert longitude and latitude to a tile (x,y) coordinates for a given zoom +#[must_use] +#[allow(clippy::cast_possible_truncation)] +#[allow(clippy::cast_sign_loss)] +pub fn tile_index(lng: f64, lat: f64, zoom: u8) -> (u32, u32) { + 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] { + assert!(zoom <= MAX_ZOOM, "zoom {zoom} must be <= {MAX_ZOOM}"); + + let tile_length = EARTH_CIRCUMFERENCE / f64::from(1_u32 << zoom); + + let left_down_bbox = tile_bbox(min_x, max_y, tile_length); + let right_top_bbox = tile_bbox(max_x, min_y, 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] +pub fn bbox_to_xyz(left: f64, bottom: f64, right: f64, top: f64, zoom: u8) -> (u32, u32, u32, u32) { + let (min_col, min_row) = tile_index(left, top, zoom); + let (max_col, max_row) = tile_index(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 +#[must_use] +#[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)] +pub fn get_zoom_precision(zoom: u8) -> usize { + assert!(zoom < MAX_ZOOM, "zoom {zoom} must be <= {MAX_ZOOM}"); + 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) +} + +/// transform WGS84 to `WebMercator` +// from https://github.com/Esri/arcgis-osm-editor/blob/e4b9905c264aa22f8eeb657efd52b12cdebea69a/src/OSMWeb10_1/Utils/WebMercator.cs +#[must_use] +pub fn wgs84_to_webmercator(lon: f64, lat: f64) -> (f64, f64) { + let x = lon * PI / 180.0 * EARTH_RADIUS; + + let rad = lat * PI / 180.0; + let sin = rad.sin(); + let y = EARTH_RADIUS / 2.0 * ((1.0 + sin) / (1.0 - sin)).ln(); + + (x, y) +} + +#[cfg(test)] +mod tests { + #![allow(clippy::unreadable_literal)] + + use std::fs::read; + + use approx::assert_relative_eq; + use insta::assert_snapshot; + + use super::*; + + #[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[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); + + 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], 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.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.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] + 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); + } +} \ No newline at end of file diff --git a/wxbox-web/src/routes/+page.svelte b/wxbox-web/src/routes/+page.svelte index 1ed3882..89814b6 100644 --- a/wxbox-web/src/routes/+page.svelte +++ b/wxbox-web/src/routes/+page.svelte @@ -15,7 +15,7 @@ } ).addTo(m); L.tileLayer( - 'http://100.64.0.1:8080/mrms_cref/{z}/{x}/{y}.png', + 'http://localhost:8080/mrms_cref/{z}/{x}/{y}.png', { attribution: '© NOAA', maxZoom: 19