src/gribnc.rs
//! **NetCDF/GRIB weather-model collector** — Copernicus Climate Data Store (CDS)
//! + Deutscher Wetterdienst (DWD) ICON open data. Numerical-weather-prediction
//! and reanalysis fields ship as **GRIB2** (DWD ICON, CDS default) or **NetCDF**
//! (CDS `format: netcdf`), two dense binary array containers.
//!
//! ## What this module lands: metadata CATALOG + grid-VALUE array unpack
//!
//! **Catalog (metadata index).** A pure-Rust, zero-dependency decoder walks the
//! GRIB2 message SECTION HEADERS (Section 0 indicator, Section 1 reference time,
//! Section 3 grid definition, Section 4 product definition) and the NetCDF
//! **classic** header (dimensions + variables) and emits one catalog row per GRIB
//! message / NetCDF variable (`format, dataset, variable, reference_time, ni, nj,
//! num_points, lat/lon bbox`) — the honest, testable seam a warehouse joins on to
//! know **what** a retrieval contains before paying to unpack it. See
//! [`parse_weather_bytes`].
//!
//! **Values (data array).** [`parse_weather_values`] / [`decode_grib2_values`] /
//! [`decode_netcdf_variable`] now unpack the ACTUAL gridded values, still pure-Rust
//! and zero-dep (matching the mainframe COBOL-copybook and gtfs-rt protobuf
//! hand-rolled binary readers). The GRIB2 Section 5/6/7 decode covers:
//! * **5.0 simple packing** — `Y = (R + X·2^E)/10^D` per grid point;
//! * **5.2 complex packing** — per-group reference / width / length;
//! * **5.3 complex packing + spatial differencing** (order 1 & 2) — the group
//! data are differences, un-differenced by cumulative sum from the Section-7
//! initial value(s) + overall minimum descriptors;
//! * **Section 6 bitmap** — a 0-bit is a MISSING point (`None`, never a fake 0.0).
//!
//! The NetCDF-classic path reads a variable's data array (fixed contiguous OR
//! record/unlimited-dim, numeric types byte/short/int/float/double) from its
//! header `begin`/`vsize`. All value paths are round-trip tested against real
//! packed byte streams built by the in-module encoders (no committed blob, no
//! network). The collector emits values under `params.mode = "values"`.
//!
//! **Still deferred:** NetCDF-4 / HDF5 (the `\x89HDF` container) array + header
//! decode — a native `libnetcdf` (C) path would go behind an OPTIONAL `netcdf`
//! feature (documented, never default), only if a pure-Rust reader proves
//! insufficient. It is an explicit, loud error here, never a silent empty result.
use anyhow::{Result, anyhow};
use arrow_array::RecordBatch;
use arrow_schema::{DataType, Schema};
use super::arrow_util::{batch, f64_col, i64_col, schema, str_col};
use super::{
AuthConfig, AuthKind, BatchSink, Bounds, CollectStats, Collector, CollectorCfg, CollectorMeta,
CollectorPath, Direction, GdprClass, SourceEndpoint, SourceKind, SourceTarget,
};
// ═══════════════════════════════════════════════════════════════════════════
// Unified metadata-catalog row (both GRIB messages and NetCDF variables fill it)
// ═══════════════════════════════════════════════════════════════════════════
/// One row of the weather-model **metadata catalog** — a GRIB2 message or a
/// NetCDF variable, described by what its header alone reveals (the array is not
/// decoded this pass). `None` fields are simply not encoded by that container's
/// header (a NetCDF classic header carries no per-variable geographic bbox).
#[derive(Debug, Clone, PartialEq)]
pub struct CatalogRow {
/// `"grib2"` or `"netcdf3"` — which container this row came from.
pub format: String,
/// The retrieval/dataset label (from `cfg.params.dataset` / the endpoint).
pub dataset: String,
/// GRIB: a `discipline/category/number` parameter label. NetCDF: the variable name.
pub variable: String,
/// GRIB: the Section-1 reference time (ISO-8601). NetCDF classic: `None`.
pub reference_time: Option<String>,
/// Grid width (GRIB Ni; NetCDF: the last/innermost dimension length).
pub ni: Option<i64>,
/// Grid height (GRIB Nj; NetCDF: the second-innermost dimension length).
pub nj: Option<i64>,
/// Total data points (GRIB number-of-data-points; NetCDF: product of var dims).
pub num_points: Option<i64>,
pub lat_min: Option<f64>,
pub lat_max: Option<f64>,
pub lon_min: Option<f64>,
pub lon_max: Option<f64>,
/// A short free-form note (GRIB grid-template #; NetCDF nc_type + dim names).
pub detail: String,
}
/// The Arrow schema of the metadata catalog — the table both collectors write.
pub fn catalog_schema() -> Schema {
schema(&[
("format", DataType::Utf8, false),
("dataset", DataType::Utf8, false),
("variable", DataType::Utf8, false),
("reference_time", DataType::Utf8, true),
("ni", DataType::Int64, true),
("nj", DataType::Int64, true),
("num_points", DataType::Int64, true),
("lat_min", DataType::Float64, true),
("lat_max", DataType::Float64, true),
("lon_min", DataType::Float64, true),
("lon_max", DataType::Float64, true),
("detail", DataType::Utf8, false),
])
}
/// Build the catalog `RecordBatch` from decoded rows.
pub fn catalog_batch(rows: &[CatalogRow]) -> Result<RecordBatch> {
let cols = vec![
str_col(
&rows
.iter()
.map(|r| Some(r.format.clone()))
.collect::<Vec<_>>(),
),
str_col(
&rows
.iter()
.map(|r| Some(r.dataset.clone()))
.collect::<Vec<_>>(),
),
str_col(
&rows
.iter()
.map(|r| Some(r.variable.clone()))
.collect::<Vec<_>>(),
),
str_col(
&rows
.iter()
.map(|r| r.reference_time.clone())
.collect::<Vec<_>>(),
),
i64_col(&rows.iter().map(|r| r.ni).collect::<Vec<_>>()),
i64_col(&rows.iter().map(|r| r.nj).collect::<Vec<_>>()),
i64_col(&rows.iter().map(|r| r.num_points).collect::<Vec<_>>()),
f64_col(&rows.iter().map(|r| r.lat_min).collect::<Vec<_>>()),
f64_col(&rows.iter().map(|r| r.lat_max).collect::<Vec<_>>()),
f64_col(&rows.iter().map(|r| r.lon_min).collect::<Vec<_>>()),
f64_col(&rows.iter().map(|r| r.lon_max).collect::<Vec<_>>()),
str_col(
&rows
.iter()
.map(|r| Some(r.detail.clone()))
.collect::<Vec<_>>(),
),
];
batch(&catalog_schema(), cols)
}
/// Detect the container by magic and decode its header → catalog rows.
/// `\x89HDF` (NetCDF-4) is recognised but deferred (a documented error, not a
/// silent empty result).
pub fn parse_weather_bytes(bytes: &[u8], dataset: &str) -> Result<Vec<CatalogRow>> {
if bytes.len() >= 4 && &bytes[0..4] == b"GRIB" {
let msgs = parse_grib2(bytes)?;
Ok(grib_to_rows(&msgs, dataset))
} else if bytes.len() >= 3 && &bytes[0..3] == b"CDF" {
let hdr = parse_netcdf_classic_header(bytes)?;
Ok(netcdf_to_rows(&hdr, dataset))
} else if bytes.len() >= 4 && &bytes[0..4] == b"\x89HDF" {
Err(anyhow!(
"NetCDF-4/HDF5 container: array + header decode DEFERRED this pass (would go behind \
the optional native `netcdf` feature)"
))
} else {
Err(anyhow!(
"unrecognised weather container (not GRIB2, NetCDF-classic or HDF5)"
))
}
}
/// Unpack the actual gridded VALUES from a weather container (the array-decode
/// companion to the metadata-only [`parse_weather_bytes`]). Detects the container
/// by magic and yields a typed [`values_schema`] `RecordBatch`:
/// * **GRIB2** → every message's grid via [`decode_grib2_values`] (simple 5.0 +
/// complex 5.2/5.3 packing, Section-6 bitmap honoured), one row per point;
/// * **NetCDF-classic** → the named `variable` (or the first numeric variable)
/// via [`decode_netcdf_variable`], one row per element.
///
/// NetCDF-4/HDF5 stays deferred (documented, never a silent empty result).
pub fn parse_weather_values(bytes: &[u8], variable: Option<&str>) -> Result<RecordBatch> {
if bytes.len() >= 4 && &bytes[0..4] == b"GRIB" {
let vals = decode_grib2_values(bytes)?;
grib_values_batch(&vals)
} else if bytes.len() >= 3 && &bytes[0..3] == b"CDF" {
let hdr = parse_netcdf_classic_header(bytes)?;
let name = match variable {
Some(v) => v.to_string(),
None => hdr
.vars
.iter()
.find(|v| nc_esize(v.nc_type).is_ok())
.map(|v| v.name.clone())
.ok_or_else(|| anyhow!("NetCDF file has no numeric variable to unpack"))?,
};
netcdf_variable_batch(bytes, &hdr, &name)
} else if bytes.len() >= 4 && &bytes[0..4] == b"\x89HDF" {
Err(anyhow!(
"NetCDF-4/HDF5 array decode DEFERRED this pass (optional native `netcdf` feature)"
))
} else {
Err(anyhow!(
"unrecognised weather container (not GRIB2, NetCDF-classic or HDF5)"
))
}
}
// ═══════════════════════════════════════════════════════════════════════════
// GRIB2 — message-header (section) decoder
// ═══════════════════════════════════════════════════════════════════════════
/// The header-level metadata of one GRIB2 message (the data array is not decoded).
#[derive(Debug, Clone, PartialEq)]
pub struct GribMessage {
pub discipline: u8,
pub parameter_category: u8,
pub parameter_number: u8,
/// ISO-8601 reference time from Section 1.
pub reference_time: String,
/// Grid definition template number (Section 3) — `0` = lat/lon (equidistant).
pub grid_template: u16,
pub num_data_points: u32,
pub ni: u32,
pub nj: u32,
/// Bounding box in degrees (from La1/Lo1/La2/Lo2, GRIB micro-degree units).
pub lat_first: f64,
pub lon_first: f64,
pub lat_last: f64,
pub lon_last: f64,
}
fn be_u16(b: &[u8]) -> u16 {
u16::from_be_bytes([b[0], b[1]])
}
fn be_u32(b: &[u8]) -> u32 {
u32::from_be_bytes([b[0], b[1], b[2], b[3]])
}
fn be_u64(b: &[u8]) -> u64 {
u64::from_be_bytes([b[0], b[1], b[2], b[3], b[4], b[5], b[6], b[7]])
}
/// GRIB2 uses **sign-magnitude** for signed lat/lon (the high bit is the sign,
/// NOT two's complement).
fn grib_i32(b: &[u8]) -> i64 {
let raw = be_u32(b);
if raw & 0x8000_0000 != 0 {
-((raw & 0x7FFF_FFFF) as i64)
} else {
raw as i64
}
}
/// One fully-walked GRIB2 message: the decoded [`GribMessage`] header PLUS the
/// raw Section 5 (data representation), Section 6 (bitmap) and Section 7 (data)
/// bytes needed to unpack the grid VALUES. Produced by [`walk_grib2`]; the
/// header-only [`parse_grib2`] and the value-unpacking [`decode_grib2_values`]
/// both build on it.
#[derive(Debug, Clone)]
struct RawMessage {
header: GribMessage,
/// Section-3 total grid points (the value grid size incl. any bitmapped-out).
num_points: u32,
/// Section-5 data-representation (packing) parameters, when present.
sec5: Option<Sec5>,
/// Section-6 bitmap indicator (0 = bitmap present, 255 = none).
bitmap_indicator: u8,
/// Section-6 bitmap bytes (1 bit / grid point, MSB-first) when indicator is 0.
bitmap: Vec<u8>,
/// Section-7 packed data bytes (the payload the packing template unpacks).
data: Vec<u8>,
}
/// Walk every GRIB2 message in `bytes`, capturing each one's header + the raw
/// Section 5/6/7 bytes. The shared spine of [`parse_grib2`] (header only) and
/// [`decode_grib2_values`] (header + unpacked grid). Walks Section 0 → the
/// intra-message section chain (dispatching Sections 1/3/4/5/6/7, skipping the
/// rest by length) → the `7777` end marker, then advances to the next message.
fn walk_grib2(bytes: &[u8]) -> Result<Vec<RawMessage>> {
let mut out = Vec::new();
let mut base = 0usize;
while base + 16 <= bytes.len() {
if &bytes[base..base + 4] != b"GRIB" {
// Not the start of a message — stop (trailing padding / EOF).
break;
}
let discipline = bytes[base + 6];
let edition = bytes[base + 7];
if edition != 2 {
return Err(anyhow!("GRIB edition {edition} not supported (only GRIB2)"));
}
let total_len = be_u64(&bytes[base + 8..base + 16]) as usize;
if total_len < 16 || base + total_len > bytes.len() {
return Err(anyhow!(
"GRIB2 message total length {total_len} runs past the buffer"
));
}
let msg_end = base + total_len;
// Walk the section chain from just after Section 0.
let mut pos = base + 16;
let mut reference_time = String::new();
let mut grid_template = 0u16;
let mut num_data_points = 0u32;
let (mut ni, mut nj) = (0u32, 0u32);
let (mut la1, mut lo1, mut la2, mut lo2) = (0f64, 0f64, 0f64, 0f64);
let mut parameter_category = 0u8;
let mut parameter_number = 0u8;
let mut sec5: Option<Sec5> = None;
let mut bitmap_indicator = 255u8;
let mut bitmap: Vec<u8> = Vec::new();
let mut data: Vec<u8> = Vec::new();
while pos + 4 <= msg_end {
if &bytes[pos..pos + 4] == b"7777" {
break; // Section 8 (end marker)
}
if pos + 5 > msg_end {
return Err(anyhow!("GRIB2 truncated section header"));
}
let sec_len = be_u32(&bytes[pos..pos + 4]) as usize;
let sec_num = bytes[pos + 4];
if sec_len < 5 || pos + sec_len > msg_end {
return Err(anyhow!(
"GRIB2 section {sec_num} length {sec_len} out of bounds"
));
}
let sec = &bytes[pos..pos + sec_len];
match sec_num {
1 if sec_len >= 21 => {
let year = be_u16(&sec[12..14]);
let (mo, da, hr, mi, se) = (sec[14], sec[15], sec[16], sec[17], sec[18]);
reference_time = format!("{year:04}-{mo:02}-{da:02}T{hr:02}:{mi:02}:{se:02}Z");
}
3 if sec_len >= 72 => {
num_data_points = be_u32(&sec[6..10]);
grid_template = be_u16(&sec[12..14]);
// Template 3.0 (lat/lon) field offsets within the section.
ni = be_u32(&sec[30..34]);
nj = be_u32(&sec[34..38]);
la1 = grib_i32(&sec[46..50]) as f64 / 1e6;
lo1 = be_u32(&sec[50..54]) as f64 / 1e6;
la2 = grib_i32(&sec[55..59]) as f64 / 1e6;
lo2 = be_u32(&sec[59..63]) as f64 / 1e6;
}
4 if sec_len >= 11 => {
// Only trust the parameter octets for product template 4.0/4.x
// that share the category/number at offsets 9/10.
parameter_category = sec[9];
parameter_number = sec[10];
}
5 => sec5 = Some(parse_sec5(sec)?),
6 if sec_len >= 6 => {
bitmap_indicator = sec[5];
if bitmap_indicator == 0 {
bitmap = sec[6..].to_vec();
}
}
7 if sec_len >= 5 => data = sec[5..].to_vec(),
_ => {}
}
pos += sec_len;
}
out.push(RawMessage {
header: GribMessage {
discipline,
parameter_category,
parameter_number,
reference_time,
grid_template,
num_data_points,
ni,
nj,
lat_first: la1,
lon_first: lo1,
lat_last: la2,
lon_last: lo2,
},
num_points: num_data_points,
sec5,
bitmap_indicator,
bitmap,
data,
});
base = msg_end;
}
if out.is_empty() {
return Err(anyhow!("no GRIB2 message found (missing `GRIB` indicator)"));
}
Ok(out)
}
/// Decode every GRIB2 message HEADER in `bytes` (a `.grib`/`.grib2` file may
/// concatenate several). The DATA sections are captured but not unpacked here —
/// use [`decode_grib2_values`] to unpack the packed grid values.
pub fn parse_grib2(bytes: &[u8]) -> Result<Vec<GribMessage>> {
Ok(walk_grib2(bytes)?.into_iter().map(|m| m.header).collect())
}
/// Map decoded GRIB messages → catalog rows (min/max normalise the bbox so
/// scanning-order doesn't flip it).
pub fn grib_to_rows(msgs: &[GribMessage], dataset: &str) -> Vec<CatalogRow> {
msgs.iter()
.map(|m| CatalogRow {
format: "grib2".into(),
dataset: dataset.into(),
variable: format!(
"d{}.c{}.n{}",
m.discipline, m.parameter_category, m.parameter_number
),
reference_time: (!m.reference_time.is_empty()).then(|| m.reference_time.clone()),
ni: Some(m.ni as i64),
nj: Some(m.nj as i64),
num_points: Some(m.num_data_points as i64),
lat_min: Some(m.lat_first.min(m.lat_last)),
lat_max: Some(m.lat_first.max(m.lat_last)),
lon_min: Some(m.lon_first.min(m.lon_last)),
lon_max: Some(m.lon_first.max(m.lon_last)),
detail: format!("grid_template={}", m.grid_template),
})
.collect()
}
// ═══════════════════════════════════════════════════════════════════════════
// GRIB2 — Section 5/6/7 data-VALUE unpack (simple + complex packing)
// ═══════════════════════════════════════════════════════════════════════════
//
// The header walk above indexes WHAT a message contains; this unpacks the actual
// packed grid values. Three data-representation templates are decoded, pure-Rust,
// zero-dep (matching the mainframe/gtfs-rt hand-rolled binary readers):
//
// * **5.0 simple packing** — each grid point is one `nbits`-wide integer `X`;
// `Y = (R + X·2^E)/10^D` (reference `R`, binary scale `E`, decimal scale `D`).
// * **5.2 complex packing** — points are split into groups; each group carries a
// reference, a bit-width and a length, and every point is `group_ref + X`.
// * **5.3 complex packing + spatial differencing** — as 5.2, but the group data
// are 1st/2nd-order differences (with the initial values + overall minimum in
// Section-7 extra descriptors), un-differenced by cumulative sum.
//
// Section 6 (bitmap) is honoured: a 0-bit marks a MISSING grid point (decoded as
// `None`, never a fabricated 0.0). NetCDF-4/HDF5 packing stays deferred (§ module).
/// Section-5 (data representation) parameters recovered for the value unpack.
/// The simple-packing fields (`template`..`bits_per_value`) are always present;
/// the complex-packing fields are only meaningful for templates 5.2 / 5.3.
#[derive(Debug, Clone, PartialEq)]
pub struct Sec5 {
/// Data-representation template number (0 = simple, 2 = complex, 3 = complex+SD).
pub template: u16,
/// Reference value `R` (IEEE-754 32-bit) — the field's packing offset.
pub reference_value: f32,
/// Binary scale factor `E` (signed; sign-magnitude on the wire).
pub binary_scale: i32,
/// Decimal scale factor `D` (signed; sign-magnitude on the wire).
pub decimal_scale: i32,
/// Bits per packed value (simple: per grid point; complex: per group reference).
pub bits_per_value: u8,
// ── complex packing (templates 5.2 / 5.3) ──────────────────────────────
/// Number of groups (`NG`).
pub ng: u32,
/// Reference added to every group width.
pub ref_group_width: u8,
/// Bits used for each (scaled) group width.
pub bits_group_width: u8,
/// Reference added to every group length.
pub ref_group_length: u32,
/// Length increment the scaled group lengths are multiplied by.
pub length_increment: u8,
/// True length of the last group (overrides its scaled length).
pub true_last_group_length: u32,
/// Bits used for each scaled group length.
pub bits_scaled_group_length: u8,
// ── spatial differencing (template 5.3) ────────────────────────────────
/// Order of spatial differencing (1 or 2; 0 when not template 5.3).
pub spatial_order: u8,
/// Octets used for each Section-7 extra descriptor (initial values + minimum).
pub spatial_desc_octets: u8,
}
/// GRIB2 signed 2-octet field (binary/decimal scale) — **sign-magnitude** (the
/// high bit is the sign, NOT two's complement), like [`grib_i32`].
fn grib_i16(b: &[u8]) -> i32 {
let raw = be_u16(b);
if raw & 0x8000 != 0 {
-((raw & 0x7FFF) as i32)
} else {
raw as i32
}
}
/// A big-endian unsigned integer over `nbytes` (≤8) octets.
fn be_uint(b: &[u8]) -> u64 {
b.iter().fold(0u64, |a, &x| (a << 8) | x as u64)
}
/// A sign-magnitude big-endian integer over `b.len()` octets (GRIB's signed form).
fn sign_mag_int(b: &[u8]) -> i64 {
if b.is_empty() {
return 0;
}
let neg = b[0] & 0x80 != 0;
let mut v = (b[0] & 0x7F) as i64;
for &x in &b[1..] {
v = (v << 8) | x as i64;
}
if neg { -v } else { v }
}
/// Parse Section 5 (data representation) → [`Sec5`]. Bounds-checked per template.
fn parse_sec5(sec: &[u8]) -> Result<Sec5> {
if sec.len() < 21 {
return Err(anyhow!(
"GRIB2 Section 5 too short ({} octets, need ≥21)",
sec.len()
));
}
let template = be_u16(&sec[9..11]);
let reference_value = f32::from_be_bytes([sec[11], sec[12], sec[13], sec[14]]);
let mut s = Sec5 {
template,
reference_value,
binary_scale: grib_i16(&sec[15..17]),
decimal_scale: grib_i16(&sec[17..19]),
bits_per_value: sec[19],
ng: 0,
ref_group_width: 0,
bits_group_width: 0,
ref_group_length: 0,
length_increment: 0,
true_last_group_length: 0,
bits_scaled_group_length: 0,
spatial_order: 0,
spatial_desc_octets: 0,
};
if template == 2 || template == 3 {
if sec.len() < 47 {
return Err(anyhow!("GRIB2 Section 5 complex template needs ≥47 octets"));
}
s.ng = be_u32(&sec[31..35]);
s.ref_group_width = sec[35];
s.bits_group_width = sec[36];
s.ref_group_length = be_u32(&sec[37..41]);
s.length_increment = sec[41];
s.true_last_group_length = be_u32(&sec[42..46]);
s.bits_scaled_group_length = sec[46];
if template == 3 {
if sec.len() < 49 {
return Err(anyhow!(
"GRIB2 Section 5 spatial-diff template needs ≥49 octets"
));
}
s.spatial_order = sec[47];
s.spatial_desc_octets = sec[48];
}
}
Ok(s)
}
/// An MSB-first bit reader over a Section-7 payload (values pack big-endian, most
/// significant bit first, with no padding between values). Reading 0 bits → 0.
struct BitReader<'a> {
b: &'a [u8],
byte: usize,
bit: u8,
}
impl<'a> BitReader<'a> {
fn new(b: &'a [u8]) -> Self {
Self { b, byte: 0, bit: 0 }
}
fn read_bits(&mut self, n: u32) -> Result<u64> {
let mut v = 0u64;
for _ in 0..n {
if self.byte >= self.b.len() {
return Err(anyhow!("GRIB2 Section 7 bitstream underrun"));
}
let bitval = (self.b[self.byte] >> (7 - self.bit)) & 1;
v = (v << 1) | bitval as u64;
self.bit += 1;
if self.bit == 8 {
self.bit = 0;
self.byte += 1;
}
}
Ok(v)
}
}
/// An MSB-first bit writer (the [`BitReader`] inverse) — used by the in-module
/// encoders to build real packed Section-7 payloads for round-trip tests.
struct BitWriter {
out: Vec<u8>,
cur: u8,
nbit: u8,
}
impl BitWriter {
fn new() -> Self {
Self {
out: Vec::new(),
cur: 0,
nbit: 0,
}
}
fn write_bits(&mut self, val: u64, n: u32) {
for i in (0..n).rev() {
let bit = ((val >> i) & 1) as u8;
self.cur = (self.cur << 1) | bit;
self.nbit += 1;
if self.nbit == 8 {
self.out.push(self.cur);
self.cur = 0;
self.nbit = 0;
}
}
}
fn finish(mut self) -> Vec<u8> {
if self.nbit > 0 {
self.cur <<= 8 - self.nbit; // pad the final octet with 0 bits
self.out.push(self.cur);
}
self.out
}
}
/// Apply the Section-5 scaling `Y = (R + X·2^E) / 10^D` to a packed integer `X`.
fn apply_scale(s5: &Sec5, x: f64) -> f64 {
(s5.reference_value as f64 + x * 2f64.powi(s5.binary_scale)) / 10f64.powi(s5.decimal_scale)
}
/// The present-point mask from Section 6: `255` = no bitmap (all present), `0` =
/// an explicit 1-bit-per-point bitmap (1 = present, 0 = missing). A pre-defined
/// bitmap (any other indicator) is an explicit unsupported error, not a silent gap.
fn present_mask(indicator: u8, bitmap: &[u8], num_points: u32) -> Result<Vec<bool>> {
match indicator {
255 => Ok(vec![true; num_points as usize]),
0 => {
let mut br = BitReader::new(bitmap);
let mut m = Vec::with_capacity(num_points as usize);
for _ in 0..num_points {
m.push(br.read_bits(1)? == 1);
}
Ok(m)
}
other => Err(anyhow!(
"GRIB2 bitmap indicator {other} (pre-defined bitmap) unsupported"
)),
}
}
/// Decode simple packing (template 5.0): `npresent` values, each `bits_per_value`
/// wide, then the `(R + X·2^E)/10^D` scale. `bits_per_value == 0` = a constant
/// field (every present point is `R/10^D`).
fn decode_simple(s5: &Sec5, data: &[u8], npresent: usize) -> Result<Vec<f64>> {
let mut br = BitReader::new(data);
let nbits = s5.bits_per_value as u32;
let mut out = Vec::with_capacity(npresent);
for _ in 0..npresent {
let x = if nbits == 0 { 0 } else { br.read_bits(nbits)? };
out.push(apply_scale(s5, x as f64));
}
Ok(out)
}
/// Decode complex packing (template 5.2) and complex + spatial differencing
/// (5.3). The group structure (per-group reference / width / length) is read from
/// the Section-7 bitstream; with `spatial`, the leading octet-aligned extra
/// descriptors (initial value(s) + overall minimum) drive a cumulative-sum
/// reconstruction of the un-differenced field.
fn decode_complex(s5: &Sec5, data: &[u8], npresent: usize, spatial: bool) -> Result<Vec<f64>> {
let ng = s5.ng as usize;
if ng == 0 {
return Err(anyhow!("GRIB2 complex packing with 0 groups"));
}
// Spatial-differencing extra descriptors (octet-aligned at the data start).
let (mut ival, mut gmin, mut off) = (Vec::<i64>::new(), 0i64, 0usize);
if spatial {
let nb = s5.spatial_desc_octets as usize;
let order = s5.spatial_order as usize;
if nb == 0 {
return Err(anyhow!(
"GRIB2 spatial differencing with 0 descriptor octets"
));
}
if order == 0 || order > 2 {
return Err(anyhow!(
"GRIB2 spatial-diff order {} unsupported (1 or 2)",
s5.spatial_order
));
}
if data.len() < nb * (order + 1) {
return Err(anyhow!(
"GRIB2 Section 7 too short for spatial-diff descriptors"
));
}
for k in 0..order {
ival.push(be_uint(&data[k * nb..(k + 1) * nb]) as i64);
}
gmin = sign_mag_int(&data[order * nb..(order + 1) * nb]);
off = nb * (order + 1);
}
let mut br = BitReader::new(&data[off..]);
// Group references, then widths, then (scaled) lengths — three parallel runs.
let mut gref = Vec::with_capacity(ng);
for _ in 0..ng {
gref.push(br.read_bits(s5.bits_per_value as u32)? as i64);
}
let mut width = Vec::with_capacity(ng);
for _ in 0..ng {
width.push(s5.ref_group_width as u32 + br.read_bits(s5.bits_group_width as u32)? as u32);
}
let mut length = Vec::with_capacity(ng);
for _ in 0..ng {
let scaled = br.read_bits(s5.bits_scaled_group_length as u32)? as u32;
length.push(s5.ref_group_length + scaled * s5.length_increment as u32);
}
if let Some(last) = length.last_mut() {
*last = s5.true_last_group_length;
}
// Group data: each point is group_ref + a width-wide offset.
let mut vals: Vec<i64> = Vec::with_capacity(npresent);
for k in 0..ng {
for _ in 0..length[k] {
let x = if width[k] == 0 {
0
} else {
br.read_bits(width[k])? as i64
};
vals.push(gref[k] + x);
}
}
// Un-difference (template 5.3): restore the field by cumulative sum.
if spatial {
let order = s5.spatial_order as usize;
for v in vals.iter_mut() {
*v += gmin;
}
if order >= 1 && !vals.is_empty() {
vals[0] = ival[0];
}
if order == 2 && vals.len() >= 2 {
vals[1] = ival[1];
}
match order {
1 => {
for i in 1..vals.len() {
vals[i] += vals[i - 1];
}
}
2 => {
for i in 2..vals.len() {
vals[i] += 2 * vals[i - 1] - vals[i - 2];
}
}
_ => {}
}
}
if vals.len() != npresent {
return Err(anyhow!(
"GRIB2 complex packing decoded {} values, expected {} present points",
vals.len(),
npresent
));
}
Ok(vals
.into_iter()
.map(|x| apply_scale(s5, x as f64))
.collect())
}
/// One GRIB2 message's unpacked grid: the decoded [`GribMessage`] header + the
/// full value vector (length = the Section-3 grid size). A `None` entry is a
/// MISSING point (a Section-6 bitmap 0-bit), never a fabricated 0.0.
#[derive(Debug, Clone)]
pub struct GribValues {
pub header: GribMessage,
pub values: Vec<Option<f64>>,
}
impl GribValues {
/// Count / min / max / mean over the PRESENT (non-missing) points — a cheap
/// sanity summary the collector and tests assert against.
pub fn stats(&self) -> (usize, Option<f64>, Option<f64>, Option<f64>) {
let present: Vec<f64> = self.values.iter().flatten().copied().collect();
if present.is_empty() {
return (0, None, None, None);
}
let min = present.iter().copied().fold(f64::INFINITY, f64::min);
let max = present.iter().copied().fold(f64::NEG_INFINITY, f64::max);
let mean = present.iter().sum::<f64>() / present.len() as f64;
(present.len(), Some(min), Some(max), Some(mean))
}
}
/// Unpack every GRIB2 message's grid VALUES in `bytes` (all present templates:
/// 5.0 simple, 5.2 complex, 5.3 complex + spatial differencing; Section-6 bitmap
/// honoured). This is the array-decode companion to the header-only [`parse_grib2`].
pub fn decode_grib2_values(bytes: &[u8]) -> Result<Vec<GribValues>> {
walk_grib2(bytes)?.into_iter().map(decode_raw).collect()
}
fn decode_raw(m: RawMessage) -> Result<GribValues> {
let np = m.num_points as usize;
let sec5 = m
.sec5
.ok_or_else(|| anyhow!("GRIB2 message has no Section 5 (cannot unpack values)"))?;
let mask = present_mask(m.bitmap_indicator, &m.bitmap, m.num_points)?;
let npresent = mask.iter().filter(|&&p| p).count();
let packed: Vec<f64> = match sec5.template {
0 => decode_simple(&sec5, &m.data, npresent)?,
2 => decode_complex(&sec5, &m.data, npresent, false)?,
3 => decode_complex(&sec5, &m.data, npresent, true)?,
t => {
return Err(anyhow!(
"GRIB2 data representation template 5.{t} unsupported (only 5.0 simple, \
5.2/5.3 complex)"
));
}
};
// Scatter present values back into the full grid (absent point = missing).
let mut values = Vec::with_capacity(np);
let mut it = packed.into_iter();
for present in &mask {
values.push(if *present { it.next() } else { None });
}
Ok(GribValues {
header: m.header,
values,
})
}
/// The lat/lon of grid point `p` for a regular lat/lon grid (template 3.0),
/// linearly interpolated across the La1/Lo1 → La2/Lo2 corners (row-major, `ni`
/// wide). `None` for any other grid template (no assumed geometry).
fn grid_latlon(h: &GribMessage, p: i64, ni: i64, nj: i64) -> (Option<f64>, Option<f64>) {
if h.grid_template != 0 || ni <= 0 || nj <= 0 {
return (None, None);
}
let (row, col) = (p / ni, p % ni);
let lat = if nj > 1 {
h.lat_first + (h.lat_last - h.lat_first) * (row as f64) / ((nj - 1) as f64)
} else {
h.lat_first
};
let lon = if ni > 1 {
h.lon_first + (h.lon_last - h.lon_first) * (col as f64) / ((ni - 1) as f64)
} else {
h.lon_first
};
(Some(lat), Some(lon))
}
/// The Arrow schema of the unpacked-VALUES table (one row per grid point).
pub fn values_schema() -> Schema {
schema(&[
("message", DataType::Int64, false),
("variable", DataType::Utf8, false),
("point", DataType::Int64, false),
("lat", DataType::Float64, true),
("lon", DataType::Float64, true),
("value", DataType::Float64, true),
])
}
/// Build the unpacked-VALUES `RecordBatch` — one row per grid point across every
/// decoded message: `(message index, variable, point index, lat, lon, value)`.
/// Missing points land as a null `value` (bitmap-honoured); lat/lon are filled
/// for regular lat/lon grids, null otherwise.
pub fn grib_values_batch(msgs: &[GribValues]) -> Result<RecordBatch> {
let (mut m_i, mut var, mut pt) = (Vec::new(), Vec::new(), Vec::new());
let (mut lat, mut lon, mut val) = (Vec::new(), Vec::new(), Vec::new());
for (mi, gv) in msgs.iter().enumerate() {
let h = &gv.header;
let vname = format!(
"d{}.c{}.n{}",
h.discipline, h.parameter_category, h.parameter_number
);
let (ni, nj) = (h.ni as i64, h.nj as i64);
for (p, v) in gv.values.iter().enumerate() {
m_i.push(Some(mi as i64));
var.push(Some(vname.clone()));
pt.push(Some(p as i64));
let (la, lo) = grid_latlon(h, p as i64, ni, nj);
lat.push(la);
lon.push(lo);
val.push(*v);
}
}
batch(
&values_schema(),
vec![
i64_col(&m_i),
str_col(&var),
i64_col(&pt),
f64_col(&lat),
f64_col(&lon),
f64_col(&val),
],
)
}
// ═══════════════════════════════════════════════════════════════════════════
// NetCDF classic (CDF-1/CDF-2) — header decoder (dimensions + variables)
// ═══════════════════════════════════════════════════════════════════════════
const NC_DIMENSION: u32 = 0x0A;
const NC_VARIABLE: u32 = 0x0B;
const NC_ATTRIBUTE: u32 = 0x0C;
/// A NetCDF classic dimension (name + length; a length of 0 = the record dim).
#[derive(Debug, Clone, PartialEq)]
pub struct NetcdfDim {
pub name: String,
pub length: u64,
}
/// A NetCDF classic variable header (name, nc_type, the dimension ids it spans,
/// and the on-disk layout — `vsize` bytes at file offset `begin`) needed to read
/// its data array.
#[derive(Debug, Clone, PartialEq)]
pub struct NetcdfVar {
pub name: String,
pub nc_type: u32,
pub dim_ids: Vec<u32>,
/// The variable's size in bytes (one record's slab for a record variable),
/// padded to a 4-byte boundary — the record stride for the record section.
pub vsize: u64,
/// The byte offset of the variable's data in the file (`begin` field).
pub begin: u64,
}
/// The decoded NetCDF classic header (no data arrays).
#[derive(Debug, Clone, PartialEq)]
pub struct NetcdfHeader {
pub num_records: u32,
pub dims: Vec<NetcdfDim>,
pub vars: Vec<NetcdfVar>,
}
/// The nc_type → human label.
fn nc_type_name(t: u32) -> &'static str {
match t {
1 => "byte",
2 => "char",
3 => "short",
4 => "int",
5 => "float",
6 => "double",
_ => "unknown",
}
}
/// A minimal big-endian cursor over the NetCDF header bytes.
struct Cur<'a> {
b: &'a [u8],
p: usize,
}
impl<'a> Cur<'a> {
fn u32(&mut self) -> Result<u32> {
if self.p + 4 > self.b.len() {
return Err(anyhow!("NetCDF header truncated (u32 at {})", self.p));
}
let v = be_u32(&self.b[self.p..self.p + 4]);
self.p += 4;
Ok(v)
}
/// A big-endian unsigned integer over `nbytes` (4 for CDF-1 `begin`, 8 for CDF-2).
fn uint(&mut self, nbytes: usize) -> Result<u64> {
if self.p + nbytes > self.b.len() {
return Err(anyhow!(
"NetCDF header truncated ({nbytes}-octet int at {})",
self.p
));
}
let v = be_uint(&self.b[self.p..self.p + nbytes]);
self.p += nbytes;
Ok(v)
}
/// A netCDF name/string: `nelems` (u32) + bytes, padded to a 4-byte boundary.
fn name(&mut self) -> Result<String> {
let n = self.u32()? as usize;
if self.p + n > self.b.len() {
return Err(anyhow!("NetCDF string of {n} bytes runs past the header"));
}
let s = String::from_utf8_lossy(&self.b[self.p..self.p + n]).into_owned();
let padded = n.div_ceil(4) * 4;
self.p += padded;
Ok(s)
}
fn skip(&mut self, n: usize) -> Result<()> {
if self.p + n > self.b.len() {
return Err(anyhow!("NetCDF header skip of {n} runs past the buffer"));
}
self.p += n;
Ok(())
}
}
/// Consume an attribute list (`ABSENT | NC_ATTRIBUTE nelems att*`) — we don't
/// keep the values, but we must walk them to reach the next list.
fn skip_att_list(c: &mut Cur) -> Result<()> {
let tag = c.u32()?;
let nelems = c.u32()?;
if tag == 0 && nelems == 0 {
return Ok(()); // ABSENT
}
if tag != NC_ATTRIBUTE {
return Err(anyhow!("expected NC_ATTRIBUTE (0x0C), got {tag:#x}"));
}
for _ in 0..nelems {
let _name = c.name()?;
let nc_type = c.u32()?;
let n = c.u32()? as usize;
let elem = match nc_type {
1 | 2 => 1,
3 => 2,
4 | 5 => 4,
6 => 8,
other => return Err(anyhow!("attribute nc_type {other} unknown")),
};
let bytes = n * elem;
c.skip(bytes.div_ceil(4) * 4)?; // values, padded to 4-byte boundary
}
Ok(())
}
/// Decode a NetCDF **classic** (CDF-1/CDF-2) header: the dimension + variable
/// lists (attribute lists are walked but not retained). NetCDF-4/HDF5 is a
/// different container (see [`parse_weather_bytes`]).
pub fn parse_netcdf_classic_header(bytes: &[u8]) -> Result<NetcdfHeader> {
if bytes.len() < 4 || &bytes[0..3] != b"CDF" {
return Err(anyhow!("not a NetCDF classic file (missing `CDF` magic)"));
}
let version = bytes[3];
if version != 1 && version != 2 {
return Err(anyhow!(
"NetCDF classic version {version} unsupported (only CDF-1/CDF-2)"
));
}
// offset_size: begin field is u32 in CDF-1, u64 in CDF-2.
let offset_size = if version == 2 { 8 } else { 4 };
let mut c = Cur { b: bytes, p: 4 };
let num_records = c.u32()?;
// dim_list = ABSENT | NC_DIMENSION nelems dim*
let mut dims = Vec::new();
let tag = c.u32()?;
let ndims = c.u32()?;
if tag == NC_DIMENSION {
for _ in 0..ndims {
let name = c.name()?;
let length = c.u32()? as u64;
dims.push(NetcdfDim { name, length });
}
} else if !(tag == 0 && ndims == 0) {
return Err(anyhow!(
"expected NC_DIMENSION (0x0A) or ABSENT, got {tag:#x}"
));
}
// gatt_list (global attributes) — walked, not retained.
skip_att_list(&mut c)?;
// var_list = ABSENT | NC_VARIABLE nelems var*
let mut vars = Vec::new();
let tag = c.u32()?;
let nvars = c.u32()?;
if tag == NC_VARIABLE {
for _ in 0..nvars {
let name = c.name()?;
let nvdims = c.u32()? as usize;
let mut dim_ids = Vec::with_capacity(nvdims);
for _ in 0..nvdims {
dim_ids.push(c.u32()?);
}
skip_att_list(&mut c)?; // per-variable attributes
let nc_type = c.u32()?;
let vsize = c.u32()? as u64;
let begin = c.uint(offset_size)?; // begin (u32 CDF-1 / u64 CDF-2)
vars.push(NetcdfVar {
name,
nc_type,
dim_ids,
vsize,
begin,
});
}
} else if !(tag == 0 && nvars == 0) {
return Err(anyhow!(
"expected NC_VARIABLE (0x0B) or ABSENT, got {tag:#x}"
));
}
Ok(NetcdfHeader {
num_records,
dims,
vars,
})
}
/// Map a decoded NetCDF classic header → catalog rows (one per variable). The
/// var's shape is resolved from its `dim_ids` against the file's dimension list;
/// `ni`/`nj` take the two innermost dims (the horizontal grid, by NetCDF's
/// row-major convention). No bbox — a classic header carries no coordinate bounds.
pub fn netcdf_to_rows(hdr: &NetcdfHeader, dataset: &str) -> Vec<CatalogRow> {
hdr.vars
.iter()
.map(|v| {
let shape: Vec<u64> = v
.dim_ids
.iter()
.map(|&id| {
hdr.dims
.get(id as usize)
.map(|d| d.length.max(hdr.num_records as u64))
.unwrap_or(0)
})
.collect();
// Innermost two dims = the horizontal grid (NetCDF is row-major).
let nj = shape.iter().rev().nth(1).copied();
let ni = shape.last().copied();
let num_points: u64 = if shape.is_empty() {
0
} else {
shape.iter().product()
};
let dim_names: Vec<&str> = v
.dim_ids
.iter()
.map(|&id| {
hdr.dims
.get(id as usize)
.map(|d| d.name.as_str())
.unwrap_or("?")
})
.collect();
CatalogRow {
format: "netcdf3".into(),
dataset: dataset.into(),
variable: v.name.clone(),
reference_time: None,
ni: ni.map(|x| x as i64),
nj: nj.map(|x| x as i64),
num_points: Some(num_points as i64),
lat_min: None,
lat_max: None,
lon_min: None,
lon_max: None,
detail: format!("{}[{}]", nc_type_name(v.nc_type), dim_names.join(",")),
}
})
.collect()
}
// ═══════════════════════════════════════════════════════════════════════════
// NetCDF classic — variable data-ARRAY read (numeric types)
// ═══════════════════════════════════════════════════════════════════════════
/// The byte size of one element of a NetCDF numeric type; `char`/text and unknown
/// types are an explicit error (this reader decodes numeric arrays only).
fn nc_esize(nc_type: u32) -> Result<usize> {
match nc_type {
1 => Ok(1), // byte (int8)
3 => Ok(2), // short
4 => Ok(4), // int
5 => Ok(4), // float
6 => Ok(8), // double
2 => Err(anyhow!(
"NetCDF `char` variable is text, not a numeric array"
)),
t => Err(anyhow!("NetCDF nc_type {t} is not a readable numeric type")),
}
}
/// Read one numeric element at `pos` (big-endian) as f64, per `nc_type`.
fn nc_read_elem(b: &[u8], pos: usize, nc_type: u32) -> f64 {
match nc_type {
1 => b[pos] as i8 as f64,
3 => i16::from_be_bytes([b[pos], b[pos + 1]]) as f64,
4 => i32::from_be_bytes([b[pos], b[pos + 1], b[pos + 2], b[pos + 3]]) as f64,
5 => f32::from_be_bytes([b[pos], b[pos + 1], b[pos + 2], b[pos + 3]]) as f64,
6 => f64::from_be_bytes([
b[pos],
b[pos + 1],
b[pos + 2],
b[pos + 3],
b[pos + 4],
b[pos + 5],
b[pos + 6],
b[pos + 7],
]),
_ => f64::NAN,
}
}
/// The record (unlimited) dimension's id, if any — the dim whose header length is 0.
fn record_dim_id(hdr: &NetcdfHeader) -> Option<usize> {
hdr.dims.iter().position(|d| d.length == 0)
}
/// The record stride = the sum of every record variable's `vsize` (the header's
/// authoritative per-record slab size); 0 when the file has no record variables.
fn record_size(hdr: &NetcdfHeader) -> u64 {
let rec = record_dim_id(hdr);
hdr.vars
.iter()
.filter(|v| v.dim_ids.first().map(|&id| Some(id as usize)) == Some(rec) && rec.is_some())
.map(|v| v.vsize)
.sum()
}
fn read_slab(b: &[u8], begin: usize, count: usize, nc_type: u32, esize: usize) -> Result<Vec<f64>> {
let end = begin
.checked_add(count * esize)
.ok_or_else(|| anyhow!("NetCDF slab length overflow"))?;
if end > b.len() {
return Err(anyhow!(
"NetCDF variable data runs past the file ({end} > {} bytes)",
b.len()
));
}
Ok((0..count)
.map(|i| nc_read_elem(b, begin + i * esize, nc_type))
.collect())
}
/// Read a NetCDF **classic** variable's data ARRAY (numeric types) as row-major
/// `f64`. Handles both a fixed (contiguous) variable and a record (unlimited-dim)
/// variable — for the latter, each record's slab is gathered at `begin + r·recsize`
/// across `num_records` records. `char`/text and NetCDF-4/HDF5 stay unsupported.
pub fn decode_netcdf_variable(bytes: &[u8], hdr: &NetcdfHeader, name: &str) -> Result<Vec<f64>> {
let v = hdr
.vars
.iter()
.find(|v| v.name == name)
.ok_or_else(|| anyhow!("NetCDF variable `{name}` not found"))?;
let esize = nc_esize(v.nc_type)?;
let rec = record_dim_id(hdr);
let is_record = matches!((v.dim_ids.first(), rec), (Some(&id), Some(r)) if id as usize == r);
// Resolve each dim id to a concrete length (guarding out-of-range ids).
let dim_len = |id: u32| -> Result<usize> {
hdr.dims
.get(id as usize)
.map(|d| d.length as usize)
.ok_or_else(|| anyhow!("NetCDF variable `{name}` references unknown dim {id}"))
};
if !is_record {
let mut count = 1usize;
for &id in &v.dim_ids {
count *= dim_len(id)?;
}
read_slab(bytes, v.begin as usize, count, v.nc_type, esize)
} else {
// Innermost (non-record) dims give one record's element count.
let mut inner = 1usize;
for &id in v.dim_ids.iter().skip(1) {
inner *= dim_len(id)?;
}
let recsize = record_size(hdr) as usize;
let mut out = Vec::with_capacity(hdr.num_records as usize * inner);
for r in 0..hdr.num_records as usize {
let begin = v.begin as usize + r * recsize;
out.extend(read_slab(bytes, begin, inner, v.nc_type, esize)?);
}
Ok(out)
}
}
/// Build the unpacked-VALUES `RecordBatch` for one NetCDF variable — one row per
/// element `(message=0, variable, point index, null lat/lon, value)`, reusing the
/// [`values_schema`] the GRIB value path emits (a classic header carries no
/// per-point coordinates, so lat/lon are null).
pub fn netcdf_variable_batch(bytes: &[u8], hdr: &NetcdfHeader, name: &str) -> Result<RecordBatch> {
let vals = decode_netcdf_variable(bytes, hdr, name)?;
let n = vals.len();
let cols = vec![
i64_col(&vec![Some(0i64); n]),
str_col(&vec![Some(name.to_string()); n]),
i64_col(&(0..n as i64).map(Some).collect::<Vec<_>>()),
f64_col(&vec![None; n]),
f64_col(&vec![None; n]),
f64_col(&vals.iter().map(|v| Some(*v)).collect::<Vec<_>>()),
];
batch(&values_schema(), cols)
}
// ═══════════════════════════════════════════════════════════════════════════
// Collectors — Copernicus CDS + DWD ICON
// ═══════════════════════════════════════════════════════════════════════════
/// Which weather-model open-data source a [`WeatherCollector`] targets.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum WeatherSource {
/// Copernicus Climate Data Store retrieve API (ERA5 / reanalysis; GRIB or NetCDF).
CopernicusCds,
/// DWD ICON NWP open data (GRIB2 forecast fields on opendata.dwd.de).
DwdIcon,
/// Copernicus **ERA5** climate reanalysis — the fixed CDS `retrieve/v1`
/// endpoint, value-decoded (mode=values by default). Key from env `CDSAPI_KEY`.
CopernicusEra5,
/// Copernicus **CAMS** atmosphere / air-quality reanalysis — the fixed ADS
/// (Atmosphere Data Store) `retrieve/v1` endpoint, value-decoded. Key from env
/// `ADSAPI_KEY`.
CopernicusCams,
}
/// The NetCDF/GRIB weather-model collector. `cfg.base_url` is the file/retrieval
/// URL (the CDS download location or a DWD ICON GRIB2 file) supplied at collect
/// time; the collector fetches the bytes and decodes the container HEADER into a
/// metadata-catalog batch (the array unpack is deferred — see the module docs).
pub struct WeatherCollector {
source: WeatherSource,
}
impl WeatherCollector {
/// **Copernicus CDS** — new retrieve endpoint
/// `https://cds.climate.copernicus.eu/api/retrieve/v1/`. The API key is read
/// from env `CDSAPI_KEY` at collect time (never hardcoded), sent as
/// `Authorization: <key>` via [`AuthConfig::HeaderKey`] wired by [`auth_config`].
///
/// [`auth_config`]: WeatherCollector::auth_config
pub fn copernicus_cds() -> Self {
Self {
source: WeatherSource::CopernicusCds,
}
}
/// **DWD ICON** open data — GRIB2 forecast fields on `opendata.dwd.de`, no auth.
pub fn dwd_icon() -> Self {
Self {
source: WeatherSource::DwdIcon,
}
}
/// **Copernicus ERA5** climate reanalysis — the fixed CDS `retrieve/v1`
/// endpoint feeding the GRIB/NetCDF value decode. Defaults to `mode=values`
/// (per-point grid rows). Key from env `CDSAPI_KEY`.
pub fn era5() -> Self {
Self {
source: WeatherSource::CopernicusEra5,
}
}
/// **Copernicus CAMS** atmosphere / air-quality reanalysis — the fixed ADS
/// `retrieve/v1` endpoint feeding the GRIB/NetCDF value decode. Defaults to
/// `mode=values`. Key from env `ADSAPI_KEY`.
pub fn cams() -> Self {
Self {
source: WeatherSource::CopernicusCams,
}
}
fn name(&self) -> &'static str {
match self.source {
WeatherSource::CopernicusCds => "copernicus-cds",
WeatherSource::DwdIcon => "dwd-icon",
WeatherSource::CopernicusEra5 => "copernicus-era5",
WeatherSource::CopernicusCams => "copernicus-cams",
}
}
/// The representative endpoint (discovery + the sample seam).
pub fn endpoint(&self) -> &'static str {
match self.source {
WeatherSource::CopernicusCds => "https://cds.climate.copernicus.eu/api/retrieve/v1/",
WeatherSource::DwdIcon => "https://opendata.dwd.de/weather/nwp/icon/grib/",
WeatherSource::CopernicusEra5 => "https://cds.climate.copernicus.eu/api/retrieve/v1/",
WeatherSource::CopernicusCams => {
"https://ads.atmosphere.copernicus.eu/api/retrieve/v1/"
}
}
}
/// The env-driven [`AuthConfig`] for this source — CDS/ERA5 read `CDSAPI_KEY`,
/// CAMS reads `ADSAPI_KEY` (never hardcoded); DWD ICON open data is
/// unauthenticated. A collect wires this into `cfg.auth`.
pub fn auth_config(&self) -> AuthConfig {
match self.source {
WeatherSource::CopernicusCds | WeatherSource::CopernicusEra5 => AuthConfig::HeaderKey {
header: "Authorization".into(),
env: "CDSAPI_KEY".into(),
},
WeatherSource::CopernicusCams => AuthConfig::HeaderKey {
header: "Authorization".into(),
env: "ADSAPI_KEY".into(),
},
WeatherSource::DwdIcon => AuthConfig::None,
}
}
/// Whether this collect should unpack grid VALUES (per-point rows) rather than
/// the metadata catalog. ERA5 / CAMS are value-first reanalysis retrievals so
/// they default to values; CDS / DWD default to the catalog. Either way an
/// explicit `params.mode == "values"` forces the value unpack.
fn wants_values(&self, cfg: &CollectorCfg) -> bool {
matches!(
self.source,
WeatherSource::CopernicusEra5 | WeatherSource::CopernicusCams
) || cfg.params.get("mode").and_then(|v| v.as_str()) == Some("values")
}
}
/// **The quirks the pure-Rust GRIB2/NetCDF decoder imposes on every weather source.**
/// All four services in this module — Copernicus CDS, DWD ICON, ERA5 and CAMS — go
/// through the same [`parse_weather_bytes`] / [`parse_weather_values`] pair, so these
/// belong to the DECODER, not to any one authority. Written once and appended to each of
/// the four cards rather than twinned (LAW 5).
fn with_container_quirks(doc: crate::human_doc::HumanDoc) -> crate::human_doc::HumanDoc {
doc.with_quirk(
"This decoder reads **GRIB2 and NetCDF-CLASSIC only**. A NetCDF-4 / HDF5 container \
(the `\\x89HDF` magic) is a LOUD, named error, never a silent empty result — where \
a dataset lets you choose, ask for GRIB (or NetCDF-3). That choice is made in the \
retrieval request, not here.",
)
.with_quirk(
"`ni`, `nj` and the lat/lon bbox are only recovered for GRIB **grid definition \
template 3.0** (a regular lat/lon grid). Any other grid lands those columns NULL \
with the template number in `detail`, and the per-point value rows then carry no \
coordinates at all — the decoder refuses to assume a geometry it cannot read.",
)
.with_quirk(
"Section 6 bitmap: a 0-bit is a MISSING grid point and lands as NULL, never as a \
fabricated 0.0. A weather field full of legitimate zeroes and a field full of \
holes are not the same table, and this is the difference.",
)
.with_quirk(
"Only data-representation templates 5.0 (simple), 5.2 (complex) and 5.3 (complex + \
spatial differencing) unpack. Anything else is a named error naming the template, \
not a wrong number — if a retrieval fails here, the template number in the message \
is what to quote when you ask the provider for another packing.",
)
}
/// **The half of the card the three Copernicus data-store services share** — one
/// programme, one API family, one licence mechanism, three different portals and
/// catalogues. `store` is the portal an operator logs in to, `key_env` the environment
/// variable this collector reads its personal access token from.
///
/// Verified 2026-09-04 against ECMWF's own *CDSAPI setup* page (`…/how-to-api`) and the
/// CDS catalogue's licence vocabulary. The per-dataset licence line is quoted from that
/// page because it is the single thing that most often stops a correctly-keyed retrieval.
fn copernicus_quirks(
doc: crate::human_doc::HumanDoc,
store: &str,
key_env: &str,
) -> crate::human_doc::HumanDoc {
doc.with_quirk(
"**A registered account is not enough: you must accept the licence of EACH \
DATASET, by hand, on that dataset's page.** ECMWF states it plainly — \"One must \
agree to the Terms of Use of a dataset before downloading any data out of it. \
This step must be done manually from the dataset page (at the bottom of the \
download form).\" Until you have, a perfectly valid, correctly-keyed retrieve \
request is refused, and no amount of re-checking the key will fix it. Some \
datasets need a SECOND licence on top (the \"Additional licence to use non \
European contributions\"), so accepting one is not accepting all.",
)
.with_quirk(format!(
"korp reads the personal access token from env `{key_env}` and sends it as \
`Authorization: <key>`. The `cdsapi` Python client reads the SAME token from \
`$HOME/.cdsapirc` — copy it from {store}/how-to-api once logged in, and do not \
expect a separate \"API key\" page."
))
.with_quirk(
"The store's retrieve API is JOB-BASED: a request creates a job and the bytes come \
from a download location the job hands back. This collector fetches BYTES from \
`cfg.base_url`, so point it at that resolved download location — the \
`…/retrieve/v1/` endpoint the collector declares is a discovery label, not \
something to GET.",
)
}
#[async_trait::async_trait]
impl Collector for WeatherCollector {
/// **Four services, four cards — never one.** They share a struct and this module's
/// decoder, but Copernicus CDS, DWD ICON, ERA5 and CAMS are four different bodies with
/// four different registrations, four different licences and (for CAMS) a different
/// portal and a different key. Giving them one card would tell an operator with an
/// `ADSAPI_KEY` problem to go and read the CDS pages.
///
/// Verified 2026-09-04: ECMWF's *CDSAPI setup* page, the CDS/ADS Terms-of-use pages,
/// the CDS catalogue licence vocabulary, and — for DWD — `opendata.dwd.de/README.txt`
/// plus a live listing of the ICON GRIB tree.
fn human_doc(&self) -> Option<crate::human_doc::HumanDoc> {
use crate::human_doc::{Contact, HumanDoc, RateLimit, Registration};
// The ECMWF Support Portal link the CDS/ADS Help pages themselves publish.
const ECMWF_SUPPORT: &str = "https://jira.ecmwf.int/plugins/servlet/desk/portal/1/create/202";
let doc = match self.source {
WeatherSource::CopernicusCds => copernicus_quirks(
HumanDoc::new(
"ECMWF — which operates the Climate Data Store for the EU's Copernicus \
Climate Change Service (C3S)",
"the European Union's Copernicus programme, and the individual dataset \
producers (ECMWF and third-party contributors) whose products the \
catalogue carries — each with its own licence",
"https://cds.climate.copernicus.eu/disclaimer-privacy",
"https://cds.climate.copernicus.eu/how-to-api",
)
.with_contact(Contact::Form(ECMWF_SUPPORT.into()))
.with_registration(Registration::SelfService {
url: "https://cds.climate.copernicus.eu/".into(),
})
.with_rate_limit(RateLimit::Unknown)
.with_quirk(
"This collector defaults to the metadata CATALOG — one row per GRIB \
message or NetCDF variable, so a warehouse can see what a retrieval \
contains before paying to unpack it. Pass `params.mode = \"values\"` to \
unpack per-point grid rows instead; `copernicus-era5` and \
`copernicus-cams` do that by default.",
),
"https://cds.climate.copernicus.eu",
"CDSAPI_KEY",
),
WeatherSource::DwdIcon => HumanDoc::new(
"Deutscher Wetterdienst (DWD) — the German national meteorological service, \
on its Open Data server",
"Deutscher Wetterdienst — ICON is DWD's own numerical weather prediction \
model, so the service that serves it is also the body that made it",
"https://www.dwd.de/EN/service/legal_notice/legal_notice_node.html",
"https://www.dwd.de/EN/ourservices/opendata/opendata.html",
)
.with_contact(Contact::Form(
"https://www.dwd.de/EN/service/contact/contact_node.html".into(),
))
.with_registration(Registration::NotNeeded)
.with_rate_limit(RateLimit::NonePublished)
.with_quirk(
"**ICON files on the open-data server are bzip2-compressed.** Measured \
2026-09-04, `…/icon/grib/00/t_2m/` lists \
`icon_global_icosahedral_single-level_2026090400_000_T_2M.grib2.bz2`. This \
decoder dispatches on the container's MAGIC BYTES and rejects `.bz2` as an \
unrecognised container — decompress before handing it the bytes.",
)
.with_quirk(
"The global ICON grid is ICOSAHEDRAL, not a regular lat/lon grid, so its \
GRIB Section 3 is not template 3.0: expect `ni`, `nj` and the bbox NULL, \
with the real template number in `detail`. That is the decoder being \
honest, not a decode failure.",
)
.with_quirk(
"Filenames carry the run's date-hour (`…_2026090400_000_T_2M.grib2.bz2`) and \
the tree is keyed by run hour (`00/`, `06/`, `12/`, `18/`) then variable. A \
file URL is only valid for that run — resolve the listing at collect time \
rather than hard-coding a path.",
)
.with_quirk(
"DWD's own README states access is granted WITHOUT registration, and that \
DWD stores the requesting IP address for up to 7 days to secure and \
optimise the server. Worth knowing before you point a fleet at it — there \
is no key to rate-limit you, only DWD's patience.",
),
WeatherSource::CopernicusEra5 => copernicus_quirks(
HumanDoc::new(
"ECMWF — which operates the Climate Data Store for the EU's Copernicus \
Climate Change Service (C3S); ERA5 is ECMWF's own reanalysis",
"the European Union's Copernicus programme (C3S); ERA5 itself is produced \
by ECMWF",
"https://cds.climate.copernicus.eu/disclaimer-privacy",
"https://cds.climate.copernicus.eu/how-to-api",
)
.with_contact(Contact::Form(ECMWF_SUPPORT.into()))
.with_registration(Registration::SelfService {
url: "https://cds.climate.copernicus.eu/".into(),
})
.with_rate_limit(RateLimit::Unknown)
.with_quirk(
"This collector defaults to VALUES — per-point grid rows out of the GRIB \
Section 5/6/7 unpack or the NetCDF variable array. That is one row per \
grid point per message, so a full ERA5 field is enormous: constrain the \
retrieval (area, variable, time) at the CDS, not afterwards. The \
metadata catalog is what `copernicus-cds` lands.",
),
"https://cds.climate.copernicus.eu",
"CDSAPI_KEY",
),
WeatherSource::CopernicusCams => copernicus_quirks(
HumanDoc::new(
"ECMWF — which operates the Atmosphere Data Store for the Copernicus \
Atmosphere Monitoring Service (CAMS)",
"the European Union's Copernicus programme (CAMS); the reanalysis itself \
is produced by ECMWF",
"https://ads.atmosphere.copernicus.eu/disclaimer-privacy",
"https://ads.atmosphere.copernicus.eu/how-to-api",
)
.with_contact(Contact::Form(ECMWF_SUPPORT.into()))
.with_registration(Registration::SelfService {
url: "https://ads.atmosphere.copernicus.eu/".into(),
})
.with_rate_limit(RateLimit::Unknown)
.with_quirk(
"CAMS is served by the **Atmosphere** Data Store \
(`ads.atmosphere.copernicus.eu`), a SEPARATE portal from the Climate \
Data Store with its own catalogue and its own per-dataset licences to \
accept. korp reads its key from `ADSAPI_KEY`, not `CDSAPI_KEY` — a \
working CDS setup does not carry over by itself.",
)
.with_quirk(
"Like ERA5, this collector defaults to VALUES (per-point grid rows), not \
the metadata catalog.",
),
"https://ads.atmosphere.copernicus.eu",
"ADSAPI_KEY",
),
};
Some(with_container_quirks(doc))
}
fn meta(&self) -> CollectorMeta {
let (license, auth, description) = match self.source {
WeatherSource::CopernicusCds => (
"Copernicus/ECMWF licence (free, registration)".to_string(),
AuthKind::ApiKey,
"Copernicus CDS retrieve API (ERA5/reanalysis GRIB & NetCDF) → variable/time/bbox \
metadata catalog (header decode; array unpack deferred). Key from env CDSAPI_KEY."
.to_string(),
),
WeatherSource::DwdIcon => (
"GeoNutzV / CC-BY (DWD open data)".to_string(),
AuthKind::None,
"DWD ICON NWP open data (GRIB2 forecast fields) → variable/time/bbox metadata \
catalog (GRIB2 header decode; array unpack deferred)."
.to_string(),
),
WeatherSource::CopernicusEra5 => (
"Copernicus/ECMWF licence (free, registration)".to_string(),
AuthKind::ApiKey,
"Copernicus ERA5 climate reanalysis (CDS retrieve/v1 GRIB & NetCDF) → per-point \
grid VALUE rows (Section 5/6/7 unpack / NetCDF variable arrays). Key from env \
CDSAPI_KEY."
.to_string(),
),
WeatherSource::CopernicusCams => (
"Copernicus/ECMWF licence (free, registration)".to_string(),
AuthKind::ApiKey,
"Copernicus CAMS atmosphere / air-quality reanalysis (ADS retrieve/v1 GRIB & \
NetCDF) → per-point grid VALUE rows. Key from env ADSAPI_KEY."
.to_string(),
),
};
CollectorMeta {
name: self.name().into(),
direction: Direction::Source,
source_kind: SourceKind::BulkFile,
path: CollectorPath::RustNative,
formats: vec!["grib2".into(), "netcdf".into()],
license,
gdpr: GdprClass::Open,
auth,
join_keys: vec!["variable".into(), "reference_time".into()],
description,
}
}
fn source_target(&self) -> SourceTarget {
let (source, table) = match self.source {
WeatherSource::CopernicusCds => {
(SourceEndpoint::CopernicusCds, "copernicus_cds_catalog")
}
WeatherSource::DwdIcon => (SourceEndpoint::DwdIcon, "dwd_icon_catalog"),
WeatherSource::CopernicusEra5 => {
(SourceEndpoint::CopernicusEra5, "copernicus_era5_values")
}
WeatherSource::CopernicusCams => {
(SourceEndpoint::CopernicusCams, "copernicus_cams_values")
}
};
SourceTarget::warehouse(source, self.endpoint(), table)
}
/// **The identity differs between the two shapes this collector lands.**
///
/// A CATALOG row is one request's summary — `(dataset, variable, reference_time)`.
/// A VALUES row is one grid point of one message, so the point index is part of the
/// identity and the derived `["variable", "reference_time"]` (which the values
/// schema does not even carry a `reference_time` for) would have collapsed a whole
/// grid into one row.
fn business_key_for(&self, cfg: &CollectorCfg) -> crate::history::BusinessKey {
if self.wants_values(cfg) {
crate::history::BusinessKey::of(
&["message", "variable", "point"],
"one grid POINT of one GRIB message — the point index is part of the \
identity, or a whole grid collapses into a single row",
)
} else {
crate::history::BusinessKey::of(
&["dataset", "variable", "reference_time"],
"one catalog row per (dataset, variable, reference time)",
)
}
}
fn schema(&self, cfg: &CollectorCfg) -> Result<Schema> {
// ERA5 / CAMS (or an explicit `params.mode == "values"`) unpack the grid
// array (values_schema); the catalog sources land the metadata catalog.
if self.wants_values(cfg) {
Ok(values_schema())
} else {
Ok(catalog_schema())
}
}
async fn collect(
&self,
cfg: &CollectorCfg,
_bounds: Bounds,
sink: &BatchSink,
) -> Result<CollectStats> {
let http = cfg.http(10);
// Prefer the collector's env-driven auth (CDSAPI_KEY) unless the caller
// overrode `cfg.auth`.
let auth_cfg = match &cfg.auth {
AuthConfig::None => self.auth_config(),
other => other.clone(),
};
let auth = auth_cfg.resolve_async().await?;
let bytes = http.get_bytes(&cfg.base_url, auth.as_ref()).await?;
let dataset = cfg
.params
.get("dataset")
.and_then(|v| v.as_str())
.unwrap_or(self.name())
.to_string();
// ERA5 / CAMS (or `params.mode == "values"`) unpack the actual grid VALUES
// (per-point rows, Section 5/6/7 for GRIB, variable arrays for NetCDF); the
// catalog sources default to the metadata catalog.
if self.wants_values(cfg) {
let variable = cfg.params.get("variable").and_then(|v| v.as_str());
let rb = parse_weather_values(&bytes, variable)?;
let n = rb.num_rows() as u64;
if n == 0 {
return Ok(CollectStats {
reason: Some(format!("{} container carried no grid values", self.name())),
..Default::default()
});
}
sink.push(rb)?;
return Ok(CollectStats {
rows: n,
batches: 1,
table: cfg.table.clone(),
reason: None,
});
}
let rows = parse_weather_bytes(&bytes, &dataset)?;
if rows.is_empty() {
return Ok(CollectStats {
reason: Some(format!(
"{} container carried no messages/variables",
self.name()
)),
..Default::default()
});
}
let rb = catalog_batch(&rows)?;
let n = rb.num_rows() as u64;
sink.push(rb)?;
Ok(CollectStats {
rows: n,
batches: 1,
table: cfg.table.clone(),
reason: None,
})
}
}
// ═══════════════════════════════════════════════════════════════════════════
// In-module encoders (synthetic fixtures — no committed binary blobs)
// ═══════════════════════════════════════════════════════════════════════════
/// Spec for [`build_grib2_message`] — the fields the header decoder recovers.
#[derive(Debug, Clone)]
pub struct GribSpec {
pub discipline: u8,
pub parameter_category: u8,
pub parameter_number: u8,
pub year: u16,
pub month: u8,
pub day: u8,
pub hour: u8,
pub ni: u32,
pub nj: u32,
/// Bbox micro-degree corners (La1/Lo1/La2/Lo2), degrees.
pub lat_first: f64,
pub lon_first: f64,
pub lat_last: f64,
pub lon_last: f64,
}
fn grib_encode_i32(deg: f64) -> [u8; 4] {
let micro = (deg.abs() * 1e6).round() as u32 & 0x7FFF_FFFF;
let raw = if deg < 0.0 {
micro | 0x8000_0000
} else {
micro
};
raw.to_be_bytes()
}
/// Encode a minimal but VALID single-message GRIB2 (Sections 0,1,3-template-3.0,
/// 4-template-4.0, a stub 5 + 7, and the `7777` end) carrying `spec`'s header
/// metadata. Used to drive the decoder offline (round-trip, no committed blob).
pub fn build_grib2_message(spec: &GribSpec) -> Vec<u8> {
// Section 1 (Identification) — length 21.
let mut s1 = vec![0u8; 21];
s1[0..4].copy_from_slice(&21u32.to_be_bytes());
s1[4] = 1;
s1[9] = 2; // master tables version
s1[12..14].copy_from_slice(&spec.year.to_be_bytes());
s1[14] = spec.month;
s1[15] = spec.day;
s1[16] = spec.hour;
// Section 3 (Grid definition, template 3.0) — length 72.
let mut s3 = vec![0u8; 72];
s3[0..4].copy_from_slice(&72u32.to_be_bytes());
s3[4] = 3;
s3[6..10].copy_from_slice(&(spec.ni * spec.nj).to_be_bytes()); // number of data points
s3[12..14].copy_from_slice(&0u16.to_be_bytes()); // template 3.0
s3[14] = 6; // shape of earth
s3[30..34].copy_from_slice(&spec.ni.to_be_bytes());
s3[34..38].copy_from_slice(&spec.nj.to_be_bytes());
s3[46..50].copy_from_slice(&grib_encode_i32(spec.lat_first));
s3[50..54].copy_from_slice(&((spec.lon_first * 1e6).round() as u32).to_be_bytes());
s3[54] = 0x30; // resolution/component flags
s3[55..59].copy_from_slice(&grib_encode_i32(spec.lat_last));
s3[59..63].copy_from_slice(&((spec.lon_last * 1e6).round() as u32).to_be_bytes());
s3[71] = 0x40; // scanning mode
// Section 4 (Product definition, template 4.0) — length 34.
let mut s4 = vec![0u8; 34];
s4[0..4].copy_from_slice(&34u32.to_be_bytes());
s4[4] = 4;
s4[7..9].copy_from_slice(&0u16.to_be_bytes()); // template 4.0
s4[9] = spec.parameter_category;
s4[10] = spec.parameter_number;
// Stub Section 5 (data representation, length 21) + Section 7 (data, length 6)
// — the walker must SKIP these by length to reach the end marker.
let mut s5 = vec![0u8; 21];
s5[0..4].copy_from_slice(&21u32.to_be_bytes());
s5[4] = 5;
let mut s7 = vec![0u8; 6];
s7[0..4].copy_from_slice(&6u32.to_be_bytes());
s7[4] = 7;
let body_len = s1.len() + s3.len() + s4.len() + s5.len() + s7.len();
let total = 16 + body_len + 4; // + Section 0 + `7777`
let mut msg = Vec::with_capacity(total);
msg.extend_from_slice(b"GRIB");
msg.extend_from_slice(&[0, 0]); // reserved
msg.push(spec.discipline);
msg.push(2); // edition
msg.extend_from_slice(&(total as u64).to_be_bytes());
msg.extend_from_slice(&s1);
msg.extend_from_slice(&s3);
msg.extend_from_slice(&s4);
msg.extend_from_slice(&s5);
msg.extend_from_slice(&s7);
msg.extend_from_slice(b"7777");
debug_assert_eq!(msg.len(), total);
msg
}
/// Bits needed to represent `v` (0 → 0, a constant field).
fn bits_needed(mut v: u64) -> u8 {
let mut n = 0u8;
while v > 0 {
n += 1;
v >>= 1;
}
n
}
/// Append `val` as `nbytes` big-endian octets.
fn push_be(out: &mut Vec<u8>, val: u64, nbytes: usize) {
let b = val.to_be_bytes();
out.extend_from_slice(&b[8 - nbytes..]);
}
/// Append `val` as `nbytes` **sign-magnitude** big-endian octets (GRIB's signed form).
fn push_sign_mag(out: &mut Vec<u8>, val: i64, nbytes: usize) {
let neg = val < 0;
let mag = val.unsigned_abs().to_be_bytes();
let mut slice = mag[8 - nbytes..].to_vec();
if neg {
slice[0] |= 0x80;
}
out.extend_from_slice(&slice);
}
/// Build the Section 1 / 3 (template 3.0) / 4 (template 4.0) header trio from a
/// [`GribSpec`], with `ndp` grid points — shared by the value-carrying encoders.
fn grib_header_sections(spec: &GribSpec, ndp: u32) -> (Vec<u8>, Vec<u8>, Vec<u8>) {
let mut s1 = vec![0u8; 21];
s1[0..4].copy_from_slice(&21u32.to_be_bytes());
s1[4] = 1;
s1[9] = 2;
s1[12..14].copy_from_slice(&spec.year.to_be_bytes());
s1[14] = spec.month;
s1[15] = spec.day;
s1[16] = spec.hour;
let mut s3 = vec![0u8; 72];
s3[0..4].copy_from_slice(&72u32.to_be_bytes());
s3[4] = 3;
s3[6..10].copy_from_slice(&ndp.to_be_bytes());
s3[12..14].copy_from_slice(&0u16.to_be_bytes());
s3[14] = 6;
s3[30..34].copy_from_slice(&spec.ni.to_be_bytes());
s3[34..38].copy_from_slice(&spec.nj.to_be_bytes());
s3[46..50].copy_from_slice(&grib_encode_i32(spec.lat_first));
s3[50..54].copy_from_slice(&((spec.lon_first * 1e6).round() as u32).to_be_bytes());
s3[54] = 0x30;
s3[55..59].copy_from_slice(&grib_encode_i32(spec.lat_last));
s3[59..63].copy_from_slice(&((spec.lon_last * 1e6).round() as u32).to_be_bytes());
s3[71] = 0x40;
let mut s4 = vec![0u8; 34];
s4[0..4].copy_from_slice(&34u32.to_be_bytes());
s4[4] = 4;
s4[7..9].copy_from_slice(&0u16.to_be_bytes());
s4[9] = spec.parameter_category;
s4[10] = spec.parameter_number;
(s1, s3, s4)
}
/// Assemble a full GRIB2 message: Section 0 + the header trio + the caller's
/// Sections 5/6/7 + the `7777` end marker.
fn grib_assemble(spec: &GribSpec, ndp: u32, s5: &[u8], s6: &[u8], s7: &[u8]) -> Vec<u8> {
let (s1, s3, s4) = grib_header_sections(spec, ndp);
let body = s1.len() + s3.len() + s4.len() + s5.len() + s6.len() + s7.len();
let total = 16 + body + 4;
let mut msg = Vec::with_capacity(total);
msg.extend_from_slice(b"GRIB");
msg.extend_from_slice(&[0, 0]);
msg.push(spec.discipline);
msg.push(2);
msg.extend_from_slice(&(total as u64).to_be_bytes());
msg.extend_from_slice(&s1);
msg.extend_from_slice(&s3);
msg.extend_from_slice(&s4);
msg.extend_from_slice(s5);
msg.extend_from_slice(s6);
msg.extend_from_slice(s7);
msg.extend_from_slice(b"7777");
msg
}
/// Encode a **simple-packed** (template 5.0) GRIB2 message carrying integer
/// `values` (one per grid point, `R = 0`, `E = 0`, `D = 0` so the decode returns
/// them exactly), with a Section-6 bitmap when any point is absent (`present`
/// false). A REAL packed byte stream — the [`decode_grib2_values`] round-trip.
pub fn build_grib2_simple(spec: &GribSpec, values: &[i64], present: &[bool]) -> Vec<u8> {
let ndp = (spec.ni * spec.nj) as usize;
assert_eq!(values.len(), ndp, "values must cover every grid point");
assert_eq!(
present.len(),
ndp,
"present mask must cover every grid point"
);
let pres_vals: Vec<i64> = values
.iter()
.zip(present)
.filter(|&(_, &p)| p)
.map(|(v, _)| *v)
.collect();
let maxv = pres_vals.iter().copied().max().unwrap_or(0).max(0);
let nbits = bits_needed(maxv as u64);
let mut s5 = vec![0u8; 21];
s5[0..4].copy_from_slice(&21u32.to_be_bytes());
s5[4] = 5;
s5[5..9].copy_from_slice(&(pres_vals.len() as u32).to_be_bytes());
s5[9..11].copy_from_slice(&0u16.to_be_bytes()); // template 5.0
s5[19] = nbits;
let s6 = if present.iter().all(|&p| p) {
let mut s = vec![0u8; 6];
s[0..4].copy_from_slice(&6u32.to_be_bytes());
s[4] = 6;
s[5] = 255; // no bitmap
s
} else {
let mut bw = BitWriter::new();
for &p in present {
bw.write_bits(p as u64, 1);
}
let bm = bw.finish();
let len = 6 + bm.len();
let mut s = vec![0u8; len];
s[0..4].copy_from_slice(&(len as u32).to_be_bytes());
s[4] = 6;
s[5] = 0; // bitmap present
s[6..].copy_from_slice(&bm);
s
};
let mut bw = BitWriter::new();
for v in &pres_vals {
if nbits > 0 {
bw.write_bits(*v as u64, nbits as u32);
}
}
let payload = bw.finish();
let mut s7 = vec![0u8; 5 + payload.len()];
s7[0..4].copy_from_slice(&((5 + payload.len()) as u32).to_be_bytes());
s7[4] = 7;
s7[5..].copy_from_slice(&payload);
grib_assemble(spec, ndp as u32, &s5, &s6, &s7)
}
/// Encode a **complex-packed** GRIB2 message (template 5.2, or 5.3 when
/// `order > 0` for spatial differencing) carrying integer `values`, split into
/// `ngroups` contiguous groups (`R = 0`, `E = 0`, `D = 0`; no bitmap). Exercises
/// the group reference/width/length machinery + the difference reconstruction —
/// the [`decode_grib2_values`] round-trip for the complex paths.
pub fn build_grib2_complex(spec: &GribSpec, values: &[i64], ngroups: usize, order: u8) -> Vec<u8> {
let ndp = (spec.ni * spec.nj) as usize;
assert_eq!(values.len(), ndp, "values must cover every grid point");
assert!(ngroups >= 1 && ngroups <= ndp, "ngroups in 1..=ndp");
let n = values.len();
let nbytes = 4usize;
// Build the packed sequence P (differences for template 5.3).
let mut ival: Vec<i64> = Vec::new();
let mut gmin = 0i64;
let mut p = vec![0i64; n];
match order {
0 => p = values.to_vec(),
1 => {
ival.push(values[0]);
let mut diff = vec![0i64; n];
for i in 1..n {
diff[i] = values[i] - values[i - 1];
}
gmin = diff[1..].iter().copied().min().unwrap_or(0);
for i in 1..n {
p[i] = diff[i] - gmin;
}
}
2 => {
ival.push(values[0]);
if n > 1 {
ival.push(values[1]);
}
let mut diff = vec![0i64; n];
for i in 2..n {
diff[i] = values[i] - 2 * values[i - 1] + values[i - 2];
}
gmin = if n > 2 {
diff[2..].iter().copied().min().unwrap_or(0)
} else {
0
};
for i in 2..n {
p[i] = diff[i] - gmin;
}
}
o => panic!("spatial-diff order {o} unsupported (0, 1 or 2)"),
}
// Split P into ngroups contiguous groups; per-group ref/width.
let (base, rem) = (n / ngroups, n % ngroups);
let mut lengths = Vec::with_capacity(ngroups);
for g in 0..ngroups {
lengths.push(base + if g < rem { 1 } else { 0 });
}
let (mut idx, mut grefs, mut widths, mut groups) = (0usize, Vec::new(), Vec::new(), Vec::new());
for &l in &lengths {
let slice = &p[idx..idx + l];
idx += l;
let gr = slice.iter().copied().min().unwrap_or(0);
let mx = slice.iter().copied().max().unwrap_or(0);
widths.push(bits_needed((mx - gr) as u64));
grefs.push(gr);
groups.push(slice.to_vec());
}
let nbits_ref = bits_needed(grefs.iter().copied().max().unwrap_or(0).max(0) as u64);
let bits_gw = bits_needed(*widths.iter().max().unwrap_or(&0) as u64);
let bits_gl = bits_needed(*lengths.iter().max().unwrap_or(&0) as u64);
let template: u16 = if order == 0 { 2 } else { 3 };
let s5len = if order == 0 { 47 } else { 49 };
let mut s5 = vec![0u8; s5len];
s5[0..4].copy_from_slice(&(s5len as u32).to_be_bytes());
s5[4] = 5;
s5[5..9].copy_from_slice(&(n as u32).to_be_bytes());
s5[9..11].copy_from_slice(&template.to_be_bytes());
s5[19] = nbits_ref;
s5[21] = 1; // group splitting method: general
s5[31..35].copy_from_slice(&(ngroups as u32).to_be_bytes());
s5[36] = bits_gw;
s5[41] = 1; // length increment
s5[42..46].copy_from_slice(&(*lengths.last().unwrap() as u32).to_be_bytes());
s5[46] = bits_gl;
if order != 0 {
s5[47] = order;
s5[48] = nbytes as u8;
}
// No bitmap.
let mut s6 = vec![0u8; 6];
s6[0..4].copy_from_slice(&6u32.to_be_bytes());
s6[4] = 6;
s6[5] = 255;
// Section 7: extra descriptors (5.3) then the group bitstream.
let mut payload = Vec::new();
if order != 0 {
for iv in &ival {
push_be(&mut payload, *iv as u64, nbytes);
}
push_sign_mag(&mut payload, gmin, nbytes);
}
let mut bw = BitWriter::new();
for gr in &grefs {
bw.write_bits(*gr as u64, nbits_ref as u32);
}
for w in &widths {
bw.write_bits(*w as u64, bits_gw as u32);
}
for l in &lengths {
bw.write_bits(*l as u64, bits_gl as u32);
}
for (gi, group) in groups.iter().enumerate() {
let (w, gr) = (widths[gi], grefs[gi]);
for v in group {
if w > 0 {
bw.write_bits((*v - gr) as u64, w as u32);
}
}
}
payload.extend_from_slice(&bw.finish());
let mut s7 = vec![0u8; 5 + payload.len()];
s7[0..4].copy_from_slice(&((5 + payload.len()) as u32).to_be_bytes());
s7[4] = 7;
s7[5..].copy_from_slice(&payload);
grib_assemble(spec, ndp as u32, &s5, &s6, &s7)
}
/// A netCDF name-string field: `nelems` (u32) + bytes padded to 4.
fn nc_name_bytes(name: &str) -> Vec<u8> {
let mut out = Vec::new();
out.extend_from_slice(&(name.len() as u32).to_be_bytes());
out.extend_from_slice(name.as_bytes());
while out.len() % 4 != 0 {
out.push(0);
}
out
}
/// Encode a minimal NetCDF **classic** (CDF-1) header: the given dims + vars,
/// empty attribute lists. Used to drive the header decoder offline.
pub fn build_netcdf_classic_header(
dims: &[(&str, u32)],
vars: &[(&str, u32, Vec<u32>)],
) -> Vec<u8> {
let mut b = Vec::new();
b.extend_from_slice(b"CDF");
b.push(1); // version 1
b.extend_from_slice(&0u32.to_be_bytes()); // numrecs
// dim_list
if dims.is_empty() {
b.extend_from_slice(&0u32.to_be_bytes()); // ABSENT tag
b.extend_from_slice(&0u32.to_be_bytes());
} else {
b.extend_from_slice(&NC_DIMENSION.to_be_bytes());
b.extend_from_slice(&(dims.len() as u32).to_be_bytes());
for (name, len) in dims {
b.extend_from_slice(&nc_name_bytes(name));
b.extend_from_slice(&len.to_be_bytes());
}
}
// gatt_list = ABSENT
b.extend_from_slice(&0u32.to_be_bytes());
b.extend_from_slice(&0u32.to_be_bytes());
// var_list
if vars.is_empty() {
b.extend_from_slice(&0u32.to_be_bytes());
b.extend_from_slice(&0u32.to_be_bytes());
} else {
b.extend_from_slice(&NC_VARIABLE.to_be_bytes());
b.extend_from_slice(&(vars.len() as u32).to_be_bytes());
for (name, nc_type, dim_ids) in vars {
b.extend_from_slice(&nc_name_bytes(name));
b.extend_from_slice(&(dim_ids.len() as u32).to_be_bytes());
for id in dim_ids {
b.extend_from_slice(&id.to_be_bytes());
}
// vatt_list = ABSENT
b.extend_from_slice(&0u32.to_be_bytes());
b.extend_from_slice(&0u32.to_be_bytes());
b.extend_from_slice(&nc_type.to_be_bytes());
b.extend_from_slice(&0u32.to_be_bytes()); // vsize
b.extend_from_slice(&0u32.to_be_bytes()); // begin (CDF-1 u32)
}
}
b
}
/// Round a byte count up to the next 4-byte boundary (NetCDF `vsize` padding).
fn pad4(n: usize) -> usize {
n.div_ceil(4) * 4
}
/// Write one numeric element (big-endian) at `pos`, per `nc_type`.
fn nc_write_elem(buf: &mut [u8], pos: usize, val: f64, nc_type: u32) {
match nc_type {
1 => buf[pos] = val as i8 as u8,
3 => buf[pos..pos + 2].copy_from_slice(&(val as i16).to_be_bytes()),
4 => buf[pos..pos + 4].copy_from_slice(&(val as i32).to_be_bytes()),
5 => buf[pos..pos + 4].copy_from_slice(&(val as f32).to_be_bytes()),
6 => buf[pos..pos + 8].copy_from_slice(&val.to_be_bytes()),
_ => {}
}
}
/// Encode a NetCDF **classic** (CDF-1) file WITH variable DATA — the header (as
/// [`build_netcdf_classic_header`]) plus a laid-out data section, so
/// [`decode_netcdf_variable`] round-trips. Fixed variables are stored contiguously
/// after the header; record (unlimited-dim, header length 0) variables follow in
/// the interleaved record section (`begin + r·recsize`). A dim of length 0 is the
/// record dim; a record variable's `data` is `numrecs × product(inner dims)`.
pub fn build_netcdf_classic_with_data(
numrecs: u32,
dims: &[(&str, u32)],
vars: &[(&str, u32, Vec<u32>, Vec<f64>)],
) -> Vec<u8> {
let rec_dim = dims.iter().position(|(_, l)| *l == 0);
// Per-variable layout metadata: element size, record-ness, per-slab count, vsize.
struct VMeta {
esize: usize,
is_rec: bool,
vsize: usize,
slab: usize,
}
let meta: Vec<VMeta> = vars
.iter()
.map(|(_, nc_type, dim_ids, _)| {
let esize = nc_esize(*nc_type).expect("numeric nc_type");
let is_rec = rec_dim.is_some()
&& dim_ids
.first()
.map(|&id| Some(id as usize) == rec_dim)
.unwrap_or(false);
let mut slab = 1usize;
for (i, &id) in dim_ids.iter().enumerate() {
if is_rec && i == 0 {
continue; // the record dim is not part of one slab
}
slab *= dims[id as usize].1 as usize;
}
VMeta {
esize,
is_rec,
vsize: pad4(slab * esize),
slab,
}
})
.collect();
// Header with placeholder `begin` fields (offsets recorded for patching).
let mut b = Vec::new();
b.extend_from_slice(b"CDF");
b.push(1);
b.extend_from_slice(&numrecs.to_be_bytes());
if dims.is_empty() {
b.extend_from_slice(&0u32.to_be_bytes());
b.extend_from_slice(&0u32.to_be_bytes());
} else {
b.extend_from_slice(&NC_DIMENSION.to_be_bytes());
b.extend_from_slice(&(dims.len() as u32).to_be_bytes());
for (name, len) in dims {
b.extend_from_slice(&nc_name_bytes(name));
b.extend_from_slice(&len.to_be_bytes());
}
}
b.extend_from_slice(&0u32.to_be_bytes()); // gatt_list ABSENT
b.extend_from_slice(&0u32.to_be_bytes());
let mut begin_field_at = Vec::with_capacity(vars.len());
if vars.is_empty() {
b.extend_from_slice(&0u32.to_be_bytes());
b.extend_from_slice(&0u32.to_be_bytes());
} else {
b.extend_from_slice(&NC_VARIABLE.to_be_bytes());
b.extend_from_slice(&(vars.len() as u32).to_be_bytes());
for (i, (name, nc_type, dim_ids, _)) in vars.iter().enumerate() {
b.extend_from_slice(&nc_name_bytes(name));
b.extend_from_slice(&(dim_ids.len() as u32).to_be_bytes());
for id in dim_ids {
b.extend_from_slice(&id.to_be_bytes());
}
b.extend_from_slice(&0u32.to_be_bytes()); // vatt_list ABSENT
b.extend_from_slice(&0u32.to_be_bytes());
b.extend_from_slice(&nc_type.to_be_bytes());
b.extend_from_slice(&(meta[i].vsize as u32).to_be_bytes());
begin_field_at.push(b.len());
b.extend_from_slice(&0u32.to_be_bytes()); // begin (CDF-1 u32) placeholder
}
}
let header_len = b.len();
// Fixed vars contiguous from header end; then the interleaved record section.
let mut begins = vec![0usize; vars.len()];
let mut off = header_len;
for i in 0..vars.len() {
if !meta[i].is_rec {
begins[i] = off;
off += meta[i].vsize;
}
}
let rec_begin = off;
let recsize: usize = (0..vars.len())
.filter(|&i| meta[i].is_rec)
.map(|i| meta[i].vsize)
.sum();
let mut rec_off = rec_begin;
for i in 0..vars.len() {
if meta[i].is_rec {
begins[i] = rec_off;
rec_off += meta[i].vsize;
}
}
for i in 0..vars.len() {
let o = begin_field_at[i];
b[o..o + 4].copy_from_slice(&(begins[i] as u32).to_be_bytes());
}
let total = if recsize > 0 {
rec_begin + numrecs as usize * recsize
} else {
rec_begin
};
let mut file = b;
file.resize(total.max(header_len), 0);
for (i, (_, nc_type, _, data)) in vars.iter().enumerate() {
let esize = meta[i].esize;
if !meta[i].is_rec {
for (j, val) in data.iter().enumerate() {
nc_write_elem(&mut file, begins[i] + j * esize, *val, *nc_type);
}
} else {
let inner = meta[i].slab;
for r in 0..numrecs as usize {
for j in 0..inner {
let idx = r * inner + j;
if idx < data.len() {
nc_write_elem(
&mut file,
begins[i] + r * recsize + j * esize,
data[idx],
*nc_type,
);
}
}
}
}
}
file
}
#[cfg(test)]
mod tests {
use super::*;
use arrow_array::{Array, Float64Array, Int64Array, StringArray};
fn sample_spec() -> GribSpec {
GribSpec {
discipline: 0,
parameter_category: 0,
parameter_number: 0, // temperature
year: 2026,
month: 7,
day: 18,
hour: 12,
ni: 360,
nj: 181,
lat_first: 90.0,
lon_first: 0.0,
lat_last: -90.0,
lon_last: 359.0,
}
}
/// A hand-built GRIB2 message round-trips through the section-header decoder:
/// discipline/parameter, reference time, grid Ni/Nj and the lat/lon bbox all
/// recover, and the stub Sections 5/7 are skipped to the `7777` end. RED if
/// any section offset or the sign-magnitude lat decode breaks.
#[test]
fn grib2_message_header_round_trips() {
let spec = sample_spec();
let bytes = build_grib2_message(&spec);
let msgs = parse_grib2(&bytes).expect("parse grib2");
assert_eq!(msgs.len(), 1, "one message");
let m = &msgs[0];
assert_eq!(m.discipline, 0);
assert_eq!(m.parameter_number, 0, "temperature parameter");
assert_eq!(
m.reference_time, "2026-07-18T12:00:00Z",
"Section-1 reference time"
);
assert_eq!(m.grid_template, 0, "lat/lon grid template");
assert_eq!(m.ni, 360);
assert_eq!(m.nj, 181);
assert_eq!(m.num_data_points, 360 * 181);
assert!((m.lat_first - 90.0).abs() < 1e-6, "La1");
assert!(
(m.lat_last + 90.0).abs() < 1e-6,
"La2 (negative, sign-magnitude)"
);
assert!((m.lon_last - 359.0).abs() < 1e-6, "Lo2");
}
/// The decoded GRIB message maps to a typed catalog batch (bbox min/max
/// normalised, num_points carried). The offline DATA contract.
#[test]
fn grib_catalog_batch_is_typed() {
let bytes = build_grib2_message(&sample_spec());
let msgs = parse_grib2(&bytes).unwrap();
let rows = grib_to_rows(&msgs, "reanalysis-era5-single-levels");
let rb = catalog_batch(&rows).expect("catalog batch");
assert_eq!(rb.num_rows(), 1);
let fmt = rb
.column_by_name("format")
.and_then(|a| a.as_any().downcast_ref::<StringArray>())
.unwrap();
assert_eq!(fmt.value(0), "grib2");
let latmin = rb
.column_by_name("lat_min")
.and_then(|a| a.as_any().downcast_ref::<Float64Array>())
.unwrap();
assert!(
(latmin.value(0) + 90.0).abs() < 1e-6,
"bbox min-normalised to -90"
);
let np = rb
.column_by_name("num_points")
.and_then(|a| a.as_any().downcast_ref::<Int64Array>())
.unwrap();
assert_eq!(np.value(0), 360 * 181);
crate::functional_status(
"korp-collectors/dwd-icon",
"grib_header_catalog",
true,
&format!("{} GRIB2 message → metadata catalog", rb.num_rows()),
);
}
/// Two concatenated GRIB2 messages (as a `.grib2` file holds) both decode.
#[test]
fn grib2_multi_message_file_decodes_all() {
let mut spec2 = sample_spec();
spec2.parameter_number = 1; // a second variable
let mut file = build_grib2_message(&sample_spec());
file.extend_from_slice(&build_grib2_message(&spec2));
let msgs = parse_grib2(&file).expect("parse multi");
assert_eq!(msgs.len(), 2, "both concatenated messages decoded");
assert_eq!(msgs[1].parameter_number, 1);
}
/// A NetCDF classic header round-trips: dimensions + variables (name, nc_type,
/// dim spans) decode, attribute lists are walked, and the variable shape maps
/// to the catalog rows (innermost two dims = ni/nj, product = num_points).
#[test]
fn netcdf_classic_header_round_trips() {
// dims: time=1, lat=181, lon=360 ; var t2m(time,lat,lon) float.
let dims = [("time", 1u32), ("lat", 181u32), ("lon", 360u32)];
let vars = vec![("t2m", 5u32, vec![0u32, 1, 2])]; // nc_type 5 = float
let bytes = build_netcdf_classic_header(&dims, &vars);
let hdr = parse_netcdf_classic_header(&bytes).expect("parse netcdf header");
assert_eq!(hdr.dims.len(), 3);
assert_eq!(hdr.dims[2].name, "lon");
assert_eq!(hdr.dims[2].length, 360);
assert_eq!(hdr.vars.len(), 1);
assert_eq!(hdr.vars[0].name, "t2m");
assert_eq!(hdr.vars[0].nc_type, 5);
assert_eq!(hdr.vars[0].dim_ids, vec![0, 1, 2]);
let rows = netcdf_to_rows(&hdr, "era5-netcdf");
assert_eq!(rows.len(), 1);
assert_eq!(rows[0].format, "netcdf3");
assert_eq!(rows[0].ni, Some(360), "innermost dim = lon = ni");
assert_eq!(rows[0].nj, Some(181), "next dim = lat = nj");
assert_eq!(
rows[0].num_points,
Some(181 * 360),
"time(1) × lat(181) × lon(360)"
);
assert_eq!(rows[0].variable, "t2m");
assert!(
rows[0].detail.starts_with("float["),
"detail carries nc_type + dims: {}",
rows[0].detail
);
// Route via the magic-byte dispatcher too.
let rows2 = parse_weather_bytes(&bytes, "era5-netcdf").expect("dispatch netcdf");
assert_eq!(rows2, rows, "dispatcher matches direct parse");
}
/// The magic-byte dispatcher routes GRIB vs NetCDF, and DEFERS NetCDF-4/HDF5
/// with a loud error (not a silent empty result).
#[test]
fn dispatcher_routes_and_defers_hdf5() {
let grib = build_grib2_message(&sample_spec());
let rows = parse_weather_bytes(&grib, "ds").expect("grib route");
assert_eq!(rows[0].format, "grib2");
let hdf5 = b"\x89HDF\r\n\x1a\n";
let err = parse_weather_bytes(hdf5, "ds").unwrap_err().to_string();
assert!(
err.contains("DEFERRED"),
"HDF5 array decode deferred loudly: {err}"
);
}
/// The CDS collector reads its key from env `CDSAPI_KEY` (never hardcoded);
/// DWD ICON is unauthenticated. The typed source endpoints are declared.
#[test]
fn collectors_declare_auth_and_endpoints() {
let cds = WeatherCollector::copernicus_cds();
match cds.auth_config() {
AuthConfig::HeaderKey { env, .. } => assert_eq!(env, "CDSAPI_KEY", "CDS key from env"),
other => panic!("CDS should be header-key auth, got {other:?}"),
}
assert_eq!(cds.source_target().source, SourceEndpoint::CopernicusCds);
assert!(
cds.endpoint().contains("api/retrieve/v1"),
"new CDS retrieve endpoint"
);
let dwd = WeatherCollector::dwd_icon();
assert!(
matches!(dwd.auth_config(), AuthConfig::None),
"DWD open data, no auth"
);
assert_eq!(dwd.source_target().source, SourceEndpoint::DwdIcon);
assert_eq!(
dwd.schema(&CollectorCfg::default()).unwrap(),
catalog_schema()
);
}
// ── GRIB2 Section 5/6/7 data-VALUE unpack (the array-decode contract) ────
/// A 4×3 grid spec for the value-unpack round trips.
fn grid_spec() -> GribSpec {
GribSpec {
discipline: 0,
parameter_category: 0,
parameter_number: 0,
year: 2026,
month: 7,
day: 18,
hour: 0,
ni: 4,
nj: 3,
lat_first: 60.0,
lon_first: 10.0,
lat_last: 58.0,
lon_last: 13.0,
}
}
/// Assert every present decoded value equals the source integer (as f64).
fn assert_values_eq(got: &[Option<f64>], want: &[i64]) {
assert_eq!(got.len(), want.len(), "value count");
for (i, (g, w)) in got.iter().zip(want).enumerate() {
let g = g.unwrap_or_else(|| panic!("value {i} missing"));
assert!((g - *w as f64).abs() < 1e-9, "value {i}: got {g}, want {w}");
}
}
/// SIMPLE packing (template 5.0): a full 12-point grid round-trips exactly
/// through the encode → `decode_grib2_values` path, and the stats summarise
/// the present points. RED if the bit reader/writer or the R/E/D scale breaks.
#[test]
fn grib2_simple_packing_round_trips_values() {
let spec = grid_spec();
let values: Vec<i64> = vec![0, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144];
let present = vec![true; 12];
let bytes = build_grib2_simple(&spec, &values, &present);
let decoded = decode_grib2_values(&bytes).expect("decode simple grib2");
assert_eq!(decoded.len(), 1, "one message");
assert_values_eq(&decoded[0].values, &values);
let (count, min, max, _mean) = decoded[0].stats();
assert_eq!(count, 12, "all present");
assert_eq!(min, Some(0.0));
assert_eq!(max, Some(144.0));
crate::functional_status(
"korp-collectors/dwd-icon",
"grib_simple_value_unpack",
true,
&format!("{} grid values simple-packed → f64", count),
);
}
/// SIMPLE packing with a Section-6 BITMAP: absent points decode to `None`
/// (missing), never a fabricated 0.0, and present points still recover.
#[test]
fn grib2_simple_packing_bitmap_marks_missing() {
let spec = grid_spec();
let values: Vec<i64> = vec![10, 0, 20, 0, 30, 40, 0, 50, 60, 70, 0, 80];
// Mark the four zero-holes as MISSING via the bitmap.
let present: Vec<bool> = values.iter().map(|v| *v != 0).collect();
let bytes = build_grib2_simple(&spec, &values, &present);
let decoded = decode_grib2_values(&bytes).expect("decode bitmapped grib2");
let vals = &decoded[0].values;
assert_eq!(vals.len(), 12, "full grid width preserved");
for (i, (v, p)) in vals.iter().zip(&present).enumerate() {
if *p {
assert_eq!(*v, Some(values[i] as f64), "present point {i}");
} else {
assert!(v.is_none(), "bitmap-absent point {i} is missing, not 0.0");
}
}
let (count, _, _, _) = decoded[0].stats();
assert_eq!(count, 8, "eight present points (four bitmapped out)");
}
/// COMPLEX packing (template 5.2), split across TWO groups, round-trips — the
/// per-group reference/width/length machinery decodes correctly.
#[test]
fn grib2_complex_packing_round_trips_values() {
let spec = grid_spec();
let values: Vec<i64> = vec![100, 101, 103, 106, 110, 115, 121, 128, 136, 145, 155, 166];
let bytes = build_grib2_complex(&spec, &values, 2, 0);
let decoded = decode_grib2_values(&bytes).expect("decode complex grib2");
assert_values_eq(&decoded[0].values, &values);
crate::functional_status(
"korp-collectors/dwd-icon",
"grib_complex_value_unpack",
true,
"12 grid values complex-packed (2 groups) → f64",
);
}
/// COMPLEX + first-order SPATIAL DIFFERENCING (template 5.3): the initial
/// value + overall minimum extra descriptors drive the cumulative-sum
/// reconstruction. RED if the un-differencing or the sign-magnitude minimum breaks.
#[test]
fn grib2_complex_spatial_diff_order1_round_trips() {
let spec = grid_spec();
// A field whose first differences vary in sign (exercises the signed gmin).
let values: Vec<i64> = vec![50, 52, 51, 55, 60, 58, 65, 70, 68, 72, 80, 79];
let bytes = build_grib2_complex(&spec, &values, 3, 1);
let decoded = decode_grib2_values(&bytes).expect("decode grib2 spatial-diff o1");
assert_values_eq(&decoded[0].values, &values);
}
/// COMPLEX + second-order SPATIAL DIFFERENCING (template 5.3, order 2).
#[test]
fn grib2_complex_spatial_diff_order2_round_trips() {
let spec = grid_spec();
let values: Vec<i64> = vec![0, 3, 8, 15, 24, 35, 48, 63, 80, 99, 120, 143]; // ~quadratic
let bytes = build_grib2_complex(&spec, &values, 2, 2);
let decoded = decode_grib2_values(&bytes).expect("decode grib2 spatial-diff o2");
assert_values_eq(&decoded[0].values, &values);
}
/// The unpacked values map to a typed Arrow `RecordBatch` — one row per grid
/// point, lat/lon filled for the regular lat/lon grid, missing → null value.
#[test]
fn grib_values_batch_is_typed() {
use arrow_array::{Array, Float64Array, Int64Array};
let spec = grid_spec();
let values: Vec<i64> = vec![7, 0, 9, 11, 13, 15, 0, 19, 21, 23, 25, 27];
let present: Vec<bool> = values.iter().map(|v| *v != 0).collect();
let bytes = build_grib2_simple(&spec, &values, &present);
let decoded = decode_grib2_values(&bytes).expect("decode");
let rb = grib_values_batch(&decoded).expect("values batch");
assert_eq!(rb.num_rows(), 12, "one row per grid point");
let val = rb
.column_by_name("value")
.and_then(|a| a.as_any().downcast_ref::<Float64Array>())
.unwrap();
assert!((val.value(0) - 7.0).abs() < 1e-9, "first value");
assert!(val.is_null(1), "bitmap-absent point → null value");
let pt = rb
.column_by_name("point")
.and_then(|a| a.as_any().downcast_ref::<Int64Array>())
.unwrap();
assert_eq!(pt.value(11), 11, "point index");
let lat = rb
.column_by_name("lat")
.and_then(|a| a.as_any().downcast_ref::<Float64Array>())
.unwrap();
assert!(!lat.is_null(0), "regular lat/lon grid fills lat");
assert!((lat.value(0) - 60.0).abs() < 1e-9, "point 0 sits at La1");
}
// ── NetCDF classic variable data-ARRAY read ─────────────────────────────
/// A FIXED (contiguous) NetCDF-classic variable round-trips: the header's
/// `begin`/`vsize` locate the slab and every element decodes to f64. RED if
/// the `begin` offset or the big-endian element read breaks.
#[test]
fn netcdf_fixed_variable_array_round_trips() {
// dims lat=3, lon=4 ; var t2m(lat,lon) float = 12 values.
let dims = [("lat", 3u32), ("lon", 4u32)];
let data: Vec<f64> = (0..12).map(|i| 273.15 + i as f64).collect();
let vars = vec![("t2m", 5u32, vec![0u32, 1], data.clone())];
let bytes = build_netcdf_classic_with_data(0, &dims, &vars);
let hdr = parse_netcdf_classic_header(&bytes).expect("parse header");
assert_eq!(hdr.vars[0].name, "t2m");
assert!(hdr.vars[0].vsize >= 48, "12 floats × 4 bytes"); // padded to ≥48
let got = decode_netcdf_variable(&bytes, &hdr, "t2m").expect("decode t2m");
assert_eq!(got.len(), 12, "3×4 grid");
for (i, (g, w)) in got.iter().zip(&data).enumerate() {
assert!((g - w).abs() < 1e-4, "element {i}: {g} vs {w}");
}
crate::functional_status(
"korp-collectors/copernicus-cds",
"netcdf_fixed_variable_read",
true,
&format!("{} NetCDF-classic values → f64", got.len()),
);
}
/// A RECORD (unlimited-dim) variable round-trips across records: a `double`
/// mixed with a fixed `float` coordinate variable exercises the interleaved
/// record section (`begin + r·recsize`) alongside a contiguous fixed var.
#[test]
fn netcdf_record_variable_array_round_trips() {
// dims time=UNLIMITED(0), lat=2, lon=2 ; fixed lat_coord(lat), record temp(time,lat,lon).
let dims = [("time", 0u32), ("lat", 2u32), ("lon", 2u32)];
let lat_coord = vec![58.0f64, 60.0];
// 3 records × (2×2) = 12 doubles.
let temp: Vec<f64> = (0..12).map(|i| 10.0 + i as f64 * 0.5).collect();
let vars = vec![
("lat_coord", 5u32, vec![1u32], lat_coord.clone()), // fixed float on lat
("temp", 6u32, vec![0u32, 1, 2], temp.clone()), // record double
];
let bytes = build_netcdf_classic_with_data(3, &dims, &vars);
let hdr = parse_netcdf_classic_header(&bytes).expect("parse header");
assert_eq!(hdr.num_records, 3, "three records");
let got_lat = decode_netcdf_variable(&bytes, &hdr, "lat_coord").expect("decode lat");
assert_eq!(got_lat.len(), 2);
assert!((got_lat[1] - 60.0).abs() < 1e-4, "fixed coord var read");
let got_temp = decode_netcdf_variable(&bytes, &hdr, "temp").expect("decode temp");
assert_eq!(got_temp.len(), 12, "3 records × 4 grid points");
for (i, (g, w)) in got_temp.iter().zip(&temp).enumerate() {
assert!((g - w).abs() < 1e-9, "record element {i}: {g} vs {w}");
}
}
/// A `char`/text variable is rejected (this reader decodes numeric arrays);
/// an unknown variable name is a hard error, not a silent empty vector.
#[test]
fn netcdf_non_numeric_and_missing_variable_rejected() {
// A `char` variable (nc_type 2) via the header-only builder — the numeric
// data path can't encode text, so the reader must reject it up front.
let dims = [("n", 4u32)];
let vars = vec![("label", 2u32, vec![0u32])]; // nc_type 2 = char
let bytes = build_netcdf_classic_header(&dims, &vars);
let hdr = parse_netcdf_classic_header(&bytes).expect("parse header");
assert!(
decode_netcdf_variable(&bytes, &hdr, "label").is_err(),
"char var rejected"
);
assert!(
decode_netcdf_variable(&bytes, &hdr, "nope").is_err(),
"missing var rejected"
);
}
/// The NetCDF variable values map to the same typed VALUES batch the GRIB path
/// emits (message=0, null lat/lon, one row per element).
#[test]
fn netcdf_variable_batch_is_typed() {
let dims = [("lat", 2u32), ("lon", 2u32)];
let data = vec![1.5f64, 2.5, 3.5, 4.5];
let vars = vec![("q", 6u32, vec![0u32, 1], data)];
let bytes = build_netcdf_classic_with_data(0, &dims, &vars);
let hdr = parse_netcdf_classic_header(&bytes).expect("hdr");
let rb = netcdf_variable_batch(&bytes, &hdr, "q").expect("nc values batch");
assert_eq!(rb.num_rows(), 4);
let val = rb
.column_by_name("value")
.and_then(|a| a.as_any().downcast_ref::<Float64Array>())
.unwrap();
assert!((val.value(2) - 3.5).abs() < 1e-9);
let lat = rb
.column_by_name("lat")
.and_then(|a| a.as_any().downcast_ref::<Float64Array>())
.unwrap();
assert!(
lat.is_null(0),
"classic header carries no per-point coordinates → null lat"
);
let var = rb
.column_by_name("variable")
.and_then(|a| a.as_any().downcast_ref::<StringArray>())
.unwrap();
assert_eq!(var.value(0), "q");
}
/// **ERA5 / CAMS** — the climate-reanalysis + air-quality presets declare their
/// own typed `source → target` (fixed retrieve/v1 endpoint, values table, new
/// SourceEndpoint variant), DEFAULT to the grid-VALUE unpack (schema =
/// values_schema, no `params.mode` needed), and a real simple-packed GRIB2
/// message round-trips through `parse_weather_values` to the RIGHT per-point
/// values — RED if any preset wiring or the value decode breaks.
#[test]
fn era5_and_cams_default_to_grid_value_rows() {
// Each preset: name, SourceEndpoint, values table, endpoint host, auth env.
let era5 = WeatherCollector::era5();
let st = era5.source_target();
assert_eq!(st.source, SourceEndpoint::CopernicusEra5);
assert_eq!(st.target_table, "copernicus_era5_values");
assert!(
st.source_endpoint.contains("retrieve/v1"),
"ERA5 fixed retrieve/v1 endpoint"
);
assert!(matches!(
era5.auth_config(),
AuthConfig::HeaderKey { ref env, .. } if env == "CDSAPI_KEY"
));
let cams = WeatherCollector::cams();
let st = cams.source_target();
assert_eq!(st.source, SourceEndpoint::CopernicusCams);
assert_eq!(st.target_table, "copernicus_cams_values");
assert!(
st.source_endpoint.contains("ads.atmosphere"),
"CAMS ADS retrieve/v1 endpoint"
);
assert!(matches!(
cams.auth_config(),
AuthConfig::HeaderKey { ref env, .. } if env == "ADSAPI_KEY"
));
// Default (no params.mode) → the VALUES schema, not the catalog.
let cfg = CollectorCfg::new(era5.endpoint(), "copernicus_era5_values");
let schema = era5.schema(&cfg).expect("era5 schema");
let cols: Vec<&str> = schema.fields().iter().map(|f| f.name().as_str()).collect();
assert_eq!(
cols,
vec!["message", "variable", "point", "lat", "lon", "value"],
"values schema"
);
// A real 2×2 simple-packed GRIB2 field unpacks to 4 per-point value rows.
let spec = GribSpec {
discipline: 0,
parameter_category: 0,
parameter_number: 0,
year: 2026,
month: 7,
day: 18,
hour: 0,
ni: 2,
nj: 2,
lat_first: 60.0,
lon_first: 10.0,
lat_last: 59.0,
lon_last: 11.0,
};
let grib = build_grib2_simple(&spec, &[10, 20, 30, 40], &[true; 4]);
let rb = parse_weather_values(&grib, None).expect("era5 value unpack");
assert_eq!(rb.num_rows(), 4, "2×2 grid → four per-point value rows");
let val = rb
.column_by_name("value")
.and_then(|a| a.as_any().downcast_ref::<Float64Array>())
.expect("value f64");
assert!((val.value(0) - 10.0).abs() < 1e-9);
assert!((val.value(3) - 40.0).abs() < 1e-9);
crate::functional_status(
"korp-collectors/copernicus-era5",
"reanalysis_grid_values",
true,
&format!(
"{} ERA5/CAMS per-point grid value rows from a GRIB2 retrieval",
rb.num_rows()
),
);
}
}