213 lines
6.5 KiB
Rust
213 lines
6.5 KiB
Rust
//! IEM n0q CONUS composite reflectivity fetch + per-POI extraction.
|
||
//!
|
||
//! Port of `Microwaveprop.Weather.NexradClient.fetch_frame/2` — the
|
||
//! subset used by f00's chain step. n0q is a palettized 8-bit PNG,
|
||
//! 12200 × 5400 pixels, 0.005°/pixel, covering (-126W..-65W, 23N..50N)
|
||
//! every 5 min. For each point of interest the extractor returns the
|
||
//! max dBZ in a ~25 km box, matching Elixir's
|
||
//! `extract_box/5` output field `max_reflectivity_dbz`.
|
||
|
||
use bytes::Bytes;
|
||
use chrono::{DateTime, Datelike, Timelike, Utc};
|
||
|
||
pub const LON_MIN: f64 = -126.0;
|
||
pub const LAT_MAX: f64 = 50.0;
|
||
pub const DEG_PER_PIXEL: f64 = 0.005;
|
||
pub const DEFAULT_BOX_HALF: i32 = 25;
|
||
pub const BASE_URL: &str = "https://mesonet.agron.iastate.edu/archive/data";
|
||
|
||
#[derive(Debug, thiserror::Error)]
|
||
pub enum NexradError {
|
||
#[error("http: {0}")]
|
||
Http(#[from] reqwest::Error),
|
||
#[error("decode: {0}")]
|
||
Decode(String),
|
||
#[error("unexpected status {0}")]
|
||
Status(u16),
|
||
}
|
||
|
||
#[derive(Debug, Clone, Copy)]
|
||
pub struct NexradObservation {
|
||
pub lat: f64,
|
||
pub lon: f64,
|
||
pub max_reflectivity_dbz: f64,
|
||
}
|
||
|
||
/// Round a timestamp down to the nearest 5-min boundary (IEM archives
|
||
/// 5-min cadence).
|
||
pub fn round_to_5min(ts: DateTime<Utc>) -> DateTime<Utc> {
|
||
let minute = ts.minute();
|
||
let rounded = (minute / 5) * 5;
|
||
ts.with_minute(rounded)
|
||
.and_then(|t| t.with_second(0))
|
||
.and_then(|t| t.with_nanosecond(0))
|
||
.expect("valid rounded timestamp")
|
||
}
|
||
|
||
pub fn frame_url(rounded: DateTime<Utc>) -> String {
|
||
format!(
|
||
"{}/{:04}/{:02}/{:02}/GIS/uscomp/n0q_{:04}{:02}{:02}{:02}{:02}.png",
|
||
BASE_URL,
|
||
rounded.year(),
|
||
rounded.month(),
|
||
rounded.day(),
|
||
rounded.year(),
|
||
rounded.month(),
|
||
rounded.day(),
|
||
rounded.hour(),
|
||
rounded.minute()
|
||
)
|
||
}
|
||
|
||
pub fn latlon_to_pixel(lat: f64, lon: f64) -> (i32, i32) {
|
||
let x = ((lon - LON_MIN) / DEG_PER_PIXEL).round() as i32;
|
||
let y = ((LAT_MAX - lat) / DEG_PER_PIXEL).round() as i32;
|
||
(x, y)
|
||
}
|
||
|
||
/// Value 0 = no echo. Values 1..=255 map linearly to -30..=+95 dBZ.
|
||
pub fn pixel_to_dbz(v: u8) -> f64 {
|
||
if v == 0 {
|
||
0.0
|
||
} else {
|
||
-30.0 + ((v - 1) as f64) / 254.0 * 125.0
|
||
}
|
||
}
|
||
|
||
/// Maximum dBZ within a 2·half+1 × 2·half+1 box around (cx, cy).
|
||
/// Pixels outside the image or zero (no echo) are skipped.
|
||
pub fn max_dbz_in_box(pixels: &[u8], width: i32, cx: i32, cy: i32, half: i32) -> f64 {
|
||
let height = (pixels.len() as i32) / width;
|
||
let x_min = (cx - half).max(0);
|
||
let x_max = (cx + half).min(width - 1);
|
||
let y_min = (cy - half).max(0);
|
||
let y_max = (cy + half).min(height - 1);
|
||
|
||
let mut max_val = 0u8;
|
||
for y in y_min..=y_max {
|
||
for x in x_min..=x_max {
|
||
let off = (y * width + x) as usize;
|
||
if off < pixels.len() {
|
||
let v = pixels[off];
|
||
if v > max_val {
|
||
max_val = v;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
pixel_to_dbz(max_val)
|
||
}
|
||
|
||
/// Decode an IEM n0q PNG to an indexed-pixel buffer.
|
||
pub fn decode_n0q_png(bytes: &[u8]) -> Result<(Vec<u8>, i32), NexradError> {
|
||
let decoder = png::Decoder::new(std::io::Cursor::new(bytes));
|
||
let mut reader = decoder
|
||
.read_info()
|
||
.map_err(|e| NexradError::Decode(format!("read_info: {e}")))?;
|
||
let info = reader.info();
|
||
let width = info.width as i32;
|
||
// IEM n0q is indexed-color 8-bit. `next_frame` fills the buffer
|
||
// with raw palette indices, which is exactly what pixel_to_dbz
|
||
// expects.
|
||
let mut buf = vec![0u8; reader.output_buffer_size()];
|
||
reader
|
||
.next_frame(&mut buf)
|
||
.map_err(|e| NexradError::Decode(format!("next_frame: {e}")))?;
|
||
Ok((buf, width))
|
||
}
|
||
|
||
/// Fetch + decode the n0q frame nearest `timestamp`, then return one
|
||
/// `max_reflectivity_dbz` per point of interest.
|
||
pub async fn fetch_frame(
|
||
client: &reqwest::Client,
|
||
timestamp: DateTime<Utc>,
|
||
points: &[(f64, f64)],
|
||
) -> Result<Vec<NexradObservation>, NexradError> {
|
||
let rounded = round_to_5min(timestamp);
|
||
let url = frame_url(rounded);
|
||
let resp = client
|
||
.get(&url)
|
||
.timeout(std::time::Duration::from_secs(60))
|
||
.send()
|
||
.await?;
|
||
if !resp.status().is_success() {
|
||
return Err(NexradError::Status(resp.status().as_u16()));
|
||
}
|
||
let body: Bytes = resp.bytes().await?;
|
||
let (pixels, width) = decode_n0q_png(&body)?;
|
||
|
||
let obs = points
|
||
.iter()
|
||
.map(|&(lat, lon)| {
|
||
let (cx, cy) = latlon_to_pixel(lat, lon);
|
||
let dbz = max_dbz_in_box(&pixels, width, cx, cy, DEFAULT_BOX_HALF);
|
||
NexradObservation {
|
||
lat,
|
||
lon,
|
||
max_reflectivity_dbz: dbz,
|
||
}
|
||
})
|
||
.collect();
|
||
Ok(obs)
|
||
}
|
||
|
||
#[cfg(test)]
|
||
mod tests {
|
||
use super::*;
|
||
use chrono::TimeZone;
|
||
|
||
#[test]
|
||
fn rounds_to_five_minute_boundary() {
|
||
let t = Utc.with_ymd_and_hms(2026, 4, 19, 14, 37, 45).unwrap();
|
||
let r = round_to_5min(t);
|
||
assert_eq!(r.minute(), 35);
|
||
assert_eq!(r.second(), 0);
|
||
assert_eq!(r.hour(), 14);
|
||
}
|
||
|
||
#[test]
|
||
fn frame_url_format() {
|
||
let t = Utc.with_ymd_and_hms(2026, 4, 19, 14, 35, 0).unwrap();
|
||
let u = frame_url(t);
|
||
assert_eq!(
|
||
u,
|
||
"https://mesonet.agron.iastate.edu/archive/data/2026/04/19/GIS/uscomp/n0q_202604191435.png"
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn pixel_to_dbz_endpoints() {
|
||
assert_eq!(pixel_to_dbz(0), 0.0);
|
||
assert!((pixel_to_dbz(1) - -30.0).abs() < 1e-9);
|
||
assert!((pixel_to_dbz(255) - 95.0).abs() < 1e-9);
|
||
}
|
||
|
||
#[test]
|
||
fn latlon_to_pixel_reasonable() {
|
||
// Corner sanity: (50N, -126W) → (0, 0)
|
||
assert_eq!(latlon_to_pixel(50.0, -126.0), (0, 0));
|
||
// Mid-CONUS: (32.9, -97.0)
|
||
let (x, y) = latlon_to_pixel(32.9, -97.0);
|
||
assert!(x > 5000 && x < 7000, "x={x}");
|
||
assert!(y > 2000 && y < 4000, "y={y}");
|
||
}
|
||
|
||
#[test]
|
||
fn max_dbz_picks_hottest_pixel_in_box() {
|
||
// 10x10 image. Non-zero pixel at (5, 5) = 200 (→ ~68 dBZ).
|
||
let width = 10;
|
||
let mut pixels = vec![0u8; 100];
|
||
pixels[5 * 10 + 5] = 200;
|
||
let dbz = max_dbz_in_box(&pixels, width, 5, 5, 2);
|
||
assert!((dbz - pixel_to_dbz(200)).abs() < 1e-9);
|
||
}
|
||
|
||
#[test]
|
||
fn box_clipped_at_edges() {
|
||
let width = 10;
|
||
let pixels = vec![100u8; 100];
|
||
// Center outside image — clipped to 0.
|
||
let dbz = max_dbz_in_box(&pixels, width, 50, 50, 5);
|
||
assert_eq!(dbz, 0.0);
|
||
}
|
||
}
|