core/wxbox
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

new tiler w/homerolled grib2

Browse source
This commit is contained in:
core 2025-01-01 19:53:23 -05:00
parent fa8f5d71df
commit 31065832b7
Signed by: core
GPG key ID: FDBF740DADDCEECF
16 changed files with 1393 additions and 335 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

799
Cargo.lock generated
View file

File diff suppressed because it is too large Load diff

7
Cargo.toml
View file

@ -1,7 +1,12 @@
[workspace]
resolver = "2"
members = ["wxbox-pal","wxbox-tiler"]
members = [ "wxbox-grib2","wxbox-pal","wxbox-tiler"]
[profile.release]
codegen-units = 1
lto = "fat"
[profile.dev.package.image]
opt-level = 3
[profile.dev.package.png]
opt-level = 3

BIN
data.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 1.8 MiB

10
wxbox-grib2/Cargo.toml Normal file
View file

@ -0,0 +1,10 @@
[package]
name = "wxbox-grib2"
version = "0.1.0"
edition = "2021"
[dependencies]
thiserror = "2"
tracing = "0.1"
png = "0.17"
image = "0.25"

42
wxbox-grib2/src/error.rs Normal file
View file

@ -0,0 +1,42 @@
use std::io;
use image::ImageError;
use png::DecodingError;
use thiserror::Error;
#[derive(Error, Debug)]
pub enum GribError {
#[error("i/o error")]
IoError(#[from] io::Error),
#[error("incorrect indicator: found {0}")]
IncorrectIndicator(u32),
#[error("incorrect edition: found {0}")]
IncorrectEdition(u8),
#[error("this library only supports global grids in table 3.1")]
NonTable31,
#[error("this library does not support reading octets 15-xx in section 3")]
ListOfNumbersNotSupported,
#[error("unsupported grid definition template {0}")]
UnsupportedGrid(u16),
#[error("unsupported product definition template {0}")]
UnsupportedProductDefTmpl(u16),
#[error("unsupported data representation template {0}")]
UnsupportedDataRepresentationTemplate(u16),
#[error("unsupported bitmap")]
UnsupportedBitmap(u8),
#[error("missing data representation")]
MissingDataRepresentation,
#[error("missing identification")]
MissingIdentification,
#[error("missing grid definition")]
MissingGridDefinition,
#[error("missing product definition")]
MissingProductDefinition,
#[error("missing bitmap")]
MissingBitmap,
#[error("missing data")]
MissingData,
#[error("image error")]
PNGImageError(#[from] DecodingError),
#[error("image error")]
ImageError(#[from] ImageError)
}

587
wxbox-grib2/src/lib.rs Normal file
View file

@ -0,0 +1,587 @@
pub mod error;
mod nommer;
pub mod wgs84;
use std::fmt::{Debug, Formatter};
use std::io::{Cursor, Read};
use std::time::Instant;
use image::codecs::png::PngDecoder;
use image::{DynamicImage, ImageBuffer, ImageDecoder, ImageFormat, ImageReader, Luma};
use tracing::{debug, warn};
use crate::error::GribError;
use crate::LatLongVectorRelativity::{EasterlyAndNortherly, IncreasingXY};
use crate::nommer::NomReader;
use crate::wgs84::LatLong;
pub const INDICATOR: u32 = u32::from_be_bytes(*b"GRIB");
pub const EDITION: u8 = 2;
#[derive(Debug)]
pub struct GribMessage {
pub indicator: u32,
pub reserved: u16,
pub discipline: u8,
pub edition: u8,
pub message_length: u64,
pub identification: IdentificationSection,
pub grid_definition: GridDefinitionSection,
pub product_definition: ProductDefinitionSection,
pub data_representation: DataRepresentationSection,
pub bitmap: BitmapSection,
pub data: DataSection,
}
impl GribMessage {
pub fn new<R: Read>(reader: R) -> Result<Self, GribError> {
let mut nommer = NomReader::new(reader);
let indicator = nommer.read_u32()?;
if indicator != INDICATOR {
return Err(GribError::IncorrectIndicator(indicator));
}
let reserved = nommer.read_u16()?;
let discipline = nommer.read_u8()?;
let edition = nommer.read_u8()?;
if edition != EDITION {
return Err(GribError::IncorrectEdition(edition));
}
let message_length = nommer.read_u64()?;
let mut identification: Option<IdentificationSection> = None;
let mut grid_definition: Option<GridDefinitionSection> = None;
let mut product_definition: Option<ProductDefinitionSection> = None;
let mut data_representation: Option<DataRepresentationSection> = None;
let mut bitmap: Option<BitmapSection> = None;
let mut data: Option<DataSection> = None;
loop {
let length = match nommer.read_u32() {
Ok(l) => l,
Err(e) => return match e.kind() {
std::io::ErrorKind::UnexpectedEof => break,
_ => return Err(e.into())
}
};
if length == u32::from_be_bytes(*b"7777") {
break;
}
let section_number = nommer.read_u8()?;
match section_number {
1 => {
identification = Some(IdentificationSection::parse(length, &mut nommer)?);
},
3 => {
grid_definition = Some(GridDefinitionSection::parse(length, &mut nommer)?);
},
4 => {
product_definition = Some(ProductDefinitionSection::parse(length, &mut nommer)?);
},
5 => {
data_representation = Some(DataRepresentationSection::parse(length, &mut nommer)?);
},
6 => {
bitmap = Some(BitmapSection::parse(length, &mut nommer)?);
},
7 => {
data = Some(DataSection::parse(length, &mut nommer)?)
},
_ => {
let _ = nommer.read_n((length - 5) as usize);
warn!("unimplemented section # {} len {}", section_number, length)
}
}
}
let identification = identification.ok_or(GribError::MissingIdentification)?;
let grid_definition = grid_definition.ok_or(GribError::MissingGridDefinition)?;
let product_definition = product_definition.ok_or(GribError::MissingProductDefinition)?;
let mut data_representation = data_representation.ok_or(GribError::MissingDataRepresentation)?;
let bitmap = bitmap.ok_or(GribError::MissingBitmap)?;
let data = data.ok_or(GribError::MissingData)?;
data_representation.load_data(data.data.clone())?;
Ok(Self {
indicator,
reserved,
discipline,
edition,
message_length,
identification,
grid_definition,
product_definition,
data_representation,
bitmap,
data
})
}
pub fn value_for(&self, coord: LatLong) -> Option<f32> {
match self.grid_definition.image_coordinates(coord) {
Some((i, j)) => self.data_representation.get_image_coordinate(i as u32, j as u32),
None => None
}
}
}
#[derive(Debug)]
pub struct IdentificationSection {
pub originating_center_id: u16,
pub originating_subcenter_id: u16,
pub master_tables_version: u8,
pub local_tables_version: u8,
pub reference_time_significance: u8,
pub year: u16,
pub month: u8,
pub day: u8,
pub hour: u8,
pub minute: u8,
pub second: u8,
pub production_status: u8,
pub processed_data_type: u8,
pub reserved: Vec<u8>
}
impl IdentificationSection {
fn parse<R: Read>(length: u32, nommer: &mut NomReader<R>) -> Result<Self, GribError> {
let originating_center_id = nommer.read_u16()?;
let originating_subcenter_id = nommer.read_u16()?;
let master_tables_version = nommer.read_u8()?;
let local_tables_version = nommer.read_u8()?;
let reference_time_significance = nommer.read_u8()?;
let year = nommer.read_u16()?;
let month = nommer.read_u8()?;
let day = nommer.read_u8()?;
let hour = nommer.read_u8()?;
let minute = nommer.read_u8()?;
let second = nommer.read_u8()?;
let production_status = nommer.read_u8()?;
let processed_data_type = nommer.read_u8()?;
let reserved = nommer.read_n((length - 21) as usize)?;
Ok(Self {
originating_center_id,
originating_subcenter_id,
master_tables_version,
local_tables_version,
reference_time_significance,
year,
month,
day,
hour,
minute,
second,
production_status,
processed_data_type,
reserved
})
}
}
#[derive(Debug)]
pub struct GridDefinitionSection {
pub source: u8,
pub number_of_data_points: u32,
pub length_of_optional_number_list: u8,
pub interpretation_of_number_list: u8,
pub grid_definition_template_number: u16,
pub grid_definition: GridDefinition
}
impl GridDefinitionSection {
fn parse<R: Read>(length: u32, nommer: &mut NomReader<R>) -> Result<Self, GribError> {
let source = nommer.read_u8()?;
if source != 0 {
return Err(GribError::NonTable31);
}
let number_of_data_points = nommer.read_u32()?;
let length_of_optional_number_list = nommer.read_u8()?;
if length_of_optional_number_list != 0 {
return Err(GribError::ListOfNumbersNotSupported);
}
let interpretation_of_number_list = nommer.read_u8()?;
let grid_definition_template_number = nommer.read_u16()?;
let grid_definition = match grid_definition_template_number {
0 => GridDefinition::LatitudeLongitude(LatLongGrid::parse(length, nommer)?),
_ => return Err(GribError::UnsupportedGrid(grid_definition_template_number)),
};
Ok(Self {
source,
number_of_data_points,
length_of_optional_number_list,
interpretation_of_number_list,
grid_definition_template_number,
grid_definition
})
}
fn image_coordinates(&self, at: LatLong) -> Option<(u64, u64)> {
self.grid_definition.image_coordinates(at)
}
}
#[derive(Debug)]
pub struct ProductDefinitionSection {
pub number_of_coordinate_values_after_template: u16,
pub product_definition_template_number: u16,
pub product_definition: ProductDefinition,
}
impl ProductDefinitionSection {
fn parse<R: Read>(length: u32, nommer: &mut NomReader<R>) -> Result<Self, GribError> {
let number_of_coordinate_values_after_template = nommer.read_u16()?;
let product_definition_template_number = nommer.read_u16()?;
if number_of_coordinate_values_after_template != 0 {
return Err(GribError::ListOfNumbersNotSupported);
}
let product_definition = match product_definition_template_number {
0 => ProductDefinition::HorizontalLayerAtPointInTime(HorizontalLayerInstantaneousProductDef::parse(length, nommer)?),
_ => return Err(GribError::UnsupportedProductDefTmpl(product_definition_template_number))
};
Ok(Self {
number_of_coordinate_values_after_template,
product_definition_template_number,
product_definition
})
}
}
#[derive(Debug)]
pub struct DataRepresentationSection {
pub number_of_data_points: u32,
pub data_representation_template_number: u16,
pub data_representation_template: DataRepresentationTemplate
}
impl DataRepresentationSection {
fn parse<R: Read>(length: u32, nommer: &mut NomReader<R>) -> Result<Self, GribError> {
let number_of_data_points = nommer.read_u32()?;
let data_representation_template_number = nommer.read_u16()?;
let data_representation_template = match data_representation_template_number {
41 => DataRepresentationTemplate::GridpointPNG(GridpointPNGDataRepresentation::parse(length, nommer)?),
_ => return Err(GribError::UnsupportedDataRepresentationTemplate(data_representation_template_number))
};
Ok(Self {
number_of_data_points,
data_representation_template_number,
data_representation_template
})
}
fn load_data(&mut self, data: Vec<u8>) -> Result<(), GribError> {
self.data_representation_template.load_data(data)
}
fn get_image_coordinate(&self, x: u32, y: u32) -> Option<f32> {
self.data_representation_template.get_image_coordinate(x, y)
}
}
#[derive(Debug)]
pub enum DataRepresentationTemplate {
GridpointPNG(GridpointPNGDataRepresentation)
}
impl DataRepresentationTemplate {
fn load_data(&mut self, data: Vec<u8>) -> Result<(), GribError> {
match self {
Self::GridpointPNG(png) => png.load_data(data)
}
}
fn get_image_coordinate(&self, x: u32, y: u32) -> Option<f32> {
match self {
Self::GridpointPNG(png) => png.get_image_coordinate(x, y)
}
}
}
#[derive(Debug)]
pub struct GridpointPNGDataRepresentation {
pub reference_value: f32,
pub binary_scale_factor: u16,
pub decimal_scale_factor: u16,
pub depth: u8,
pub type_of_values: u8,
pub image: Option<ImageBuffer<Luma<u16>, Vec<u16>>>
}
impl GridpointPNGDataRepresentation {
fn parse<R: Read>(_length: u32, nommer: &mut NomReader<R>) -> Result<Self, GribError> {
let reference_value = nommer.read_f32()?;
let binary_scale_factor = nommer.read_u16()?;
let decimal_scale_factor = nommer.read_u16()?;
let depth = nommer.read_u8()?;
let type_of_values = nommer.read_u8()?;
Ok(Self {
reference_value,
binary_scale_factor,
decimal_scale_factor,
depth,
type_of_values,
image: None
})
}
fn load_data(&mut self, data: Vec<u8>) -> Result<(), GribError> {
let mut image_reader = ImageReader::new(Cursor::new(data));
image_reader.set_format(ImageFormat::Png);
let image = image_reader.decode()?;
self.image = Some(image.to_luma16());
Ok(())
}
fn get_image_coordinate(&self, x: u32, y: u32) -> Option<f32> {
match &self.image {
Some(i) => {
if x >= i.width() || y >= i.height() {
None
} else {
let datapoint = i.get_pixel(x, y).0[0] as f32;
let diff = datapoint * 2.0_f32.powi(self.binary_scale_factor as i32);
let dig_factor = 10_f32.powi(-(self.decimal_scale_factor as i32));
let value = (self.reference_value + diff) * dig_factor;
Some(value)
}
},
None => None,
}
}
}
#[derive(Debug)]
pub struct BitmapSection {
pub indicator: u8,
pub bitmap: Vec<u8>
}
impl BitmapSection {
fn parse<R: Read>(length: u32, nommer: &mut NomReader<R>) -> Result<Self, GribError> {
let indicator = nommer.read_u8()?;
if indicator != 255 {
return Err(GribError::UnsupportedBitmap(indicator));
}
let bitmap = nommer.read_n((length - 6) as usize)?;
Ok(Self {
indicator,
bitmap
})
}
}
pub struct DataSection {
pub data: Vec<u8>
}
impl DataSection {
fn parse<R: Read>(length: u32, nommer: &mut NomReader<R>) -> Result<Self, GribError> {
let data = nommer.read_n((length - 5) as usize)?;
Ok(Self {
data
})
}
}
impl Debug for DataSection {
fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result {
write!(f, "[data - {} bytes]", self.data.len())
}
}
#[derive(Debug)]
pub enum ProductDefinition {
HorizontalLayerAtPointInTime(HorizontalLayerInstantaneousProductDef),
}
#[derive(Debug)]
pub struct HorizontalLayerInstantaneousProductDef {
pub parameter_category: u8,
pub parameter_number: u8,
pub type_of_generating_process: u8,
pub background_generating_process_id: u8,
pub type_of_generating_process_identified: u8,
pub obs_data_cutoff_hours: u16,
pub obs_data_cutoff_minutes: u8,
pub indicator_of_unit_in_time_range: u8,
pub forecast_time: u32,
pub first_fixed_surface_type: u8,
pub scale_factor_first_fixed_surface: u8,
pub scaled_value_first_fixed_surface: u32,
pub second_fixed_surface_type: u8,
pub scale_factor_second_fixed_surface: u8,
pub scaled_value_second_fixed_surface: u32
}
impl HorizontalLayerInstantaneousProductDef {
fn parse<R: Read>(_length: u32, nommer: &mut NomReader<R>) -> Result<Self, GribError> {
let parameter_category = nommer.read_u8()?;
let parameter_number = nommer.read_u8()?;
let type_of_generating_process = nommer.read_u8()?;
let background_generating_process_id = nommer.read_u8()?;
let type_of_generating_process_identified = nommer.read_u8()?;
let obs_data_cutoff_hours = nommer.read_u16()?;
let obs_data_cutoff_minutes = nommer.read_u8()?;
let indicator_of_unit_in_time_range = nommer.read_u8()?;
let forecast_time = nommer.read_u32()?;
let first_fixed_surface_type = nommer.read_u8()?;
let scale_factor_first_fixed_surface = nommer.read_u8()?;
let scaled_value_first_fixed_surface = nommer.read_u32()?;
let second_fixed_surface_type = nommer.read_u8()?;
let scale_factor_second_fixed_surface = nommer.read_u8()?;
let scaled_value_second_fixed_surface = nommer.read_u32()?;
Ok(Self {
parameter_category,
parameter_number,
type_of_generating_process,
background_generating_process_id,
type_of_generating_process_identified,
obs_data_cutoff_hours,
obs_data_cutoff_minutes,
indicator_of_unit_in_time_range,
forecast_time,
first_fixed_surface_type,
scale_factor_first_fixed_surface,
scaled_value_first_fixed_surface,
second_fixed_surface_type,
scale_factor_second_fixed_surface,
scaled_value_second_fixed_surface,
})
}
}
#[derive(Debug)]
pub enum GridDefinition {
LatitudeLongitude(LatLongGrid)
}
impl GridDefinition {
fn image_coordinates(&self, at: LatLong) -> Option<(u64, u64)> {
match self {
GridDefinition::LatitudeLongitude(grid) => grid.image_coordinates(at)
}
}
}
#[derive(Debug)]
pub struct LatLongGrid {
pub shape_of_the_earth: u8,
pub scale_factor_of_radius_of_spherical_earth: u8,
pub scale_value_of_radius_of_spherical_earth: u32,
pub scale_factor_of_major_axis_of_oblate_spheroid_earth: u8,
pub scaled_value_of_major_axis_of_oblate_spheroid_earth: u32,
pub scale_factor_of_minor_axis_of_oblate_spheroid_earth: u8,
pub scaled_value_of_minor_axis_of_oblate_spheroid_earth: u32,
pub ni: u32, // number of points along a parallel
pub nj: u32, // number of points along a meridian
pub basic_angle_of_the_initial_production_domain: u32,
pub subdivisions_of_basic_angle: u32,
pub la1: f64,
pub lo1: f64,
pub i_direction_increments_given: bool,
pub j_direction_increments_given: bool,
pub vector_relativity: LatLongVectorRelativity,
pub la2: f64,
pub lo2: f64,
pub di: f64,
pub dj: f64,
pub scanning_mode_flags: u8,
}
#[derive(Debug)]
pub enum LatLongVectorRelativity {
EasterlyAndNortherly,
IncreasingXY
}
impl LatLongGrid {
fn parse<R: Read>(length: u32, nommer: &mut NomReader<R>) -> Result<Self, GribError> {
if length != 72 {
return Err(GribError::ListOfNumbersNotSupported);
}
let shape_of_the_earth= nommer.read_u8()?;
let scale_factor_of_radius_of_spherical_earth= nommer.read_u8()?;
let scale_value_of_radius_of_spherical_earth= nommer.read_u32()?;
let scale_factor_of_major_axis_of_oblate_spheroid_earth= nommer.read_u8()?;
let scaled_value_of_major_axis_of_oblate_spheroid_earth= nommer.read_u32()?;
let scale_factor_of_minor_axis_of_oblate_spheroid_earth= nommer.read_u8()?;
let scaled_value_of_minor_axis_of_oblate_spheroid_earth= nommer.read_u32()?;
let ni= nommer.read_u32()?; // number of points along a parallel
let nj= nommer.read_u32()?; // number of points along a meridian
let basic_angle_of_the_initial_production_domain= nommer.read_u32()?;
let subdivisions_of_basic_angle= nommer.read_u32()?;
let la1= nommer.read_u32()? as f64 * 10.0_f64.powi(-6);
let lo1= convert_longitude(nommer.read_u32()? as f64 * 10.0_f64.powi(-6));
let resolution_and_component_flags= nommer.read_u8()?;
let i_direction_increments_given = resolution_and_component_flags & (1 << 5) != 0;
let j_direction_increments_given = resolution_and_component_flags & (1 << 4) != 0;
let vector_relativity = if resolution_and_component_flags & (1 << 3) != 0 { IncreasingXY } else { EasterlyAndNortherly };
let la2= nommer.read_u32()? as f64 * 10.0_f64.powi(-6);
let lo2= convert_longitude(nommer.read_u32()? as f64 * 10.0_f64.powi(-6));
let di= nommer.read_u32()? as f64 * 10.0_f64.powi(-6);
let dj= nommer.read_u32()? as f64 * 10.0_f64.powi(-6);
let scanning_mode_flags= nommer.read_u8()?;
Ok(Self {
shape_of_the_earth,
scale_factor_of_radius_of_spherical_earth,
scale_value_of_radius_of_spherical_earth,
scale_factor_of_major_axis_of_oblate_spheroid_earth,
scaled_value_of_major_axis_of_oblate_spheroid_earth,
scale_factor_of_minor_axis_of_oblate_spheroid_earth,
scaled_value_of_minor_axis_of_oblate_spheroid_earth,
ni,
nj,
basic_angle_of_the_initial_production_domain,
subdivisions_of_basic_angle,
la1,
lo1,
i_direction_increments_given,
j_direction_increments_given,
vector_relativity,
la2,
lo2,
di,
dj,
scanning_mode_flags,
})
}
fn image_coordinates(&self, at: LatLong) -> Option<(u64, u64)> {
if at.lat > self.la1 || at.lat < self.la2 { return None; }
if at.long < self.lo1 || at.long > self.lo2 { return None; }
let lat = ((at.lat - self.la1) / self.di ).round() * self.di + self.la1;
let long = ((at.long - self.lo1) / self.dj).round() * self.dj + self.lo1;
if (lat - at.lat).abs() > self.di || (long - at.long).abs() > self.dj {
return None;
}
let diff_lat = at.lat - self.la1;
let num_steps_lat = diff_lat / -self.di;
let diff_long = at.long - self.lo1;
let num_steps_long = diff_long / self.dj;
Some((num_steps_long as u64, num_steps_lat as u64))
}
}
fn convert_longitude(longitude: f64) -> f64 {
if longitude > 180.0 {
longitude - 360.0
} else {
longitude
}
}

44
wxbox-grib2/src/nommer.rs Normal file
View file

@ -0,0 +1,44 @@
use std::io::Read;
#[derive(Debug)]
pub struct NomReader<R: Read> {
inner: R
}
impl<R: Read> NomReader<R> {
pub fn new(inner: R) -> NomReader<R> {
Self {
inner
}
}
pub fn read_u8(&mut self) -> Result<u8, std::io::Error> {
let mut buf = [0u8; 1];
self.inner.read_exact(&mut buf)?;
Ok(u8::from_be_bytes(buf))
}
pub fn read_u16(&mut self) -> Result<u16, std::io::Error> {
let mut buf = [0u8; 2];
self.inner.read_exact(&mut buf)?;
Ok(u16::from_be_bytes(buf))
}
pub fn read_u32(&mut self) -> Result<u32, std::io::Error> {
let mut buf = [0u8; 4];
self.inner.read_exact(&mut buf)?;
Ok(u32::from_be_bytes(buf))
}
pub fn read_u64(&mut self) -> Result<u64, std::io::Error> {
let mut buf = [0u8; 8];
self.inner.read_exact(&mut buf)?;
Ok(u64::from_be_bytes(buf))
}
pub fn read_n(&mut self, n: usize) -> Result<Vec<u8>, std::io::Error> {
let mut buf = vec![0u8; n];
self.inner.read_exact(&mut buf)?;
Ok(buf)
}
pub fn read_f32(&mut self) -> Result<f32, std::io::Error> {
let mut buf = [0u8; 4];
self.inner.read_exact(&mut buf)?;
Ok(f32::from_be_bytes(buf))
}
}

5
wxbox-grib2/src/wgs84.rs Normal file
View file

@ -0,0 +1,5 @@
#[derive(Debug, Copy, Clone)]
pub struct LatLong {
pub lat: f64,
pub long: f64
}

5
wxbox-tiler/Cargo.toml
View file

@ -16,8 +16,9 @@ tokio = "1"
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"
wxbox-grib2 = { version = "0.1", path = "../wxbox-grib2" }
tracing = "0.1"
tracing-subscriber = "0.3"
[dev-dependencies]
approx = "0.5"

55
wxbox-tiler/src/grib2.rs
View file

@ -1,55 +0,0 @@
use std::collections::BTreeMap;
use ordered_float::OrderedFloat;
pub fn closest_key<V>(map: &BTreeMap<OrderedFloat<f64>, V>, val: f64) -> Option<(OrderedFloat<f64>, f64)> {
let mut r1 = map.range(OrderedFloat(val)..);
let r2 = map.range(..=OrderedFloat(val));
let o1 = r1.next();
let o2 = r2.last();
match (o1, o2) {
(None, None) => {
None
},
(Some(i), None) => {
Some((*i.0, (*i.0 - val).abs()))
},
(None, Some(i)) => {
Some((*i.0, (*i.0 - val).abs()))
},
(Some(i1), Some(i2)) => {
// abs(f - i)
let i1_dist = (i1.0 - val).abs();
let i2_dist = (i2.0 - val).abs();
Some(if i1_dist < i2_dist {
(*i1.0, i1_dist)
} else {
(*i2.0, i2_dist)
})
}
}
}
pub fn lookup<V>(map: &BTreeMap<OrderedFloat<f64>, V>, val: f64) -> Option<&V> {
let mut r1 = map.range(OrderedFloat(val)..);
let mut r2 = map.range(..OrderedFloat(val));
let o1 = r1.next();
let o2 = r2.next();
match (o1, o2) {
(None, None) => None,
(Some(i), None) => Some(i.1),
(None, Some(i)) => Some(i.1),
(Some(i1), Some(i2)) => {
// abs(f - i)
let i1_dist = (i1.0 - val).abs();
let i2_dist = (i2.0 - val).abs();
Some(if i1_dist < i2_dist {
i1.1
} else {
i2.1
})
}
}
}
pub type LookupTable2D = BTreeMap<OrderedFloat<f64>, BTreeMap<OrderedFloat<f64>, f64>>;

27
wxbox-tiler/src/main.rs
View file

@ -1,39 +1,28 @@
pub(crate) mod sources;
mod grib2;
mod pixmap;
use std::collections::{BTreeMap};
use std::collections::{BTreeMap, HashMap};
use std::fmt::Debug;
use std::sync::Arc;
use tokio::sync::RwLock;
use std::time::SystemTime;
use actix_web::{App, HttpServer};
use actix_web::web::Data;
use crate::grib2::{LookupTable2D};
#[global_allocator]
static GLOBAL: tikv_jemallocator::Jemalloc = tikv_jemallocator::Jemalloc;
#[repr(usize)]
#[derive(Ord, PartialOrd, Eq, PartialEq, Debug, Clone, Copy)]
pub enum LutKey {
NoaaMrmsMergedCompositeReflectivityQc = 1
}
use wxbox_grib2::GribMessage;
pub struct AppState {
lut_cache: RwLock<BTreeMap<LutKey, LookupTable2D>>,
lut_cache_timestamps: RwLock<BTreeMap<LutKey, SystemTime>>
grib2_cache: RwLock<HashMap<String, Arc<RwLock<GribMessage>>>>,
grib2_cache_timestamps: RwLock<HashMap<String, SystemTime>>
}
#[actix_web::main]
async fn main() -> std::io::Result<()> {
let f = reqwest::get("https://mrms.ncep.noaa.gov/data/2D/ReflectivityAtLowestAltitude/MRMS_ReflectivityAtLowestAltitude.latest.grib2.gz").await.unwrap();
tracing_subscriber::fmt::init();
let data = Data::new(AppState {
lut_cache: RwLock::new(BTreeMap::new()),
lut_cache_timestamps: RwLock::new(BTreeMap::new())
grib2_cache: RwLock::new(HashMap::new()),
grib2_cache_timestamps: RwLock::new(HashMap::new())
});
HttpServer::new(move || {
App::new()

127
wxbox-tiler/src/sources/grib2.rs
View file

@ -1,24 +1,20 @@
use std::collections::BTreeMap;
use std::collections::{BTreeMap, HashMap};
use std::f64::consts::PI;
use std::io::{BufWriter, Cursor, Read};
use std::time::{SystemTime, UNIX_EPOCH};
use eccodes::{CodesHandle, FallibleStreamingIterator, ProductKind};
use std::sync::Arc;
use std::time::SystemTime;
use flate2::read::GzDecoder;
use ndarray::Zip;
use ordered_float::OrderedFloat;
use png::{BitDepth, ColorType, Encoder};
use tokio::sync::RwLock;
use wxbox_grib2::GribMessage;
use wxbox_grib2::wgs84::LatLong;
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 {
eprintln!("debug: try lock needs_reload() lut_cache_timestamps start");
pub async fn needs_reload(lct: &RwLock<HashMap<String, SystemTime>>, lutkey: &String, valid_for: u64) -> bool {
let lct_reader = lct.read().await;
eprintln!("debug: try lock needs_reload() lut_cache_timestamps locked");
if let Some(t) = lct_reader.get(&lutkey) {
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;
@ -51,66 +47,31 @@ pub async fn load(url: &str, is_gzipped: bool) -> Vec<u8> {
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 {
eprintln!("debug: try lock reload_if_required() lut_cache_timestamps start WRITE");
pub async fn reload_if_required(from: &str, needs_gzip: bool, valid_for: u64, lut_key: &String, lct_cache: &RwLock<HashMap<String, SystemTime>>, lut_cache: &RwLock<HashMap<String, Arc<RwLock<GribMessage>>>>) {
if needs_reload(&lct_cache, lut_key, valid_for).await {
let mut lct_writer = lct_cache.write().await;
eprintln!("debug: try lock reload_if_required() lut_cache_timestamps locked WRITE");
let map = eccodes_remap(load(from, needs_gzip).await);
let message = load(from, needs_gzip).await;
let grib = GribMessage::new(Cursor::new(message)).unwrap();
eprintln!("debug: try lock reload_if_required() lut_cache start");
lut_cache.write().await.insert(lut_key, map);
eprintln!("debug: try lock reload_if_required() lut_cache locked");
lct_writer.insert(lut_key, SystemTime::now());
lut_cache.write().await.insert(lut_key.clone(), Arc::new(RwLock::new(grib)));
lct_writer.insert(lut_key.clone(), 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, missing: Option<f64>, range_folded: Option<f64>, no_coverage: Option<f64>) -> Vec<u8> {
pub async fn render(xtile: f64, ytile: f64, z: i32, tilesize: usize, pal: Palette, map: &Arc<RwLock<GribMessage>>, missing: Option<f64>, range_folded: Option<f64>, no_coverage: Option<f64>) -> 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;
{
let message = map.read().await;
for x in 0..tilesize {
for y in 0..tilesize {
let tx = (xtile_times_tilesize + x as f64) / denominator;
@ -119,35 +80,25 @@ pub fn render(xtile: f64, ytile: f64, z: i32, tilesize: usize, pal: Palette, map
let lon = (TWO_PI * tx - PI).to_degrees();
let lat = ((PI - TWO_PI * ty).exp().atan() * 2.0_f64 - HALF_PI).to_degrees();
let nearest = message.value_for(LatLong {
lat,
long: lon
}).map(|u| u as f64);
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 {
c if Some(*c) == no_coverage => Color { red: 0, green: 0, blue: 0, alpha: 30 },
c if Some(*c) == missing => Color { red: 0, green: 0, blue: 0, alpha: 0 },
c if Some(*c) == range_folded => Color { red: 141, green: 0, blue: 160, alpha: 0 },
value_at_pos => {
pal.colorize(*value_at_pos)
}
let color = match nearest {
Some(c) if Some(c) == no_coverage => Color { red: 0, green: 0, blue: 0, alpha: 30 },
Some(c) if Some(c) == missing => Color { red: 0, green: 0, blue: 0, alpha: 0 },
Some(c) if Some(c) == range_folded => Color { red: 141, green: 0, blue: 160, alpha: 0 },
Some(value_at_pos) => {
pal.colorize(value_at_pos)
},
None => Color { red: 0, green: 0, blue: 0, alpha: 30 }
};
image.set(y, x, color);
}
}
}
let mut buf: Vec<u8> = vec![];
// borrow checker insanity
@ -176,28 +127,26 @@ pub fn render(xtile: f64, ytile: f64, z: i32, tilesize: usize, pal: Palette, map
#[macro_export]
macro_rules! grib2_handler {
(mount $f:ident, at: $path:expr, from: $from:expr, needs_gzip: $needs_gzip:expr, valid_for: $valid_for:expr, lut_with: $lut_key:expr, palette: $pal:expr, missing: $missing:expr, range_folded: $rf:expr, no_coverage: $nc:expr) => {
(mount $f:ident, at: $path:expr, from: $from:expr, needs_gzip: $needs_gzip:expr, valid_for: $valid_for:expr, palette: $pal:expr, missing: $missing:expr, range_folded: $rf:expr, no_coverage: $nc:expr) => {
#[::actix_web::get($path)]
async fn $f(path: ::actix_web::web::Path<(i32,u32,u32)>, data: ::actix_web::web::Data<crate::AppState>) -> ::actix_web::HttpResponse {
crate::sources::grib2::reload_if_required(
$from,
$needs_gzip,
$valid_for,
$lut_key,
&data.lut_cache_timestamps,
&data.lut_cache
&String::from($path),
&data.grib2_cache_timestamps,
&data.grib2_cache
).await;
eprintln!("debug: try lock handler lut_cache_timestamps start");
let lct_reader = data.lut_cache_timestamps.read().await;
eprintln!("debug: try lock handler lut_cache_timestamps locked");
if let Some(map) = data.lut_cache.read().await.get(&$lut_key) {
let lct_reader = data.grib2_cache_timestamps.read().await;
if let Some(map) = data.grib2_cache.read().await.get($path) {
::actix_web::HttpResponse::Ok()
.insert_header(::actix_web::http::header::ContentType(::mime::IMAGE_PNG))
.insert_header(("x-wxbox-tiler-data-valid-time", lct_reader.get(&$lut_key).expect("impossible").duration_since(::std::time::UNIX_EPOCH).expect("time went backwards").as_secs().to_string()))
.insert_header(("x-wxbox-tiler-data-valid-time", lct_reader.get($path).expect("impossible").duration_since(::std::time::UNIX_EPOCH).expect("time went backwards").as_secs().to_string()))
.insert_header(("Access-Control-Allow-Origin", "*"))
.insert_header(("Access-Control-Expose-Headers", "*"))
.insert_header(("Access-Control-Allow-Headers", "*"))
.body(crate::sources::grib2::render(path.1 as f64, path.2 as f64, path.0, 256, $pal, map, $missing, $rf, $nc))
.body(crate::sources::grib2::render(path.1 as f64, path.2 as f64, path.0, 256, $pal, map, $missing, $rf, $nc).await)
} else {
::actix_web::HttpResponse::new(::actix_web::http::StatusCode::NOT_FOUND)
}

3
wxbox-tiler/src/sources/noaa/mod.rs
View file

@ -1,4 +1,4 @@
use crate::{grib2_handler, LutKey};
use crate::grib2_handler;
grib2_handler! {
mount noaa_mrms_merged_composite_reflectivity_qc,
@ -6,7 +6,6 @@ grib2_handler! {
from: "https://mrms.ncep.noaa.gov/data/2D/MergedReflectivityQCComposite/MRMS_MergedReflectivityQCComposite.latest.grib2.gz",
needs_gzip: true,
valid_for: 120,
lut_with: LutKey::NoaaMrmsMergedCompositeReflectivityQc,
palette: wxbox_pal::parser::parse(wxbox_pal::default_palettes::DEFAULT_REFLECTIVITY_PALETTE).unwrap(),
missing: Some(-99.0),
range_folded: None,

1
wxbox-web/package.json
View file

@ -24,7 +24,6 @@
"prettier": "^3.3.3",
"prettier-plugin-svelte": "^3.2.7",
"svelte": "^5.1.3",
"svelte-adapter-bun": "^0.5.2",
"svelte-check": "^4.0.5",
"typescript": "^5.6.3",
"typescript-eslint": "^8.11.0",

2
wxbox-web/src/routes/+page.svelte
View file

@ -373,7 +373,7 @@
{/if}
<div class="footer text-sm">
<p>built with &lt;3</p>
<p>u8.lc coredoes.dev :)</p>
<p>u8.lc & coredoes.dev :)</p>
</div>
</div>

2
wxbox-web/svelte.config.js
View file

@ -1,4 +1,4 @@
import adapter from 'svelte-adapter-bun';
import adapter from '@sveltejs/adapter-auto';
import { vitePreprocess } from '@sveltejs/vite-plugin-svelte';
/** @type {import('@sveltejs/kit').Config} */