From 814daf48cf3ea32dd3b700ba86163295a0abe96e Mon Sep 17 00:00:00 2001 From: core Date: Thu, 24 Oct 2024 22:02:23 -0400 Subject: [PATCH] abstracted grib2 tile handler --- wxbox-tiler/src/main.rs | 4 +- wxbox-tiler/src/sources/grib2.rs | 168 ++++++++++++++++++++++++ wxbox-tiler/src/sources/mod.rs | 3 +- wxbox-tiler/src/sources/noaa/mod.rs | 196 +++------------------------- wxbox-web/src/routes/+page.svelte | 2 +- 5 files changed, 188 insertions(+), 185 deletions(-) create mode 100644 wxbox-tiler/src/sources/grib2.rs diff --git a/wxbox-tiler/src/main.rs b/wxbox-tiler/src/main.rs index e46c1e5..be74263 100644 --- a/wxbox-tiler/src/main.rs +++ b/wxbox-tiler/src/main.rs @@ -17,9 +17,9 @@ use crate::sources::noaa::mrms_cref; static GLOBAL: tikv_jemallocator::Jemalloc = tikv_jemallocator::Jemalloc; #[repr(usize)] -#[derive(Ord, PartialOrd, Eq, PartialEq, Debug, Clone)] +#[derive(Ord, PartialOrd, Eq, PartialEq, Debug, Clone, Copy)] pub enum LutKey { - NoaaMrmsCref = 1 + NoaaMrmsMergedCompositeReflectivityQc = 1 } pub struct AppState { diff --git a/wxbox-tiler/src/sources/grib2.rs b/wxbox-tiler/src/sources/grib2.rs new file mode 100644 index 0000000..55cfcc5 --- /dev/null +++ b/wxbox-tiler/src/sources/grib2.rs @@ -0,0 +1,168 @@ +use std::collections::BTreeMap; +use std::f64::consts::PI; +use std::io::{BufWriter, Cursor, Read}; +use std::time::SystemTime; +use eccodes::{CodesHandle, FallibleStreamingIterator, ProductKind}; +use flate2::read::GzDecoder; +use ndarray::Zip; +use ordered_float::OrderedFloat; +use png::{BitDepth, ColorType, Encoder}; +use tokio::sync::RwLock; +use wxbox_pal::{Color, ColorPalette, Palette}; +use crate::grib2::{closest_key, LookupTable2D}; +use crate::LutKey; +use crate::pixmap::Pixmap; + +pub async fn needs_reload(lct: &RwLock>, lutkey: &LutKey, valid_for: u64) -> bool { + let lct_reader = lct.read().await; + + if let Some(t) = lct_reader.get(&lutkey) { + let dur = SystemTime::now().duration_since(*t).expect("time went backwards").as_secs(); + if dur > valid_for { + return true; + } + } else { + return true; + } + + false +} + +pub async fn load(url: &str, is_gzipped: bool) -> Vec { + let f = reqwest::get(url).await.unwrap(); + let data = f.bytes().await.unwrap().to_vec(); + + let out; + + if is_gzipped { + let b = Cursor::new(data); + let mut decoder = GzDecoder::new(b); + + let mut data = vec![]; + decoder.read_to_end(&mut data).expect("decoding failure"); + + out = data; + } else { + out = data; + } + + out +} + +pub fn eccodes_remap(data: Vec) -> LookupTable2D { + let product_kind = ProductKind::GRIB; + + let mut handle = CodesHandle::new_from_memory(data.clone(), product_kind).unwrap(); + let msg = handle.next().expect("empty grib2 file").expect("empty grib2 file").try_clone().unwrap(); + + let mut map = LookupTable2D::new(); + + let arr = msg.to_lons_lats_values().unwrap(); + + Zip::from(&arr.latitudes) + .and(&arr.longitudes) + .and(&arr.values) + .for_each(|&lat, &long, &value| { + let lat = OrderedFloat(lat); + + let mut long = long; + if long > 180.0 { + long -= 360.0; + } + + let long = OrderedFloat(long); + + if let Some(row) = map.get_mut(&lat) { + row.insert(long, value); + } else { + map.insert(lat, BTreeMap::from([ + (long, value) + ])); + } + }); + + map +} + +pub async fn reload_if_required(from: &str, needs_gzip: bool, valid_for: u64, lut_key: LutKey, lct_cache: &RwLock>, lut_cache: &RwLock>) { + if needs_reload(&lct_cache, &lut_key, valid_for).await { + let mut lct_writer = lct_cache.write().await; + + let map = eccodes_remap(load(from, needs_gzip).await); + + lut_cache.write().await.insert(lut_key, map); + lct_writer.insert(lut_key, SystemTime::now()); + } +} + +const TWO_PI: f64 = PI * 2.0; +const HALF_PI: f64 = PI / 2.0; + +pub fn render(xtile: f64, ytile: f64, z: i32, tilesize: usize, pal: Palette, map: &LookupTable2D) -> Vec { + let mut image: Pixmap = Pixmap::new(); + + let denominator = 2.0_f64.powi(z) * tilesize as f64; + let xtile_times_tilesize = xtile * tilesize as f64; + let ytile_times_tilesize = ytile * tilesize as f64; + + for x in 0..tilesize { + for y in 0..tilesize { + let tx = (xtile_times_tilesize + x as f64) / denominator; + let ty = (ytile_times_tilesize + y as f64) / denominator; + + let lon = (TWO_PI * tx - PI).to_degrees(); + let lat = ((PI - TWO_PI * ty).exp().atan() * 2.0_f64 - HALF_PI).to_degrees(); + + + let closest_lat_in_map = &closest_key(map, lat).unwrap(); + + if closest_lat_in_map.1 > 0.01 { + 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 { + image.set(y, x, Color { red: 0, green: 0, blue: 0, alpha: 30 }); + continue; + } + + let color = match value { + -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 => { + pal.colorize(*value_at_pos) + } + }; + image.set(y, x, color); + } + } + + let mut buf: Vec = vec![]; + // borrow checker insanity + { + let mut cur: Cursor<_> = Cursor::new(&mut buf); + let w = &mut BufWriter::new(&mut cur); + let mut encoder = Encoder::new(w, tilesize as u32, tilesize as u32); + encoder.set_color(ColorType::Rgba); + encoder.set_depth(BitDepth::Eight); + encoder.set_source_gamma(png::ScaledFloat::from_scaled(45455)); + encoder.set_source_gamma(png::ScaledFloat::new(1.0 / 2.2)); + let source_chromaticities = png::SourceChromaticities::new( + (0.31270, 0.32900), + (0.64000, 0.33000), + (0.30000, 0.60000), + (0.15000, 0.06000) + ); + encoder.set_source_chromaticities(source_chromaticities); + let mut writer = encoder.write_header().unwrap(); + writer.write_image_data(&image.to_raw()).expect("failed to encode png"); + writer.finish().unwrap(); + } + + buf +} \ No newline at end of file diff --git a/wxbox-tiler/src/sources/mod.rs b/wxbox-tiler/src/sources/mod.rs index 3e48278..26dc6ae 100644 --- a/wxbox-tiler/src/sources/mod.rs +++ b/wxbox-tiler/src/sources/mod.rs @@ -1 +1,2 @@ -pub mod noaa; \ No newline at end of file +pub mod noaa; +mod grib2; \ No newline at end of file diff --git a/wxbox-tiler/src/sources/noaa/mod.rs b/wxbox-tiler/src/sources/noaa/mod.rs index 7cfc548..0067054 100644 --- a/wxbox-tiler/src/sources/noaa/mod.rs +++ b/wxbox-tiler/src/sources/noaa/mod.rs @@ -10,196 +10,30 @@ use flate2::read::GzDecoder; use ndarray::{Zip}; use ordered_float::{FloatCore, OrderedFloat}; use png::{BitDepth, ColorType, Encoder}; -use wxbox_pal::{Color, ColorPalette, create_test_palette, Palette}; +use wxbox_pal::{Color, ColorPalette, create_test_palette}; use crate::{AppState, LutKey}; use crate::grib2::{closest_key, LookupTable2D}; use crate::pixmap::Pixmap; -use crate::tile::xyz_to_bbox; +use crate::sources::grib2::{eccodes_remap, load, needs_reload, reload_if_required, render}; -#[get("/mrms_cref/{z}/{x}/{y}.png")] -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; +#[get("/noaa_mrms_merged_composite_reflectivity_qc/{z}/{x}/{y}.png")] +async fn mrms_cref(path: web::Path<(i32, u32, u32)>, data: Data) -> HttpResponse { + let pal = create_test_palette(); - 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; - - let lct_reader = data.lut_cache_timestamps.read().await; - - if let Some(t) = lct_reader.get(&LutKey::NoaaMrmsCref) { - let dur = SystemTime::now().duration_since(*t).expect("time went backwards").as_secs(); - if dur > 120 { - needs_reload = true; - } - } else { - needs_reload = true; - } - - drop(lct_reader); - - if needs_reload { - let mut lct_writer = data.lut_cache_timestamps.write().await; - let mut lc_writer = data.lut_cache.write().await; - let f = reqwest::get("https://mrms.ncep.noaa.gov/data/2D/ReflectivityAtLowestAltitude/MRMS_ReflectivityAtLowestAltitude.latest.grib2.gz").await.unwrap(); - let data = f.bytes().await.unwrap().to_vec(); - let b = Cursor::new(data); - let mut decoder = GzDecoder::new(b); - - let mut data = vec![]; - decoder.read_to_end(&mut data).expect("decoding failure"); - - let product_kind = ProductKind::GRIB; - - let mut handle = CodesHandle::new_from_memory(data.clone(), product_kind).unwrap(); - let msg = handle.next().expect("empty grib2 file").expect("empty grib2 file").try_clone().unwrap(); - - let mut map = LookupTable2D::new(); - - let arr = msg.to_lons_lats_values().unwrap(); - - Zip::from(&arr.latitudes) - .and(&arr.longitudes) - .and(&arr.values) - .for_each(|&lat, &long, &value| { - let lat = OrderedFloat(lat); - - let mut long = long; - if long > 180.0 { - long -= 360.0; - } - - let long = OrderedFloat(long); - - if let Some(row) = map.get_mut(&lat) { - row.insert(long, value); - } else { - map.insert(lat, BTreeMap::from([ - (long, value) - ])); - } - }); - - lc_writer.insert(LutKey::NoaaMrmsCref, map); - lct_writer.insert(LutKey::NoaaMrmsCref, SystemTime::now()); - } - - if let Some(map) = data.lut_cache.read().await.get(&LutKey::NoaaMrmsCref) { -/* - let pal = Palette { - colors: BTreeMap::from([ - color!(-30, 165, 165, 165, 0), - color!(10, 0, 165, 255), - color!(20, 16, 255, 8), - color!(35, 251, 238, 0), - color!(50, 255, 0, 0), - color!(65, 247, 1, 249), - color!(75, 255, 255, 255), - color!(85, 184, 184, 184), - color!(95, 184, 184, 184) - ]), - clear_value: -99.0 - }; - - - */ - let mut image: Pixmap = Pixmap::new(); - - 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 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 { - 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 { - image.set(y, x, Color { red: 0, green: 0, blue: 0, alpha: 30 }); - continue; - } - - let pal = create_test_palette(); - - let color = match value { - -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 => { - pal.colorize(*value_at_pos) - /* - match *value_at_pos { - ..=-30.0 => color!(rgb: 165, 165, 165), - -30.0..=10.0 => color!(rgb: 0, 165, 255), - 10.0..=20.0 => color!(rgb: 16, 255, 8), - 20.0..=35.0 => color!(rgb: 251, 238, 0), - 35.0..=50.0 => color!(rgb: 255, 0, 0), - 50.0..=65.0 => color!(rgb: 247, 1, 249), - 65.0..=75.0 => color!(rgb: 255, 255, 255), - 75.0.. => color!(rgb: 184, 184, 184), - _ => unreachable!() - } - - */ - } - }; - image.set(y, x, color); - } - } - - let mut buf: Vec = vec![]; - // borrow checker insanity - { - let mut cur: Cursor<_> = Cursor::new(&mut buf); - let w = &mut BufWriter::new(&mut cur); - let mut encoder = Encoder::new(w, 256, 256); - encoder.set_color(ColorType::Rgba); - encoder.set_depth(BitDepth::Eight); - encoder.set_source_gamma(png::ScaledFloat::from_scaled(45455)); - encoder.set_source_gamma(png::ScaledFloat::new(1.0 / 2.2)); - let source_chromaticities = png::SourceChromaticities::new( - (0.31270, 0.32900), - (0.64000, 0.33000), - (0.30000, 0.60000), - (0.15000, 0.06000) - ); - encoder.set_source_chromaticities(source_chromaticities); - let mut writer = encoder.write_header().unwrap(); - writer.write_image_data(&image.to_raw()).expect("failed to encode png"); - writer.finish().unwrap(); - } + reload_if_required( + "https://mrms.ncep.noaa.gov/data/2D/MergedReflectivityQCComposite/MRMS_MergedReflectivityQCComposite.latest.grib2.gz", + true, + 120, + LutKey::NoaaMrmsMergedCompositeReflectivityQc, + &data.lut_cache_timestamps, + &data.lut_cache + ).await; + if let Some(map) = data.lut_cache.read().await.get(&LutKey::NoaaMrmsMergedCompositeReflectivityQc) { HttpResponse::Ok() .insert_header(header::ContentType(mime::IMAGE_PNG)) - .body(buf) + .body(render(path.1 as f64, path.2 as f64, path.0, 256, pal, map)) } else { - println!("gridded LUT not available for LutKey::NoaaMrmsCref while handling /mrms_cref/{z}/{x}/{y} tile request"); HttpResponse::new(StatusCode::NOT_FOUND) } } diff --git a/wxbox-web/src/routes/+page.svelte b/wxbox-web/src/routes/+page.svelte index 9ab25dd..bcd4402 100644 --- a/wxbox-web/src/routes/+page.svelte +++ b/wxbox-web/src/routes/+page.svelte @@ -19,7 +19,7 @@ "OpenStreetMap": L.tileLayer("https://tile.openstreetmap.org/{z}/{x}/{y}.png", { attribution: "©" }) }, { - "NOAA/MRMS CONUS Reflectivity at Lowest Altitude": L.tileLayer("http://localhost:8080/mrms_cref/{z}/{x}/{y}.png", { attribution: "© NOAA" }) + "NOAA/MRMS CONUS Reflectivity at Lowest Altitude": L.tileLayer("http://localhost:8080/noaa_mrms_merged_composite_reflectivity_qc/{z}/{x}/{y}.png", { attribution: "© NOAA" }) } ).addTo(map);