sidereon_core/ionex/
samples.rs

1//! Sample-backed IONEX vertical-TEC source.
2//!
3//! The canonical IONEX intermediate representation is a set of vertical-TEC
4//! maps on a strictly increasing epoch axis, with descending latitude nodes,
5//! ascending longitude nodes, signed grid steps, shell geometry, and optional
6//! RMS maps. IONEX text is one serialization of that IR; [`super::Ionex`] is
7//! the parser. This module builds the same evaluatable product directly from
8//! samples, with no text in the loop, and drives the exact same slant-delay
9//! evaluator the parsed path uses.
10//!
11//! # Byte-identical parity with the parser path
12//!
13//! [`Ionex::from_samples`] accepts the same field values the parser stores:
14//! map epochs, TEC/RMS grids in TECU, node axes, signed steps, shell geometry,
15//! and `EXPONENT`. [`Ionex::tec_grid_samples`] clones those fields out of a
16//! parsed or sample-built product. Therefore
17//! `Ionex::from_samples(ionex.tec_grid_samples())` rebuilds an equal product
18//! byte-for-byte in every stored float and epoch. Unlike SP3 samples, there is
19//! no SI reconstruction boundary here: VTEC is TECU on both sides. The only
20//! lossy boundary remains serialize-through-text, where the existing writer
21//! recovers scaled integer fields with `round(value / 10^EXPONENT)`.
22
23use super::grid::{Ionex, IonexParts};
24use super::j2000_seconds_from_instant;
25use crate::astro::time::model::{Instant, InstantRepr};
26
27const IONEX_AXIS_DEG_LIMIT: f64 = 360.0;
28const NANOS_PER_SECOND: i128 = 1_000_000_000;
29
30/// One vertical-TEC sample at one grid node.
31#[derive(Debug, Clone, Copy, PartialEq)]
32pub struct TecSample {
33    /// Map epoch.
34    pub epoch: Instant,
35    /// Latitude node in degrees.
36    pub lat_deg: f64,
37    /// Longitude node in degrees.
38    pub lon_deg: f64,
39    /// Vertical TEC in TECU.
40    pub vtec_tecu: f64,
41    /// Optional RMS value in TECU.
42    pub rms_tecu: Option<f64>,
43}
44
45/// Whole-grid IONEX vertical-TEC samples.
46#[derive(Debug, Clone, PartialEq)]
47pub struct TecGridSamples {
48    /// Map epochs as instants, strictly increasing.
49    pub map_epochs: Vec<Instant>,
50    /// Latitude node values in degrees, descending.
51    pub lat_nodes_deg: Vec<f64>,
52    /// Longitude node values in degrees, ascending.
53    pub lon_nodes_deg: Vec<f64>,
54    /// Signed latitude step in degrees.
55    pub dlat_deg: f64,
56    /// Signed longitude step in degrees.
57    pub dlon_deg: f64,
58    /// Single-layer shell height in kilometers.
59    pub shell_height_km: f64,
60    /// Mean earth radius used by the geometry, in kilometers.
61    pub base_radius_km: f64,
62    /// The IONEX `EXPONENT` header field.
63    pub exponent: i32,
64    /// Per-map vertical-TEC grids, indexed `[map][i_lat][i_lon]` (TECU).
65    pub tec_maps: Vec<Vec<Vec<f64>>>,
66    /// Per-map RMS grids, indexed `[map][i_lat][i_lon]` (TECU); empty if absent.
67    pub rms_maps: Vec<Vec<Vec<f64>>>,
68}
69
70/// Validation failure building an IONEX sample source.
71#[derive(Debug, Clone, PartialEq)]
72pub enum TecSamplesError {
73    /// No TEC samples were supplied.
74    Empty,
75    /// A latitude or longitude axis has fewer than two nodes.
76    TooFewNodes(usize),
77    /// Latitude nodes are not strictly descending.
78    NonMonotonicLat,
79    /// Longitude nodes are not strictly ascending.
80    NonMonotonicLon,
81    /// Map epochs are not strictly increasing.
82    NonMonotonicEpochs,
83    /// A map epoch cannot be expressed as an exact integer J2000 second.
84    EpochNotRepresentable,
85    /// Grid dimensions do not match the epoch or node axes.
86    ShapeMismatch,
87    /// RMS map count or node coverage does not match the TEC maps.
88    RmsCountMismatch,
89    /// A supplied float was NaN or infinite.
90    NonFiniteValue,
91    /// A signed grid step is zero or has the wrong sign for the node ordering.
92    NonPositiveStep,
93    /// An axis coordinate or step falls outside `[-360, 360]` degrees.
94    AxisOutOfRange(f64),
95}
96
97impl core::fmt::Display for TecSamplesError {
98    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
99        match self {
100            Self::Empty => write!(f, "no IONEX TEC samples supplied"),
101            Self::TooFewNodes(count) => {
102                write!(f, "IONEX grid axis has {count} nodes; need at least two")
103            }
104            Self::NonMonotonicLat => {
105                write!(f, "IONEX latitude nodes must be strictly descending")
106            }
107            Self::NonMonotonicLon => {
108                write!(f, "IONEX longitude nodes must be strictly ascending")
109            }
110            Self::NonMonotonicEpochs => {
111                write!(f, "IONEX map epochs must be strictly increasing")
112            }
113            Self::EpochNotRepresentable => {
114                write!(f, "IONEX map epoch is not an exact integer J2000 second")
115            }
116            Self::ShapeMismatch => {
117                write!(f, "IONEX TEC grid dimensions do not match the axes")
118            }
119            Self::RmsCountMismatch => write!(f, "IONEX RMS maps do not match TEC maps"),
120            Self::NonFiniteValue => write!(f, "IONEX sample value is not finite"),
121            Self::NonPositiveStep => write!(f, "IONEX grid step is not valid"),
122            Self::AxisOutOfRange(value) => {
123                write!(f, "IONEX axis value {value} is outside [-360, 360] degrees")
124            }
125        }
126    }
127}
128
129impl std::error::Error for TecSamplesError {}
130
131impl Ionex {
132    /// Build an IONEX product directly from whole-grid samples.
133    pub fn from_samples(samples: TecGridSamples) -> core::result::Result<Self, TecSamplesError> {
134        validate_grid_samples(&samples)?;
135        Self::from_parts(IonexParts {
136            lat_nodes_deg: samples.lat_nodes_deg,
137            lon_nodes_deg: samples.lon_nodes_deg,
138            dlat_deg: samples.dlat_deg,
139            dlon_deg: samples.dlon_deg,
140            shell_height_km: samples.shell_height_km,
141            base_radius_km: samples.base_radius_km,
142            exponent: samples.exponent,
143            map_epochs: samples.map_epochs,
144            tec_maps: samples.tec_maps,
145            rms_maps: samples.rms_maps,
146            skipped_records: 0,
147        })
148        .map_err(|_| {
149            // Public sample validation mirrors Ionex::from_parts. This fallback is
150            // only for future private invariants that TecSamplesError cannot yet
151            // classify more precisely.
152            TecSamplesError::ShapeMismatch
153        })
154    }
155
156    /// Build an IONEX product from a flat stream of node samples.
157    pub fn from_node_samples(
158        samples: impl IntoIterator<Item = TecSample>,
159        shell_height_km: f64,
160        base_radius_km: f64,
161        exponent: i32,
162    ) -> core::result::Result<Self, TecSamplesError> {
163        let samples: Vec<TecSample> = samples.into_iter().collect();
164        if samples.is_empty() {
165            return Err(TecSamplesError::Empty);
166        }
167        for sample in &samples {
168            validate_axis_value(sample.lat_deg)?;
169            validate_axis_value(sample.lon_deg)?;
170            validate_finite(sample.vtec_tecu)?;
171            if let Some(rms) = sample.rms_tecu {
172                validate_finite(rms)?;
173            }
174            exact_j2000_second(sample.epoch).ok_or(TecSamplesError::EpochNotRepresentable)?;
175        }
176
177        let mut map_epochs = Vec::new();
178        let mut lat_nodes_deg = Vec::new();
179        let mut lon_nodes_deg = Vec::new();
180        for sample in &samples {
181            let epoch_s =
182                exact_j2000_second(sample.epoch).ok_or(TecSamplesError::EpochNotRepresentable)?;
183            if !map_epochs
184                .iter()
185                .any(|&epoch| exact_j2000_second(epoch) == Some(epoch_s))
186            {
187                map_epochs.push(sample.epoch);
188            }
189            push_unique_bits(&mut lat_nodes_deg, sample.lat_deg);
190            push_unique_bits(&mut lon_nodes_deg, sample.lon_deg);
191        }
192
193        map_epochs.sort_by_key(|epoch| {
194            exact_j2000_second(*epoch).expect("sample epochs were already validated")
195        });
196        lat_nodes_deg.sort_by(|a, b| b.total_cmp(a));
197        lon_nodes_deg.sort_by(f64::total_cmp);
198
199        if lat_nodes_deg.len() < 2 {
200            return Err(TecSamplesError::TooFewNodes(lat_nodes_deg.len()));
201        }
202        if lon_nodes_deg.len() < 2 {
203            return Err(TecSamplesError::TooFewNodes(lon_nodes_deg.len()));
204        }
205
206        let nmap = map_epochs.len();
207        let nlat = lat_nodes_deg.len();
208        let nlon = lon_nodes_deg.len();
209        let has_rms = samples.iter().any(|sample| sample.rms_tecu.is_some());
210        let mut tec_maps = vec![vec![vec![f64::NAN; nlon]; nlat]; nmap];
211        let mut rms_maps = if has_rms {
212            vec![vec![vec![f64::NAN; nlon]; nlat]; nmap]
213        } else {
214            Vec::new()
215        };
216        let mut seen = vec![false; nmap * nlat * nlon];
217
218        for sample in samples {
219            let map_index = map_epochs
220                .iter()
221                .position(|&epoch| exact_j2000_second(epoch) == exact_j2000_second(sample.epoch))
222                .expect("sample epoch exists in the map axis");
223            let lat_index = find_bits(&lat_nodes_deg, sample.lat_deg)
224                .expect("sample latitude exists in the latitude axis");
225            let lon_index = find_bits(&lon_nodes_deg, sample.lon_deg)
226                .expect("sample longitude exists in the longitude axis");
227            let flat_index = (map_index * nlat + lat_index) * nlon + lon_index;
228            if seen[flat_index] {
229                return Err(TecSamplesError::ShapeMismatch);
230            }
231            seen[flat_index] = true;
232            tec_maps[map_index][lat_index][lon_index] = sample.vtec_tecu;
233            if has_rms {
234                let rms = sample.rms_tecu.ok_or(TecSamplesError::RmsCountMismatch)?;
235                rms_maps[map_index][lat_index][lon_index] = rms;
236            }
237        }
238        if seen.iter().any(|&value| !value) {
239            return Err(TecSamplesError::ShapeMismatch);
240        }
241
242        Self::from_samples(TecGridSamples {
243            map_epochs,
244            dlat_deg: lat_nodes_deg[1] - lat_nodes_deg[0],
245            dlon_deg: lon_nodes_deg[1] - lon_nodes_deg[0],
246            lat_nodes_deg,
247            lon_nodes_deg,
248            shell_height_km,
249            base_radius_km,
250            exponent,
251            tec_maps,
252            rms_maps,
253        })
254    }
255
256    /// Extract this product as whole-grid IONEX samples.
257    pub fn tec_grid_samples(&self) -> TecGridSamples {
258        TecGridSamples {
259            map_epochs: self.map_epochs().to_vec(),
260            lat_nodes_deg: self.lat_nodes_deg().to_vec(),
261            lon_nodes_deg: self.lon_nodes_deg().to_vec(),
262            dlat_deg: self.dlat_deg(),
263            dlon_deg: self.dlon_deg(),
264            shell_height_km: self.shell_height_km(),
265            base_radius_km: self.base_radius_km(),
266            exponent: self.exponent(),
267            tec_maps: self.tec_maps().to_vec(),
268            rms_maps: self.rms_maps().to_vec(),
269        }
270    }
271
272    /// Extract this product as one sample per grid node.
273    pub fn tec_samples(&self) -> Vec<TecSample> {
274        let nmap = self.map_epochs().len();
275        let nlat = self.lat_nodes_deg().len();
276        let nlon = self.lon_nodes_deg().len();
277        let mut out = Vec::with_capacity(nmap * nlat * nlon);
278        let has_rms = !self.rms_maps().is_empty();
279        for (map_index, &epoch) in self.map_epochs().iter().enumerate() {
280            for (lat_index, &lat_deg) in self.lat_nodes_deg().iter().enumerate() {
281                for (lon_index, &lon_deg) in self.lon_nodes_deg().iter().enumerate() {
282                    out.push(TecSample {
283                        epoch,
284                        lat_deg,
285                        lon_deg,
286                        vtec_tecu: self.tec_maps()[map_index][lat_index][lon_index],
287                        rms_tecu: if has_rms {
288                            Some(self.rms_maps()[map_index][lat_index][lon_index])
289                        } else {
290                            None
291                        },
292                    });
293                }
294            }
295        }
296        out
297    }
298}
299
300fn validate_grid_samples(samples: &TecGridSamples) -> core::result::Result<(), TecSamplesError> {
301    if samples.map_epochs.is_empty() || samples.tec_maps.is_empty() {
302        return Err(TecSamplesError::Empty);
303    }
304    validate_axis(&samples.lat_nodes_deg, true)?;
305    validate_axis(&samples.lon_nodes_deg, false)?;
306    validate_finite(samples.dlat_deg)?;
307    validate_finite(samples.dlon_deg)?;
308    validate_axis_value(samples.dlat_deg)?;
309    validate_axis_value(samples.dlon_deg)?;
310    if samples.dlat_deg >= 0.0 || samples.dlon_deg <= 0.0 {
311        return Err(TecSamplesError::NonPositiveStep);
312    }
313    validate_finite(samples.shell_height_km)?;
314    validate_finite(samples.base_radius_km)?;
315    validate_epochs(&samples.map_epochs)?;
316
317    if samples.tec_maps.len() != samples.map_epochs.len() {
318        return Err(TecSamplesError::ShapeMismatch);
319    }
320    validate_maps(
321        &samples.tec_maps,
322        samples.map_epochs.len(),
323        samples.lat_nodes_deg.len(),
324        samples.lon_nodes_deg.len(),
325        TecSamplesError::ShapeMismatch,
326    )?;
327    if !samples.rms_maps.is_empty() {
328        if samples.rms_maps.len() != samples.map_epochs.len() {
329            return Err(TecSamplesError::RmsCountMismatch);
330        }
331        validate_maps(
332            &samples.rms_maps,
333            samples.map_epochs.len(),
334            samples.lat_nodes_deg.len(),
335            samples.lon_nodes_deg.len(),
336            TecSamplesError::RmsCountMismatch,
337        )?;
338    }
339    Ok(())
340}
341
342fn validate_axis(nodes: &[f64], descending: bool) -> core::result::Result<(), TecSamplesError> {
343    if nodes.len() < 2 {
344        return Err(TecSamplesError::TooFewNodes(nodes.len()));
345    }
346    for &node in nodes {
347        validate_axis_value(node)?;
348    }
349    if descending {
350        if nodes.windows(2).any(|w| w[1] >= w[0]) {
351            return Err(TecSamplesError::NonMonotonicLat);
352        }
353    } else if nodes.windows(2).any(|w| w[1] <= w[0]) {
354        return Err(TecSamplesError::NonMonotonicLon);
355    }
356    Ok(())
357}
358
359fn validate_epochs(map_epochs: &[Instant]) -> core::result::Result<(), TecSamplesError> {
360    let mut previous_s = None;
361    for &epoch in map_epochs {
362        let seconds = exact_j2000_second(epoch).ok_or(TecSamplesError::EpochNotRepresentable)?;
363        if previous_s.is_some_and(|previous| seconds <= previous) {
364            return Err(TecSamplesError::NonMonotonicEpochs);
365        }
366        previous_s = Some(seconds);
367    }
368    Ok(())
369}
370
371fn validate_maps(
372    maps: &[Vec<Vec<f64>>],
373    expected_maps: usize,
374    expected_lat: usize,
375    expected_lon: usize,
376    dimension_error: TecSamplesError,
377) -> core::result::Result<(), TecSamplesError> {
378    if maps.len() != expected_maps {
379        return Err(dimension_error.clone());
380    }
381    for map in maps {
382        if map.len() != expected_lat {
383            return Err(dimension_error.clone());
384        }
385        for row in map {
386            if row.len() != expected_lon {
387                return Err(dimension_error.clone());
388            }
389            for &value in row {
390                validate_finite(value)?;
391            }
392        }
393    }
394    Ok(())
395}
396
397fn validate_axis_value(value: f64) -> core::result::Result<(), TecSamplesError> {
398    validate_finite(value)?;
399    if !(-IONEX_AXIS_DEG_LIMIT..=IONEX_AXIS_DEG_LIMIT).contains(&value) {
400        return Err(TecSamplesError::AxisOutOfRange(value));
401    }
402    Ok(())
403}
404
405fn validate_finite(value: f64) -> core::result::Result<(), TecSamplesError> {
406    if value.is_finite() {
407        Ok(())
408    } else {
409        Err(TecSamplesError::NonFiniteValue)
410    }
411}
412
413fn exact_j2000_second(epoch: Instant) -> Option<i64> {
414    match epoch.repr {
415        InstantRepr::Nanos(nanos) => {
416            if nanos % NANOS_PER_SECOND != 0 {
417                return None;
418            }
419            let seconds = nanos / NANOS_PER_SECOND;
420            i64::try_from(seconds).ok()
421        }
422        InstantRepr::JulianDate(_) => j2000_seconds_from_instant(epoch),
423    }
424}
425
426fn push_unique_bits(values: &mut Vec<f64>, value: f64) {
427    if !values
428        .iter()
429        .any(|&existing| same_axis_node(existing, value))
430    {
431        values.push(value);
432    }
433}
434
435fn find_bits(values: &[f64], value: f64) -> Option<usize> {
436    values
437        .iter()
438        .position(|&existing| same_axis_node(existing, value))
439}
440
441fn same_axis_node(a: f64, b: f64) -> bool {
442    a == b || a.to_bits() == b.to_bits()
443}
sidereon_core/ionex/samples.rs

sidereon_core/ionex/
samples.rs
