abstracted grib2 tile handler

This commit is contained in:
core 2024-10-24 22:02:23 -04:00
parent d292ae0193
commit 814daf48cf
Signed by: core
GPG Key ID: FDBF740DADDCEECF
5 changed files with 188 additions and 185 deletions

View File

@ -17,9 +17,9 @@ use crate::sources::noaa::mrms_cref;
static GLOBAL: tikv_jemallocator::Jemalloc = tikv_jemallocator::Jemalloc; static GLOBAL: tikv_jemallocator::Jemalloc = tikv_jemallocator::Jemalloc;
#[repr(usize)] #[repr(usize)]
#[derive(Ord, PartialOrd, Eq, PartialEq, Debug, Clone)] #[derive(Ord, PartialOrd, Eq, PartialEq, Debug, Clone, Copy)]
pub enum LutKey { pub enum LutKey {
NoaaMrmsCref = 1 NoaaMrmsMergedCompositeReflectivityQc = 1
} }
pub struct AppState { pub struct AppState {

View File

@ -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<BTreeMap<LutKey, SystemTime>>, 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<u8> {
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<u8>) -> 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<BTreeMap<LutKey, SystemTime>>, lut_cache: &RwLock<BTreeMap<LutKey, LookupTable2D>>) {
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<u8> {
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<u8> = 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
}

View File

@ -1 +1,2 @@
pub mod noaa; pub mod noaa;
mod grib2;

View File

@ -10,196 +10,30 @@ use flate2::read::GzDecoder;
use ndarray::{Zip}; use ndarray::{Zip};
use ordered_float::{FloatCore, OrderedFloat}; use ordered_float::{FloatCore, OrderedFloat};
use png::{BitDepth, ColorType, Encoder}; 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::{AppState, LutKey};
use crate::grib2::{closest_key, LookupTable2D}; use crate::grib2::{closest_key, LookupTable2D};
use crate::pixmap::Pixmap; 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<AppState>) -> HttpResponse {
let z = path.0;
let x = path.1;
let y = path.2;
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;
}
#[get("/noaa_mrms_merged_composite_reflectivity_qc/{z}/{x}/{y}.png")]
async fn mrms_cref(path: web::Path<(i32, u32, u32)>, data: Data<AppState>) -> HttpResponse {
let pal = create_test_palette(); let pal = create_test_palette();
let color = match value { reload_if_required(
-999.0 => Color { red: 0, green: 0, blue: 0, alpha: 30 }, "https://mrms.ncep.noaa.gov/data/2D/MergedReflectivityQCComposite/MRMS_MergedReflectivityQCComposite.latest.grib2.gz",
-99.0 => Color { red: 0, green: 0, blue: 0, alpha: 0 }, true,
value_at_pos => { 120,
pal.colorize(*value_at_pos) LutKey::NoaaMrmsMergedCompositeReflectivityQc,
/* &data.lut_cache_timestamps,
match *value_at_pos { &data.lut_cache
..=-30.0 => color!(rgb: 165, 165, 165), ).await;
-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<u8> = 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();
}
if let Some(map) = data.lut_cache.read().await.get(&LutKey::NoaaMrmsMergedCompositeReflectivityQc) {
HttpResponse::Ok() HttpResponse::Ok()
.insert_header(header::ContentType(mime::IMAGE_PNG)) .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 { } else {
println!("gridded LUT not available for LutKey::NoaaMrmsCref while handling /mrms_cref/{z}/{x}/{y} tile request");
HttpResponse::new(StatusCode::NOT_FOUND) HttpResponse::new(StatusCode::NOT_FOUND)
} }
} }

View File

@ -19,7 +19,7 @@
"OpenStreetMap": L.tileLayer("https://tile.openstreetmap.org/{z}/{x}/{y}.png", { attribution: "&copy;" }) "OpenStreetMap": L.tileLayer("https://tile.openstreetmap.org/{z}/{x}/{y}.png", { attribution: "&copy;" })
}, },
{ {
"NOAA/MRMS CONUS Reflectivity at Lowest Altitude": L.tileLayer("http://localhost:8080/mrms_cref/{z}/{x}/{y}.png", { attribution: "&copy; NOAA" }) "NOAA/MRMS CONUS Reflectivity at Lowest Altitude": L.tileLayer("http://localhost:8080/noaa_mrms_merged_composite_reflectivity_qc/{z}/{x}/{y}.png", { attribution: "&copy; NOAA" })
} }
).addTo(map); ).addTo(map);