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

formatting
Some checks are pending
Verify Latest Dependencies / Verify Latest Dependencies (push) Waiting to run
Details
build and test / wxbox - latest (push) Waiting to run
Details

Browse source
This commit is contained in:
core 2025-04-05 22:06:56 -04:00
parent 913ec262bd
commit 23571406ea
Signed by: core
GPG key ID: FDBF740DADDCEECF
10 changed files with 152 additions and 105 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

60
crates/ar2/src/lib.rs
View file

@ -1,8 +1,8 @@
use std::fmt::Debug;
use nexrad_data::volume::File;
use nexrad_decode::messages::digital_radar_data::{GenericDataBlock, RadialStatus};
use nexrad_decode::messages::MessageContents;
use nexrad_decode::messages::digital_radar_data::{GenericDataBlock, RadialStatus};
use nexrad_decode::result::Error;
use std::fmt::Debug;
pub mod sites;
@ -39,7 +39,7 @@ pub struct MomentData {
pub offset: f32,
pub values: Vec<u8>,
pub start_range: u16,
pub sample_interval: u16
pub sample_interval: u16,
}
impl MomentData {
@ -97,8 +97,7 @@ pub fn parse(input: Vec<u8>) -> nexrad_data::result::Result<Scan> {
if let MessageContents::DigitalRadarData(radar_data_message) = contents {
if vcp.is_none() {
if let Some(volume_block) = &radar_data_message.volume_data_block {
vcp =
Some(volume_block.volume_coverage_pattern_number);
vcp = Some(volume_block.volume_coverage_pattern_number);
}
}
radials.push(into_radial(*radar_data_message)?);
@ -112,9 +111,15 @@ pub fn parse(input: Vec<u8>) -> nexrad_data::result::Result<Scan> {
})
}
fn into_radial(message: nexrad_decode::messages::digital_radar_data::Message) -> nexrad_data::result::Result<Radial> {
fn into_radial(
message: nexrad_decode::messages::digital_radar_data::Message,
) -> nexrad_data::result::Result<Radial> {
Ok(Radial {
collection_timestamp: message.header.date_time().ok_or(Error::MessageMissingDateError)?.timestamp_millis(),
collection_timestamp: message
.header
.date_time()
.ok_or(Error::MessageMissingDateError)?
.timestamp_millis(),
azimuth_number: message.header.azimuth_number,
azimuth_angle_degrees: message.header.azimuth_angle,
azimuth_spacing_degrees: message.header.azimuth_resolution_spacing as f32 * 0.5,
@ -123,11 +128,21 @@ fn into_radial(message: nexrad_decode::messages::digital_radar_data::Message) ->
elevation_number_degrees: message.header.elevation_angle,
reflectivity: message.reflectivity_data_block.map(|u| into_moment_data(u)),
velocity: message.velocity_data_block.map(|u| into_moment_data(u)),
spectrum_width: message.spectrum_width_data_block.map(|u| into_moment_data(u)),
differential_reflectivity: message.differential_reflectivity_data_block.map(|u| into_moment_data(u)),
differential_phase: message.differential_phase_data_block.map(|u| into_moment_data(u)),
correlation_coefficient: message.correlation_coefficient_data_block.map(|u| into_moment_data(u)),
specific_differential_phase: message.specific_diff_phase_data_block.map(|u| into_moment_data(u)),
spectrum_width: message
.spectrum_width_data_block
.map(|u| into_moment_data(u)),
differential_reflectivity: message
.differential_reflectivity_data_block
.map(|u| into_moment_data(u)),
differential_phase: message
.differential_phase_data_block
.map(|u| into_moment_data(u)),
correlation_coefficient: message
.correlation_coefficient_data_block
.map(|u| into_moment_data(u)),
specific_differential_phase: message
.specific_diff_phase_data_block
.map(|u| into_moment_data(u)),
})
}
@ -137,58 +152,41 @@ fn into_moment_data(block: GenericDataBlock) -> MomentData {
offset: block.header.offset,
values: block.encoded_data,
start_range: block.header.data_moment_range,
sample_interval: block.header.data_moment_range_sample_interval
sample_interval: block.header.data_moment_range_sample_interval,
}
}
impl Sweep {
pub fn new(elevation_number: u8, radials: Vec<Radial>) -> Self {
Self {
elevation_number,
radials,
}
}
pub fn from_radials(radials: Vec<Radial>) -> Vec<Self> {
let mut sweeps = Vec::new();
let mut sweep_elevation_number = None;
let mut sweep_radials = Vec::new();
for radial in radials {
if let Some(elevation_number) = sweep_elevation_number {
if elevation_number != radial.elevation_number {
sweeps.push(Sweep::new(elevation_number, sweep_radials));
sweep_radials = Vec::new();
}
}
sweep_elevation_number = Some(radial.elevation_number);
sweep_radials.push(radial);
}
sweeps
}
}
pub const DATA_BYTES: &[u8] = include_bytes!("../KCRP20170825_235733_V06");
pub const DATA_BYTES: &[u8] = include_bytes!("../KCRP20170825_235733_V06");

15
crates/ar2/src/main.rs
View file

@ -4,5 +4,16 @@ use wxbox_ar2::parse;
fn main() {
let f = fs::read("KCRP20170825_235733_V06").unwrap();
let f = parse(f).unwrap();
println!("{:?}", f.sweeps.get(0).unwrap().radials.get(0).unwrap().reflectivity.as_ref().unwrap());
}
println!(
"{:?}",
f.sweeps
.get(0)
.unwrap()
.radials
.get(0)
.unwrap()
.reflectivity
.as_ref()
.unwrap()
);
}

12
crates/ar2/src/sites/wsr88d.rs
View file

@ -83,8 +83,10 @@ mod tests {
let radar = SITES.sites.get("KCRP").unwrap();
let radar_theta = radar.long;
let radar_phi = radar.lat;
let radar_r = ((A.powi(2) * radar_phi.cos()).powi(2) + (B.powi(2) * radar_phi.sin()).powi(2)
/ (A * radar_phi.cos()).powi(2) + (B * radar_phi.sin()).powi(2)).sqrt();
let radar_r = ((A.powi(2) * radar_phi.cos()).powi(2)
+ (B.powi(2) * radar_phi.sin()).powi(2) / (A * radar_phi.cos()).powi(2)
+ (B * radar_phi.sin()).powi(2))
.sqrt();
let radar_x = radar_r * radar_theta.cos() * radar_phi.sin();
let radar_y = radar_r * radar_theta.sin() * radar_phi.sin();
@ -94,7 +96,6 @@ mod tests {
let measurement_phi = radar_phi;
let measurement_r = radar_r;
let measurement_x = measurement_r * measurement_theta.cos() * measurement_phi.sin();
let measurement_y = measurement_r * measurement_theta.sin() * measurement_phi.sin();
let measurement_z = measurement_r * measurement_theta.cos();
@ -103,7 +104,8 @@ mod tests {
let radar_local_y = measurement_y - radar_y;
let radar_local_z = measurement_z - radar_z;
let radar_local_r = (radar_local_x.powi(2) + radar_local_y.powi(2) + radar_local_z.powi(2)).sqrt();
let radar_local_r =
(radar_local_x.powi(2) + radar_local_y.powi(2) + radar_local_z.powi(2)).sqrt();
let radar_local_theta = (radar_local_y / radar_local_x).atan();
let radar_local_phi = (radar_local_z / radar_local_r).acos();
@ -111,4 +113,4 @@ mod tests {
let elevation = radar_local_phi;
let distance = radar_local_r;
}
}
}

31
crates/client/src/map/sources/datalayers.rs
View file

@ -93,19 +93,18 @@ pub static NOAA_MRMS_MERGED_RHOHV_3KM_CONUS: LazyLock<LayerSource> =
tile_size: 512,
max_zoom: 19,
});
pub static ROC_KCRP_TEST: LazyLock<LayerSource> =
LazyLock::new(|| LayerSource {
source_id: 0x0EFAF68CB6E30B8E,
tile_url: format!(
"{}/nexrad/kcrp_ref_test/{{z}}/{{x}}/{{y}}@2x.png",
env!("TILER_BASE_URL")
),
display_name: "KCRP SSR Test".into(),
short_name: "KCRP SSR TEST (KCRP)".into(),
type_hint: LayerTypeHint::RadarData,
location: "KCRP".into(),
source: "NOAA / ROC".into(),
source_link: "https://roc.noaa.gov".into(),
tile_size: 512,
max_zoom: 19,
});
pub static ROC_KCRP_TEST: LazyLock<LayerSource> = LazyLock::new(|| LayerSource {
source_id: 0x0EFAF68CB6E30B8E,
tile_url: format!(
"{}/nexrad/kcrp_ref_test/{{z}}/{{x}}/{{y}}@2x.png",
env!("TILER_BASE_URL")
),
display_name: "KCRP SSR Test".into(),
short_name: "KCRP SSR TEST (KCRP)".into(),
type_hint: LayerTypeHint::RadarData,
location: "KCRP".into(),
source: "NOAA / ROC".into(),
source_link: "https://roc.noaa.gov".into(),
tile_size: 512,
max_zoom: 19,
});

2
crates/client/src/map/tiles.rs
View file

@ -62,7 +62,7 @@ impl LayerManager {
&*NOAA_MRMS_MERGED_CREF_QC_GUAM, // Guam
// -- NOAA MRMS - Merged RhoHV (Differential Reflectivity) @ 3km
&*NOAA_MRMS_MERGED_RHOHV_3KM_CONUS,
&*ROC_KCRP_TEST
&*ROC_KCRP_TEST,
];
for layer in layers {

2
crates/grib2/src/lib.rs
View file

@ -2,13 +2,13 @@ pub mod error;
pub mod wgs84;
use crate::error::GribError;
use wxbox_nommer::NomReader;
use crate::wgs84::LatLong;
use crate::LatLongVectorRelativity::{EasterlyAndNortherly, IncreasingXY};
use image::{DynamicImage, ImageFormat, ImageReader};
use std::fmt::{Debug, Formatter};
use std::io::{Cursor, Read};
use tracing::warn;
use wxbox_nommer::NomReader;
pub const INDICATOR: u32 = u32::from_be_bytes(*b"GRIB");
pub const EDITION: u8 = 2;

2
crates/tiler/src/config.rs
View file

@ -1,6 +1,6 @@
use crate::grib2::Grib2DataConfig;
use serde::{Deserialize, Serialize};
use crate::nexrad::NexradDataConfig;
use serde::{Deserialize, Serialize};
#[derive(Serialize, Deserialize)]
pub struct Config {

32
crates/tiler/src/grib2.rs
View file

@ -9,6 +9,7 @@ use flate2::read::GzDecoder;
use image::codecs::png::PngEncoder;
use image::{Rgba, RgbaImage};
use moka::future::Cache;
use rayon::iter::IntoParallelIterator;
use serde::{Deserialize, Serialize};
use std::collections::HashMap;
use std::f64::consts::PI;
@ -17,7 +18,6 @@ use std::io;
use std::io::{Cursor, ErrorKind};
use std::num::TryFromIntError;
use std::sync::Arc;
use rayon::iter::IntoParallelIterator;
use tokio::io::AsyncReadExt;
use tracing::{debug, info_span};
use wxbox_grib2::GribMessage;
@ -146,21 +146,29 @@ async fn render_to_png(
let tile_y_times_tilesize = tile_id.y as f64 * tile_id.size as f64;
let generate_pixels_span = info_span!("generate_pixels");
let pixels = (0..tile_id.size).into_par_iter().map(|x| {
(0..tile_id.size).into_par_iter().map(|y| {
let x_cartesian = (tile_x_times_tilesize + x as f64) / n;
let y_cartesian = (tile_y_times_tilesize + y as f64) / n;
let pixels = (0..tile_id.size)
.into_par_iter()
.map(|x| {
(0..tile_id.size)
.into_par_iter()
.map(|y| {
let x_cartesian = (tile_x_times_tilesize + x as f64) / n;
let y_cartesian = (tile_y_times_tilesize + y as f64) / n;
let long = (TWO_PI * x_cartesian - PI).to_degrees();
let lat = ((PI - TWO_PI * y_cartesian).exp().atan() * 2.0_f64 - HALF_PI).to_degrees();
let long = (TWO_PI * x_cartesian - PI).to_degrees();
let lat =
((PI - TWO_PI * y_cartesian).exp().atan() * 2.0_f64 - HALF_PI).to_degrees();
let nearest_value = data.value_for(LatLong { lat, long }).map(|u| u as f64);
let nearest_value = data.value_for(LatLong { lat, long }).map(|u| u as f64);
let color = colorize(nearest_value, &data_source).unwrap_or(Rgba::from([0,0,0,0]));
let color =
colorize(nearest_value, &data_source).unwrap_or(Rgba::from([0, 0, 0, 0]));
color
}).collect::<Vec<_>>()
}).collect::<Vec<_>>();
color
})
.collect::<Vec<_>>()
})
.collect::<Vec<_>>();
drop(generate_pixels_span);
let put_pixels_span = info_span!("put_pixels");

4
crates/tiler/src/main.rs
View file

@ -1,11 +1,12 @@
mod config;
mod error;
mod grib2;
mod tiles;
mod nexrad;
mod tiles;
use crate::config::Config;
use crate::grib2::{Grib2DataCache, Grib2TileCache, grib2_handler};
use crate::nexrad::{NexradDataCache, NexradTileCache, nexrad_handler};
use crate::tiles::{DataId, TileId};
use axum::Router;
use axum::routing::get;
@ -16,7 +17,6 @@ use std::sync::Arc;
use tracing_subscriber::fmt::format::FmtSpan;
use tracing_subscriber::util::SubscriberInitExt;
use wxbox_grib2::GribMessage;
use crate::nexrad::{nexrad_handler, NexradDataCache, NexradTileCache};
#[derive(Clone)]
pub struct AppState {

97
crates/tiler/src/nexrad.rs
View file

@ -6,9 +6,12 @@ use axum::extract::{Path, State};
use axum::http::{StatusCode, header};
use axum::response::IntoResponse;
use flate2::read::GzDecoder;
use geographiclib::Geodesic;
use image::codecs::png::PngEncoder;
use image::{Rgba, RgbaImage};
use moka::future::Cache;
use rayon::iter::IntoParallelIterator;
use rayon::iter::ParallelIterator;
use serde::{Deserialize, Serialize};
use std::collections::HashMap;
use std::f64::consts::PI;
@ -17,14 +20,11 @@ use std::io;
use std::io::{Cursor, ErrorKind};
use std::num::TryFromIntError;
use std::sync::Arc;
use geographiclib::Geodesic;
use rayon::iter::IntoParallelIterator;
use tokio::io::AsyncReadExt;
use tracing::{debug, info_span};
use wxbox_ar2::{parse, MomentValue, Scan, DATA_BYTES, Radial, Sweep};
use wxbox_ar2::sites::wsr88d::{SITES, Wsr88dSite};
use wxbox_ar2::{DATA_BYTES, MomentValue, Radial, Scan, Sweep, parse};
use wxbox_grib2::GribMessage;
use rayon::iter::ParallelIterator;
use wxbox_grib2::wgs84::LatLong;
use wxbox_pal::{ColorPalette, Palette};
@ -126,7 +126,15 @@ fn radius_at_latitude(phi: f64) -> f64 {
(top / bottom).sqrt()
}
fn calculate_inverse(tile_x_times_tilesize: f64, tile_y_times_tilesize: f64, x: usize, y: usize, g: &Geodesic, radar: &Wsr88dSite, n: f64) -> (f64, f64) {
fn calculate_inverse(
tile_x_times_tilesize: f64,
tile_y_times_tilesize: f64,
x: usize,
y: usize,
g: &Geodesic,
radar: &Wsr88dSite,
n: f64,
) -> (f64, f64) {
let x_cartesian = (tile_x_times_tilesize + x as f64) / n;
let y_cartesian = (tile_y_times_tilesize + y as f64) / n;
@ -145,7 +153,9 @@ fn locate_radial(azimuth: f64, sweep: &Sweep) -> Option<&Radial> {
for radial in &sweep.radials {
let this_dist = (radial.azimuth_angle_degrees as f64 - azimuth).abs();
if this_dist.abs() > radial.azimuth_spacing_degrees as f64 { continue; }
if this_dist.abs() > radial.azimuth_spacing_degrees as f64 {
continue;
}
if let Some(previous_best) = closest_radial {
let curr_dist = (previous_best.azimuth_angle_degrees as f64 - azimuth).abs();
@ -156,8 +166,6 @@ fn locate_radial(azimuth: f64, sweep: &Sweep) -> Option<&Radial> {
} else {
closest_radial = Some(radial);
}
}
closest_radial
@ -189,39 +197,58 @@ async fn render_to_png(
drop(prep_span);
let generate_pixels_span = info_span!("generate_pixels");
let pixels = (0..tile_id.size).into_par_iter().map(|x| {
(0..tile_id.size).into_par_iter().map(|y| {
let (d_m, azimuth) = calculate_inverse(tile_x_times_tilesize, tile_y_times_tilesize, x, y, &g, radar, n);
let Some(radial) = locate_radial(azimuth, first_sweep) else { return colorize(None, &data_source).unwrap_or(Rgba::from([0,0,0,0])) };
let pixels = (0..tile_id.size)
.into_par_iter()
.map(|x| {
(0..tile_id.size)
.into_par_iter()
.map(|y| {
let (d_m, azimuth) = calculate_inverse(
tile_x_times_tilesize,
tile_y_times_tilesize,
x,
y,
&g,
radar,
n,
);
let Some(radial) = locate_radial(azimuth, first_sweep) else {
return colorize(None, &data_source).unwrap_or(Rgba::from([0, 0, 0, 0]));
};
let reflectivity = radial.reflectivity.as_ref().unwrap();
let reflectivity = radial.reflectivity.as_ref().unwrap();
let distance_km = d_m / 1000.0;
let distance_km = d_m / 1000.0;
let first_gate_distance = reflectivity.start_range as f64 / 1000.0; // stored in meters
let gate_spacing = reflectivity.sample_interval as f64 / 1000.0; // also stored in meters
let first_gate_distance = reflectivity.start_range as f64 / 1000.0; // stored in meters
let gate_spacing = reflectivity.sample_interval as f64 / 1000.0; // also stored in meters
if distance_km < first_gate_distance {
return colorize(None, &data_source).unwrap_or(Rgba::from([0,0,0,0]))
}
if distance_km < first_gate_distance {
return colorize(None, &data_source).unwrap_or(Rgba::from([0, 0, 0, 0]));
}
let gate = ((distance_km - first_gate_distance) / gate_spacing).round() as usize;
let gate =
((distance_km - first_gate_distance) / gate_spacing).round() as usize;
let data = &reflectivity.values;
let data = &reflectivity.values;
if gate > data.len() || gate < 0 {
return colorize(None, &data_source).unwrap_or(Rgba::from([0,0,0,0]));
}
if gate > data.len() || gate < 0 {
return colorize(None, &data_source).unwrap_or(Rgba::from([0, 0, 0, 0]));
}
let color = data.get(gate).map(|raw_value| match raw_value {
0 => MomentValue::BelowThreshold,
1 => MomentValue::RangeFolded,
_ => MomentValue::Value((*raw_value as f32 - reflectivity.offset) / reflectivity.scale),
});
let color = data.get(gate).map(|raw_value| match raw_value {
0 => MomentValue::BelowThreshold,
1 => MomentValue::RangeFolded,
_ => MomentValue::Value(
(*raw_value as f32 - reflectivity.offset) / reflectivity.scale,
),
});
return colorize(color, &data_source).unwrap_or(Rgba::from([0,0,0,0]));
}).collect::<Vec<_>>()
}).collect::<Vec<_>>();
return colorize(color, &data_source).unwrap_or(Rgba::from([0, 0, 0, 0]));
})
.collect::<Vec<_>>()
})
.collect::<Vec<_>>();
drop(generate_pixels_span);
let put_pixels_span = info_span!("put_pixels");
@ -242,7 +269,6 @@ async fn render_to_png(
let output = Arc::new(result);
drop(encode_span);
drop(span);
@ -251,7 +277,10 @@ async fn render_to_png(
Ok(output)
}
fn colorize(value: Option<MomentValue>, data_source: &NexradDataSource) -> anyhow::Result<Rgba<u8>> {
fn colorize(
value: Option<MomentValue>,
data_source: &NexradDataSource,
) -> anyhow::Result<Rgba<u8>> {
Ok(match value {
Some(MomentValue::BelowThreshold) => Rgba([0, 0, 0, 0]),
Some(MomentValue::RangeFolded) => Rgba([119, 0, 125, 255]),