sidereon_core/
observation_qc.rs

1//! RINEX observation-file quality-control rollups.
2//!
3//! This module works from an already parsed [`RinexObs`] product. It does not
4//! parse, repair, or resample files; it reports the completeness and signal
5//! indicators a caller needs before choosing solver inputs.
6
7use std::collections::BTreeMap;
8
9use crate::id::{GnssSatelliteId, GnssSystem};
10use crate::rinex::observations::{ObsEpochTime, RinexObs};
11use crate::rinex_common::{dominant_obs_interval_s, obs_epoch_seconds};
12
13/// Options controlling RINEX observation QC aggregation.
14#[derive(Debug, Clone, Copy, PartialEq)]
15pub struct ObservationQcOptions {
16    /// Override the header `INTERVAL` value when detecting missing epochs.
17    pub interval_override_s: Option<f64>,
18    /// Minimum `delta / interval` ratio that is treated as a data gap.
19    pub gap_factor: f64,
20}
21
22impl Default for ObservationQcOptions {
23    fn default() -> Self {
24        Self {
25            interval_override_s: None,
26            gap_factor: 1.5,
27        }
28    }
29}
30
31/// Error returned when QC options are invalid.
32#[derive(Debug, Clone, Copy, PartialEq, thiserror::Error)]
33pub enum ObservationQcError {
34    /// The supplied nominal interval was zero, negative, or non-finite.
35    #[error("invalid observation QC interval: must be finite and positive")]
36    InvalidInterval,
37    /// The supplied gap factor was not finite and greater than one.
38    #[error("invalid observation QC gap factor: must be finite and greater than one")]
39    InvalidGapFactor,
40}
41
42/// Source of the interval used for gap detection.
43#[derive(Debug, Clone, Copy, PartialEq, Eq)]
44pub enum IntervalSource {
45    /// Caller override.
46    Override,
47    /// Header `INTERVAL`.
48    Header,
49    /// Modal positive epoch delta inferred from the body.
50    Inferred,
51    /// Not enough positive epoch deltas were available.
52    Unresolved,
53}
54
55/// Non-fatal QC note.
56#[derive(Debug, Clone, Copy, PartialEq, Eq)]
57pub enum ObservationQcNote {
58    /// Adjacent observation epochs were duplicate or out of order.
59    NonMonotonicEpoch { epoch_index: usize },
60    /// No interval could be resolved.
61    IntervalUnresolved,
62}
63
64/// Aggregate QC report for one parsed RINEX observation file.
65#[derive(Debug, Clone, PartialEq)]
66pub struct ObservationQcReport {
67    /// Total number of epoch records retained by the parser, including events.
68    pub total_epoch_records: usize,
69    /// Count of normal observation epochs (`flag == 0`) and power-failure
70    /// observation epochs (`flag == 1`).
71    pub observation_epochs: usize,
72    /// Count of non-observation event records (`flag > 1`).
73    pub event_records: usize,
74    /// Count of observation epochs marked as power-failure epochs (`flag == 1`).
75    pub power_failure_epochs: usize,
76    /// Count of malformed records skipped by the RINEX observation parser.
77    pub skipped_records: usize,
78    /// Interval used for gap detection.
79    pub interval_s: Option<f64>,
80    /// Where `interval_s` came from.
81    pub interval_source: IntervalSource,
82    /// Estimated number of missing nominal epochs across all detected gaps.
83    pub missing_epochs: usize,
84    /// Gaps detected from adjacent observation epochs and the nominal interval.
85    pub data_gaps: Vec<ObservationDataGap>,
86    /// Per-satellite observation completeness.
87    pub satellites: Vec<SatelliteObservationQc>,
88    /// Per-satellite, per-code observation completeness and SSI statistics.
89    pub satellite_signals: Vec<SatelliteSignalQc>,
90    /// Per-system, per-code observation completeness and SSI statistics.
91    pub system_signals: Vec<SystemSignalQc>,
92    /// Non-fatal QC notes.
93    pub notes: Vec<ObservationQcNote>,
94}
95
96/// One detected gap between adjacent observation epochs.
97#[derive(Debug, Clone, PartialEq)]
98pub struct ObservationDataGap {
99    /// Epoch immediately before the gap.
100    pub start_epoch: ObsEpochTime,
101    /// Epoch immediately after the gap.
102    pub end_epoch: ObsEpochTime,
103    /// Nominal interval used for the estimate.
104    pub nominal_interval_s: f64,
105    /// Observed delta between the two retained epochs.
106    pub observed_delta_s: f64,
107    /// Estimated missing nominal epochs between `start_epoch` and `end_epoch`.
108    pub missing_epochs: usize,
109}
110
111/// Per-satellite observation counts.
112#[derive(Debug, Clone, PartialEq, Eq)]
113pub struct SatelliteObservationQc {
114    /// Satellite id.
115    pub satellite: GnssSatelliteId,
116    /// Epochs where the satellite has at least one non-blank observation value.
117    pub epochs_with_observations: usize,
118    /// Non-blank observation values across all codes and observation epochs.
119    pub value_observations: usize,
120}
121
122/// Per-satellite, per-observation-code counts.
123#[derive(Debug, Clone, PartialEq)]
124pub struct SatelliteSignalQc {
125    /// Satellite id.
126    pub satellite: GnssSatelliteId,
127    /// RINEX observation code, e.g. `C1C`, `L1C`, or `S1C`.
128    pub code: String,
129    /// Non-blank values for this satellite/code pair.
130    pub value_observations: usize,
131    /// Signal-strength indicator statistics for non-blank values that carried
132    /// an SSI digit.
133    pub ssi: Option<SsiHistogram>,
134    /// Raw S-code statistics when this code is an `S*` observable.
135    pub snr: Option<SnrStats>,
136}
137
138/// Per-system, per-observation-code counts.
139#[derive(Debug, Clone, PartialEq)]
140pub struct SystemSignalQc {
141    /// GNSS constellation.
142    pub system: GnssSystem,
143    /// RINEX observation code, e.g. `C1C`, `L1C`, or `S1C`.
144    pub code: String,
145    /// Non-blank values for this system/code pair across satellites.
146    pub value_observations: usize,
147    /// Signal-strength indicator statistics for non-blank values that carried
148    /// an SSI digit.
149    pub ssi: Option<SsiHistogram>,
150    /// Raw S-code statistics when this code is an `S*` observable.
151    pub snr: Option<SnrStats>,
152}
153
154/// Histogram over RINEX SSI digits. Index 0 is blank/unknown.
155#[derive(Debug, Clone, Copy, PartialEq, Eq)]
156pub struct SsiHistogram {
157    /// Counts indexed by SSI digit.
158    pub counts: [u64; 10],
159}
160
161/// Summary statistics over raw numeric `S*` signal-strength observations.
162#[derive(Debug, Clone, Copy, PartialEq)]
163pub struct SnrStats {
164    /// Number of samples.
165    pub n: usize,
166    /// Arithmetic mean.
167    pub mean: f64,
168    /// Minimum sample.
169    pub min: f64,
170    /// Maximum sample.
171    pub max: f64,
172    /// Sample standard deviation, absent for one sample.
173    pub std: Option<f64>,
174}
175
176/// Build a QC report with default options.
177pub fn observation_qc(obs: &RinexObs) -> ObservationQcReport {
178    observation_qc_with_options(obs, ObservationQcOptions::default())
179        .expect("default observation QC options are valid")
180}
181
182/// Build a QC report with explicit options.
183pub fn observation_qc_with_options(
184    obs: &RinexObs,
185    options: ObservationQcOptions,
186) -> Result<ObservationQcReport, ObservationQcError> {
187    validate_options(options)?;
188
189    let mut satellites: BTreeMap<GnssSatelliteId, SatelliteAccum> = BTreeMap::new();
190    let mut satellite_signals: BTreeMap<(GnssSatelliteId, String), SignalAccum> = BTreeMap::new();
191    let mut system_signals: BTreeMap<(GnssSystem, String), SignalAccum> = BTreeMap::new();
192    let mut observation_epoch_times = Vec::new();
193
194    let mut observation_epochs = 0;
195    let mut event_records = 0;
196    let mut power_failure_epochs = 0;
197
198    for epoch in obs.epochs() {
199        if epoch.flag > 1 {
200            event_records += 1;
201            continue;
202        }
203
204        observation_epochs += 1;
205        if epoch.flag == 1 {
206            power_failure_epochs += 1;
207        }
208        observation_epoch_times.push(epoch.epoch);
209
210        for (satellite, values) in &epoch.sats {
211            let value_observations = values.iter().filter(|value| value.value.is_some()).count();
212            if value_observations == 0 {
213                continue;
214            }
215
216            let satellite_acc = satellites.entry(*satellite).or_default();
217            satellite_acc.epochs_with_observations += 1;
218            satellite_acc.value_observations += value_observations;
219
220            let Some(codes) = obs.header().obs_codes.get(&satellite.system) else {
221                continue;
222            };
223
224            for (index, value) in values.iter().enumerate() {
225                if value.value.is_none() {
226                    continue;
227                }
228
229                let Some(code) = codes.get(index) else {
230                    continue;
231                };
232
233                let sat_signal = satellite_signals
234                    .entry((*satellite, code.clone()))
235                    .or_default();
236                sat_signal.add(code, value.value, value.ssi);
237
238                let sys_signal = system_signals
239                    .entry((satellite.system, code.clone()))
240                    .or_default();
241                sys_signal.add(code, value.value, value.ssi);
242            }
243        }
244    }
245
246    let mut notes = non_monotonic_notes(&observation_epoch_times);
247    let (interval_s, interval_source) =
248        resolve_interval(obs, options, &observation_epoch_times, &mut notes)?;
249    let data_gaps = detect_gaps(options, &observation_epoch_times, interval_s)?;
250    let missing_epochs = data_gaps.iter().map(|gap| gap.missing_epochs).sum();
251
252    Ok(ObservationQcReport {
253        total_epoch_records: obs.epochs().len(),
254        observation_epochs,
255        event_records,
256        power_failure_epochs,
257        skipped_records: obs.skipped_records,
258        interval_s,
259        interval_source,
260        missing_epochs,
261        data_gaps,
262        satellites: satellites
263            .into_iter()
264            .map(|(satellite, acc)| SatelliteObservationQc {
265                satellite,
266                epochs_with_observations: acc.epochs_with_observations,
267                value_observations: acc.value_observations,
268            })
269            .collect(),
270        satellite_signals: satellite_signals
271            .into_iter()
272            .map(|((satellite, code), acc)| SatelliteSignalQc {
273                satellite,
274                code,
275                value_observations: acc.value_observations,
276                ssi: acc.ssi.finish(),
277                snr: acc.snr.finish(),
278            })
279            .collect(),
280        system_signals: system_signals
281            .into_iter()
282            .map(|((system, code), acc)| SystemSignalQc {
283                system,
284                code,
285                value_observations: acc.value_observations,
286                ssi: acc.ssi.finish(),
287                snr: acc.snr.finish(),
288            })
289            .collect(),
290        notes,
291    })
292}
293
294fn validate_options(options: ObservationQcOptions) -> Result<(), ObservationQcError> {
295    if !options.gap_factor.is_finite() || options.gap_factor <= 1.0 {
296        return Err(ObservationQcError::InvalidGapFactor);
297    }
298
299    if let Some(interval_s) = options.interval_override_s {
300        validate_interval(interval_s)?;
301    }
302
303    Ok(())
304}
305
306fn validate_interval(interval_s: f64) -> Result<(), ObservationQcError> {
307    if interval_s.is_finite() && interval_s > 0.0 {
308        Ok(())
309    } else {
310        Err(ObservationQcError::InvalidInterval)
311    }
312}
313
314fn resolve_interval(
315    obs: &RinexObs,
316    options: ObservationQcOptions,
317    observation_epoch_times: &[ObsEpochTime],
318    notes: &mut Vec<ObservationQcNote>,
319) -> Result<(Option<f64>, IntervalSource), ObservationQcError> {
320    let Some(interval_s) = options.interval_override_s else {
321        if let Some(interval_s) = obs.header().interval_s {
322            validate_interval(interval_s)?;
323            return Ok((Some(interval_s), IntervalSource::Header));
324        }
325        if let Some(interval_s) = dominant_obs_interval_s(observation_epoch_times) {
326            return Ok((Some(interval_s), IntervalSource::Inferred));
327        }
328        notes.push(ObservationQcNote::IntervalUnresolved);
329        return Ok((None, IntervalSource::Unresolved));
330    };
331    validate_interval(interval_s)?;
332    Ok((Some(interval_s), IntervalSource::Override))
333}
334
335fn detect_gaps(
336    options: ObservationQcOptions,
337    observation_epoch_times: &[ObsEpochTime],
338    interval_s: Option<f64>,
339) -> Result<Vec<ObservationDataGap>, ObservationQcError> {
340    let Some(interval_s) = interval_s else {
341        return Ok(Vec::new());
342    };
343
344    let mut gaps = Vec::new();
345    for window in observation_epoch_times.windows(2) {
346        let start_epoch = window[0];
347        let end_epoch = window[1];
348        let observed_delta_s = obs_epoch_seconds(end_epoch) - obs_epoch_seconds(start_epoch);
349        if observed_delta_s <= 0.0 || observed_delta_s <= interval_s * options.gap_factor {
350            continue;
351        }
352
353        let missing_epochs = ((observed_delta_s / interval_s).round() as isize - 1) as usize;
354        gaps.push(ObservationDataGap {
355            start_epoch,
356            end_epoch,
357            nominal_interval_s: interval_s,
358            observed_delta_s,
359            missing_epochs,
360        });
361    }
362
363    Ok(gaps)
364}
365
366fn non_monotonic_notes(observation_epoch_times: &[ObsEpochTime]) -> Vec<ObservationQcNote> {
367    let mut notes = Vec::new();
368    for (idx, window) in observation_epoch_times.windows(2).enumerate() {
369        if obs_epoch_seconds(window[1]) - obs_epoch_seconds(window[0]) <= 0.0 {
370            notes.push(ObservationQcNote::NonMonotonicEpoch {
371                epoch_index: idx + 1,
372            });
373        }
374    }
375    notes
376}
377
378#[derive(Debug, Default)]
379struct SatelliteAccum {
380    epochs_with_observations: usize,
381    value_observations: usize,
382}
383
384#[derive(Debug, Default)]
385struct SignalAccum {
386    value_observations: usize,
387    ssi: SsiAccum,
388    snr: SnrAccum,
389}
390
391impl SignalAccum {
392    fn add(&mut self, code: &str, value: Option<f64>, ssi: Option<u8>) {
393        self.value_observations += 1;
394        self.ssi.add(ssi);
395        if code.starts_with('S') {
396            if let Some(value) = value {
397                self.snr.add(value);
398            }
399        }
400    }
401}
402
403#[derive(Debug, Default)]
404struct SsiAccum {
405    counts: [u64; 10],
406}
407
408impl SsiAccum {
409    fn add(&mut self, value: Option<u8>) {
410        let idx = value.unwrap_or(0).min(9) as usize;
411        self.counts[idx] += 1;
412    }
413
414    fn finish(self) -> Option<SsiHistogram> {
415        if self.counts.iter().all(|count| *count == 0) {
416            return None;
417        }
418
419        Some(SsiHistogram {
420            counts: self.counts,
421        })
422    }
423}
424
425#[derive(Debug, Default)]
426struct SnrAccum {
427    samples: Vec<f64>,
428}
429
430impl SnrAccum {
431    fn add(&mut self, value: f64) {
432        self.samples.push(value);
433    }
434
435    fn finish(self) -> Option<SnrStats> {
436        if self.samples.is_empty() {
437            return None;
438        }
439        let n = self.samples.len();
440        let mean = self.samples.iter().sum::<f64>() / n as f64;
441        let min = self.samples.iter().copied().fold(f64::INFINITY, f64::min);
442        let max = self
443            .samples
444            .iter()
445            .copied()
446            .fold(f64::NEG_INFINITY, f64::max);
447        let std = (n > 1).then(|| {
448            let sum_sq = self
449                .samples
450                .iter()
451                .map(|value| {
452                    let residual = *value - mean;
453                    residual * residual
454                })
455                .sum::<f64>();
456            (sum_sq / (n - 1) as f64).sqrt()
457        });
458        Some(SnrStats {
459            n,
460            mean,
461            min,
462            max,
463            std,
464        })
465    }
466}
467
468#[cfg(test)]
469mod tests {
470    use super::*;
471    use crate::rinex::observations::{ObsEpoch, ObsHeader, ObsValue};
472    use serde_json::Value;
473    use std::collections::BTreeMap;
474    use std::path::PathBuf;
475
476    #[test]
477    fn observation_qc_counts_epochs_satellites_signals_and_ssi() {
478        let g01 = sat(1);
479        let g02 = sat(2);
480        let obs = observation_file(vec![
481            epoch(
482                0,
483                0.0,
484                0,
485                BTreeMap::from([
486                    (
487                        g01,
488                        vec![
489                            obs_value(Some(1.0), Some(5)),
490                            obs_value(Some(2.0), Some(6)),
491                            obs_value(None, None),
492                        ],
493                    ),
494                    (
495                        g02,
496                        vec![
497                            obs_value(Some(10.0), Some(4)),
498                            obs_value(None, None),
499                            obs_value(None, None),
500                        ],
501                    ),
502                ]),
503            ),
504            epoch(
505                0,
506                30.0,
507                1,
508                BTreeMap::from([(
509                    g01,
510                    vec![
511                        obs_value(Some(3.0), Some(7)),
512                        obs_value(None, None),
513                        obs_value(Some(9.0), Some(8)),
514                    ],
515                )]),
516            ),
517            epoch(1, 0.0, 2, BTreeMap::new()),
518        ]);
519
520        let report = observation_qc(&obs);
521
522        assert_eq!(report.total_epoch_records, 3);
523        assert_eq!(report.observation_epochs, 2);
524        assert_eq!(report.event_records, 1);
525        assert_eq!(report.power_failure_epochs, 1);
526        assert_eq!(report.skipped_records, 0);
527        assert_eq!(report.satellites.len(), 2);
528        assert_eq!(
529            report.satellites[0],
530            SatelliteObservationQc {
531                satellite: g01,
532                epochs_with_observations: 2,
533                value_observations: 4,
534            }
535        );
536        assert_eq!(
537            report.satellites[1],
538            SatelliteObservationQc {
539                satellite: g02,
540                epochs_with_observations: 1,
541                value_observations: 1,
542            }
543        );
544
545        let g01_c1c = report
546            .satellite_signals
547            .iter()
548            .find(|signal| signal.satellite == g01 && signal.code == "C1C")
549            .expect("G01 C1C signal");
550        assert_eq!(g01_c1c.value_observations, 2);
551        assert_eq!(
552            g01_c1c.ssi,
553            Some(SsiHistogram {
554                counts: [0, 0, 0, 0, 0, 1, 0, 1, 0, 0],
555            })
556        );
557        assert_eq!(g01_c1c.snr, None);
558
559        let gps_c1c = report
560            .system_signals
561            .iter()
562            .find(|signal| signal.system == GnssSystem::Gps && signal.code == "C1C")
563            .expect("GPS C1C signal");
564        assert_eq!(gps_c1c.value_observations, 3);
565        assert_eq!(
566            gps_c1c.ssi,
567            Some(SsiHistogram {
568                counts: [0, 0, 0, 0, 1, 1, 0, 1, 0, 0],
569            })
570        );
571
572        let gps_s1c = report
573            .system_signals
574            .iter()
575            .find(|signal| signal.system == GnssSystem::Gps && signal.code == "S1C")
576            .expect("GPS S1C signal");
577        assert_eq!(
578            gps_s1c.snr,
579            Some(SnrStats {
580                n: 1,
581                mean: 9.0,
582                min: 9.0,
583                max: 9.0,
584                std: None,
585            })
586        );
587    }
588
589    #[test]
590    fn observation_qc_detects_nominal_interval_gaps() {
591        let g01 = sat(1);
592        let obs = observation_file(vec![
593            epoch(
594                0,
595                0.0,
596                0,
597                BTreeMap::from([(g01, vec![obs_value(Some(1.0), Some(5))])]),
598            ),
599            epoch(
600                1,
601                30.0,
602                0,
603                BTreeMap::from([(g01, vec![obs_value(Some(2.0), Some(6))])]),
604            ),
605        ]);
606
607        let report = observation_qc(&obs);
608
609        assert_eq!(report.missing_epochs, 2);
610        assert_eq!(report.data_gaps.len(), 1);
611        assert_eq!(report.data_gaps[0].nominal_interval_s, 30.0);
612        assert_eq!(report.data_gaps[0].observed_delta_s, 90.0);
613        assert_eq!(report.data_gaps[0].missing_epochs, 2);
614    }
615
616    #[test]
617    fn observation_qc_infers_interval_when_header_is_absent() {
618        let g01 = sat(1);
619        let mut obs = observation_file(vec![
620            epoch(
621                0,
622                0.0,
623                0,
624                BTreeMap::from([(g01, vec![obs_value(Some(1.0), Some(5))])]),
625            ),
626            epoch(
627                0,
628                30.0,
629                0,
630                BTreeMap::from([(g01, vec![obs_value(Some(2.0), Some(6))])]),
631            ),
632            epoch(
633                2,
634                0.0,
635                0,
636                BTreeMap::from([(g01, vec![obs_value(Some(3.0), Some(7))])]),
637            ),
638        ]);
639        obs.header.interval_s = None;
640
641        let report = observation_qc(&obs);
642
643        assert_eq!(report.interval_s, Some(30.0));
644        assert_eq!(report.interval_source, IntervalSource::Inferred);
645        assert_eq!(report.missing_epochs, 2);
646    }
647
648    #[test]
649    fn observation_qc_notes_non_monotonic_epochs_and_excludes_them_from_gaps() {
650        let g01 = sat(1);
651        let obs = observation_file(vec![
652            epoch(
653                1,
654                0.0,
655                0,
656                BTreeMap::from([(g01, vec![obs_value(Some(1.0), Some(5))])]),
657            ),
658            epoch(
659                0,
660                30.0,
661                0,
662                BTreeMap::from([(g01, vec![obs_value(Some(2.0), Some(6))])]),
663            ),
664        ]);
665
666        let report = observation_qc(&obs);
667
668        assert_eq!(
669            report.notes,
670            vec![ObservationQcNote::NonMonotonicEpoch { epoch_index: 1 }]
671        );
672        assert!(report.data_gaps.is_empty());
673    }
674
675    #[test]
676    fn observation_qc_rejects_invalid_options() {
677        let obs = observation_file(Vec::new());
678
679        let err = observation_qc_with_options(
680            &obs,
681            ObservationQcOptions {
682                interval_override_s: Some(0.0),
683                gap_factor: 1.5,
684            },
685        )
686        .expect_err("invalid interval");
687        assert_eq!(err, ObservationQcError::InvalidInterval);
688
689        let err = observation_qc_with_options(
690            &obs,
691            ObservationQcOptions {
692                interval_override_s: None,
693                gap_factor: 1.0,
694            },
695        )
696        .expect_err("invalid gap factor");
697        assert_eq!(err, ObservationQcError::InvalidGapFactor);
698    }
699
700    #[test]
701    fn observation_qc_matches_independent_real_fixture_oracles() {
702        let doc = read_json_fixture("qc/observation_qc_real_oracles.json");
703        assert_eq!(
704            doc["provenance"]["generator"],
705            "crates/sidereon-core/fixtures-generators/generate_observation_qc_oracles.py"
706        );
707        for fixture in doc["fixtures"].as_array().expect("fixtures array") {
708            let rel = fixture["path"].as_str().expect("fixture path");
709            let text = std::fs::read_to_string(fixture_path(rel))
710                .unwrap_or_else(|e| panic!("read {rel}: {e}"));
711            let obs = RinexObs::parse(&text).unwrap_or_else(|e| panic!("parse {rel}: {e}"));
712            let report = observation_qc(&obs);
713
714            assert_eq!(
715                report.total_epoch_records,
716                fixture["total_epoch_records"].as_u64().unwrap() as usize,
717                "{rel}"
718            );
719            assert_eq!(
720                report.observation_epochs,
721                fixture["observation_epochs"].as_u64().unwrap() as usize,
722                "{rel}"
723            );
724            assert_eq!(
725                report.event_records,
726                fixture["event_records"].as_u64().unwrap() as usize,
727                "{rel}"
728            );
729            assert_eq!(
730                report.power_failure_epochs,
731                fixture["power_failure_epochs"].as_u64().unwrap() as usize,
732                "{rel}"
733            );
734            assert_eq!(
735                report.skipped_records,
736                fixture["skipped_records"].as_u64().unwrap() as usize,
737                "{rel}"
738            );
739            assert_close(
740                report.interval_s.expect("oracle interval"),
741                fixture["interval_s"].as_f64().unwrap(),
742                rel,
743            );
744            assert_eq!(
745                report.missing_epochs,
746                fixture["missing_epochs"].as_u64().unwrap() as usize,
747                "{rel}"
748            );
749            assert_gaps(&report.data_gaps, &fixture["data_gaps"], rel);
750            assert_satellites(&report.satellites, &fixture["satellites"], rel);
751            assert_satellite_signals(
752                &report.satellite_signals,
753                &fixture["satellite_signals"],
754                rel,
755            );
756            assert_system_signals(&report.system_signals, &fixture["system_signals"], rel);
757        }
758    }
759
760    fn observation_file(epochs: Vec<ObsEpoch>) -> RinexObs {
761        RinexObs {
762            header: ObsHeader {
763                version: 3.05,
764                approx_position_m: None,
765                antenna_delta_hen_m: None,
766                obs_codes: BTreeMap::from([(
767                    GnssSystem::Gps,
768                    vec!["C1C".to_string(), "L1C".to_string(), "S1C".to_string()],
769                )]),
770                program_run_by_date: None,
771                comments: Vec::new(),
772                marker_number: None,
773                marker_type: None,
774                observer: None,
775                agency: None,
776                receiver: None,
777                antenna: None,
778                interval_s: Some(30.0),
779                time_of_first_obs: None,
780                time_of_last_obs: None,
781                n_satellites: None,
782                prn_obs_counts: BTreeMap::new(),
783                phase_shifts: Vec::new(),
784                scale_factors: Vec::new(),
785                glonass_slots: BTreeMap::new(),
786                glonass_cod_phs_bis: None,
787                signal_strength_unit: None,
788                leap_seconds: None,
789                marker_name: None,
790                unretained_header_labels: Vec::new(),
791            },
792            epochs,
793            skipped_records: 0,
794        }
795    }
796
797    fn epoch(
798        minute: u8,
799        second: f64,
800        flag: u8,
801        sats: BTreeMap<GnssSatelliteId, Vec<ObsValue>>,
802    ) -> ObsEpoch {
803        ObsEpoch {
804            epoch: ObsEpochTime {
805                year: 2024,
806                month: 1,
807                day: 1,
808                hour: 0,
809                minute,
810                second,
811            },
812            flag,
813            rcv_clock_offset_s: None,
814            epoch_picoseconds: None,
815            declared_record_count: sats.len(),
816            special_record_count: if flag > 1 { sats.len() } else { 0 },
817            sats,
818        }
819    }
820
821    fn obs_value(value: Option<f64>, ssi: Option<u8>) -> ObsValue {
822        ObsValue {
823            value,
824            lli: None,
825            ssi,
826        }
827    }
828
829    fn sat(prn: u8) -> GnssSatelliteId {
830        GnssSatelliteId::new(GnssSystem::Gps, prn).expect("valid GPS PRN")
831    }
832
833    fn fixture_path(rel: &str) -> PathBuf {
834        PathBuf::from(env!("CARGO_MANIFEST_DIR")).join(rel)
835    }
836
837    fn read_json_fixture(rel: &str) -> Value {
838        let path = fixture_path(&format!("tests/fixtures/{rel}"));
839        let raw = std::fs::read_to_string(&path)
840            .unwrap_or_else(|e| panic!("read {}: {e}", path.display()));
841        serde_json::from_str(&raw).unwrap_or_else(|e| panic!("parse {}: {e}", path.display()))
842    }
843
844    fn assert_close(actual: f64, expected: f64, context: &str) {
845        assert!(
846            (actual - expected).abs() <= 1.0e-9,
847            "{context}: actual {actual:?}, expected {expected:?}"
848        );
849    }
850
851    fn assert_gaps(actual: &[ObservationDataGap], expected: &Value, context: &str) {
852        let expected = expected.as_array().expect("gap array");
853        assert_eq!(actual.len(), expected.len(), "{context}");
854        for (actual, expected) in actual.iter().zip(expected) {
855            assert_epoch(&actual.start_epoch, &expected["start_epoch"], context);
856            assert_epoch(&actual.end_epoch, &expected["end_epoch"], context);
857            assert_close(
858                actual.nominal_interval_s,
859                expected["nominal_interval_s"].as_f64().unwrap(),
860                context,
861            );
862            assert_close(
863                actual.observed_delta_s,
864                expected["observed_delta_s"].as_f64().unwrap(),
865                context,
866            );
867            assert_eq!(
868                actual.missing_epochs,
869                expected["missing_epochs"].as_u64().unwrap() as usize,
870                "{context}"
871            );
872        }
873    }
874
875    fn assert_epoch(actual: &ObsEpochTime, expected: &Value, context: &str) {
876        assert_eq!(
877            actual.year,
878            expected["year"].as_i64().unwrap() as i32,
879            "{context}"
880        );
881        assert_eq!(
882            actual.month,
883            expected["month"].as_u64().unwrap() as u8,
884            "{context}"
885        );
886        assert_eq!(
887            actual.day,
888            expected["day"].as_u64().unwrap() as u8,
889            "{context}"
890        );
891        assert_eq!(
892            actual.hour,
893            expected["hour"].as_u64().unwrap() as u8,
894            "{context}"
895        );
896        assert_eq!(
897            actual.minute,
898            expected["minute"].as_u64().unwrap() as u8,
899            "{context}"
900        );
901        assert_close(actual.second, expected["second"].as_f64().unwrap(), context);
902    }
903
904    fn assert_satellites(actual: &[SatelliteObservationQc], expected: &Value, context: &str) {
905        let expected = expected.as_array().expect("satellites array");
906        assert_eq!(actual.len(), expected.len(), "{context}");
907        let actual = actual
908            .iter()
909            .map(|sat| {
910                (
911                    sat.satellite.to_string(),
912                    (sat.epochs_with_observations, sat.value_observations),
913                )
914            })
915            .collect::<BTreeMap<_, _>>();
916        for expected in expected {
917            let satellite = expected["satellite"].as_str().unwrap();
918            let actual = actual
919                .get(satellite)
920                .unwrap_or_else(|| panic!("{context}: missing satellite {satellite}"));
921            assert_eq!(
922                actual.0,
923                expected["epochs_with_observations"].as_u64().unwrap() as usize,
924                "{context} {satellite}"
925            );
926            assert_eq!(
927                actual.1,
928                expected["value_observations"].as_u64().unwrap() as usize,
929                "{context} {satellite}"
930            );
931        }
932    }
933
934    fn assert_satellite_signals(actual: &[SatelliteSignalQc], expected: &Value, context: &str) {
935        let expected = expected.as_array().expect("satellite signals array");
936        assert_eq!(actual.len(), expected.len(), "{context}");
937        let actual = actual
938            .iter()
939            .map(|signal| {
940                (
941                    (signal.satellite.to_string(), signal.code.as_str()),
942                    (signal.value_observations, signal.ssi, signal.snr),
943                )
944            })
945            .collect::<BTreeMap<_, _>>();
946        for expected in expected {
947            let satellite = expected["satellite"].as_str().unwrap();
948            let code = expected["code"].as_str().unwrap();
949            let actual = actual
950                .get(&(satellite.to_string(), code))
951                .unwrap_or_else(|| panic!("{context}: missing {satellite} {code}"));
952            assert_eq!(
953                actual.0,
954                expected["value_observations"].as_u64().unwrap() as usize,
955                "{context} {satellite} {code}"
956            );
957            assert_ssi(actual.1, &expected["ssi"], context);
958            assert_snr(actual.2, &expected["snr"], context);
959        }
960    }
961
962    fn assert_system_signals(actual: &[SystemSignalQc], expected: &Value, context: &str) {
963        let expected = expected.as_array().expect("system signals array");
964        assert_eq!(actual.len(), expected.len(), "{context}");
965        let actual = actual
966            .iter()
967            .map(|signal| {
968                (
969                    (signal.system.letter().to_string(), signal.code.as_str()),
970                    (signal.value_observations, signal.ssi, signal.snr),
971                )
972            })
973            .collect::<BTreeMap<_, _>>();
974        for expected in expected {
975            let system = expected["system"].as_str().unwrap();
976            let code = expected["code"].as_str().unwrap();
977            let actual = actual
978                .get(&(system.to_string(), code))
979                .unwrap_or_else(|| panic!("{context}: missing {system} {code}"));
980            assert_eq!(
981                actual.0,
982                expected["value_observations"].as_u64().unwrap() as usize,
983                "{context} {system} {code}"
984            );
985            assert_ssi(actual.1, &expected["ssi"], context);
986            assert_snr(actual.2, &expected["snr"], context);
987        }
988    }
989
990    fn assert_ssi(actual: Option<SsiHistogram>, expected: &Value, context: &str) {
991        if expected.is_null() {
992            assert_eq!(actual, None, "{context}");
993            return;
994        }
995        let expected = expected
996            .as_array()
997            .expect("ssi array")
998            .iter()
999            .map(|value| value.as_u64().unwrap())
1000            .collect::<Vec<_>>();
1001        assert_eq!(actual.expect("ssi").counts.to_vec(), expected, "{context}");
1002    }
1003
1004    fn assert_snr(actual: Option<SnrStats>, expected: &Value, context: &str) {
1005        if expected.is_null() {
1006            assert_eq!(actual, None, "{context}");
1007            return;
1008        }
1009        let actual = actual.expect("snr");
1010        assert_eq!(
1011            actual.n,
1012            expected["n"].as_u64().unwrap() as usize,
1013            "{context}"
1014        );
1015        assert_close(actual.mean, expected["mean"].as_f64().unwrap(), context);
1016        assert_close(actual.min, expected["min"].as_f64().unwrap(), context);
1017        assert_close(actual.max, expected["max"].as_f64().unwrap(), context);
1018        if expected["std"].is_null() {
1019            assert_eq!(actual.std, None, "{context}");
1020        } else {
1021            assert_close(
1022                actual.std.expect("std"),
1023                expected["std"].as_f64().unwrap(),
1024                context,
1025            );
1026        }
1027    }
1028}
sidereon_core/observation_qc.rs

sidereon_core/
observation_qc.rs
