functioning tileserver !!

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",

wxbox-pal/src/lib.rs

@@ -49,7 +49,7 @@ impl Palette {
         }
     }
 }
-
+/*
 #[cfg(test)]
 mod tests {
     use std::collections::BTreeMap;
@@ -86,3 +86,5 @@ mod tests {
         assert_eq!(pal.color_for(15.0), color!(rgb: 0, 165, 255));
     }
 }
+
+ */

wxbox-tiler/Cargo.toml

@@ -18,3 +18,7 @@ png = "0.17"
 mime = "0.3.17"
 eccodes = { version = "0.12", features = ["message_ndarray"] }
 ndarray = "0.16"
+
+[dev-dependencies]
+approx = "0.5"
+insta = "1.40"

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;

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<AppState>) -> HttpResponse {
+async fn mrms_cref(path: web::Path<(u8, u32, u32)>, data: Data<AppState>) -> 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 [min_lng,min_lat,max_lng,max_lat] = xyz_to_bbox(z, x, y, x, y);
-    let lat = (PI - (y as f64 / 2.0_f64.powi(z)) * 2.0 * PI).sinh().atan() * (180.0 / PI);
 
+    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<AppState>) -> Ht
 
         let mut image: Pixmap = Pixmap::new();
 
-        let tile_size = 360.0 / 2.0_f64.powi(z);
+        let degrees_lat_per_px = (max_lat - min_lat) / 256.0;
-        let degrees_per_px = tile_size / 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 tx = (xtile * 256.0 + x as f64) / (2.0_f64.powi(z as i32) * 256.0);
-                let lat = lat - degrees_per_px * y as f64;
+                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),

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

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