From be4b8b5f9b32979d6c4f601d8457beb1e29de1be Mon Sep 17 00:00:00 2001 From: gwbres Date: Thu, 7 Sep 2023 15:05:55 +0200 Subject: [PATCH] Observation Data QC improvements (#158) --------- Signed-off-by: Guillaume W. Bres --- rinex-cli/config/advanced_repair.json | 12 - rinex-cli/config/advanced_study.json | 12 - rinex-cli/config/basic_manual_gap.json | 7 - .../config/{basic.json => gnss_snr30db.json} | 0 rinex-cli/config/sv_manual_gap.json | 7 + rinex-cli/doc/qc.md | 121 ++-- rinex-cli/src/cli.rs | 1 + rinex-cli/src/main.rs | 26 +- rinex/Cargo.toml | 4 +- rinex/src/algorithm/averaging.rs | 67 +- rinex/src/algorithm/derivative.rs | 38 ++ rinex/src/algorithm/mod.rs | 5 + rinex/src/hardware.rs | 4 +- rinex/src/header.rs | 30 +- rinex/src/lib.rs | 256 ++++---- rinex/src/meteo/record.rs | 66 -- rinex/src/observation/mod.rs | 113 ---- rinex/src/observation/record.rs | 305 +-------- rinex/src/qc/analysis/mod.rs | 29 +- rinex/src/qc/analysis/obs.rs | 610 ++++++++++-------- rinex/src/qc/analysis/sampling.rs | 38 +- rinex/src/qc/analysis/sv.rs | 2 +- rinex/src/qc/mod.rs | 88 ++- rinex/src/qc/opts.rs | 24 +- rinex/src/qc/sampling.rs | 128 ---- rinex/src/tests/mod.rs | 1 - rinex/src/tests/obs.rs | 2 +- rinex/src/tests/stats.rs | 347 ---------- 28 files changed, 804 insertions(+), 1539 deletions(-) delete mode 100644 rinex-cli/config/advanced_repair.json delete mode 100644 rinex-cli/config/advanced_study.json delete mode 100644 rinex-cli/config/basic_manual_gap.json rename rinex-cli/config/{basic.json => gnss_snr30db.json} (100%) create mode 100644 rinex-cli/config/sv_manual_gap.json create mode 100644 rinex/src/algorithm/derivative.rs delete mode 100644 rinex/src/qc/sampling.rs delete mode 100644 rinex/src/tests/stats.rs diff --git a/rinex-cli/config/advanced_repair.json b/rinex-cli/config/advanced_repair.json deleted file mode 100644 index d5d4027be..000000000 --- a/rinex-cli/config/advanced_repair.json +++ /dev/null @@ -1,12 +0,0 @@ -{ - "reporting": "PerConstellation", - "statistical": { - "window": "2 hours" - }, - "processing": { - "cs": "study", - "iono_rate_tolerance": 400E-2, - "iono_rate_tolerance_dt": "1 min", - "clk_drift_window": "10 mins" - } -} diff --git a/rinex-cli/config/advanced_study.json b/rinex-cli/config/advanced_study.json deleted file mode 100644 index d5d4027be..000000000 --- a/rinex-cli/config/advanced_study.json +++ /dev/null @@ -1,12 +0,0 @@ -{ - "reporting": "PerConstellation", - "statistical": { - "window": "2 hours" - }, - "processing": { - "cs": "study", - "iono_rate_tolerance": 400E-2, - "iono_rate_tolerance_dt": "1 min", - "clk_drift_window": "10 mins" - } -} diff --git a/rinex-cli/config/basic_manual_gap.json b/rinex-cli/config/basic_manual_gap.json deleted file mode 100644 index 82618e5d0..000000000 --- a/rinex-cli/config/basic_manual_gap.json +++ /dev/null @@ -1,7 +0,0 @@ -{ - "reporting": "PerConstellation", - "statistical": { - "window": "10%" - }, - "manual_gap": "1O seconds" -} diff --git a/rinex-cli/config/basic.json b/rinex-cli/config/gnss_snr30db.json similarity index 100% rename from rinex-cli/config/basic.json rename to rinex-cli/config/gnss_snr30db.json diff --git a/rinex-cli/config/sv_manual_gap.json b/rinex-cli/config/sv_manual_gap.json new file mode 100644 index 000000000..d7b4f68f6 --- /dev/null +++ b/rinex-cli/config/sv_manual_gap.json @@ -0,0 +1,7 @@ +{ + "classification": "Sv", + "gap_tolerance": { + "centuries": 0, + "nanoseconds": 7200000000000 + } +} diff --git a/rinex-cli/doc/qc.md b/rinex-cli/doc/qc.md index 241c9a22b..dd560ca9f 100644 --- a/rinex-cli/doc/qc.md +++ b/rinex-cli/doc/qc.md @@ -1,10 +1,10 @@ Quality Check (QC) ================== -RINEX quality check is a special mode of this tool, activated with the `--qc` option. +RINEX quality check is a special mode, activated with `--qc`. QC is first developed for Observation files analysis, but this tool -will accept other RINEX files, for which it will compute basic statistical anlysis. +will accept other RINEX files, for which it will compute basic statistical analysis. QC is a well known procedure in data preprocessing. There are a few differences people who are very familiar with `teqc` must @@ -17,7 +17,7 @@ to advanced "teqc" users. Among them: * Quick GNSS filters (`-G`, `-R`, ...) still exist * In augmented mode, `-no_orbit X` is feasible if you know how to operate the [preprocessor](doc/preprocessing.md) -* Similar position anlysis and reporting +* Similar position analysis and reporting * Similar signals analysis and reporting ## Differences with `teqc` @@ -31,13 +31,14 @@ Although we support as many Navigation files to be specified, when the _augmente This allows for example Glonass and other NAV context to be correctly defined. 1. `--fp [FILE]` for the Observation file -2. `--nav [FILE1] [FILE2]..` : pass NAV contexts +2. `--nav [FILE1] --nav [FILE2]..` : define NAV context +3. `--sp3 [FILE1] --sp3 [FILE2]..` : pass SP3 data We will generate products the [workspace](https://github.com/gwbres/rinex/tree/rinex-cli/workspace) folder, that includes QC reports. -Unlike teqc, we do not support BINEX nor SP3 input data/files as of today. +Unlike teqc, we do not support BINEX input data/files as of today. Unlike teqc we do not limit ourselves to the analysis of GPS and Glonass constellations. @@ -57,10 +58,10 @@ For example, this would apply to ## QC specific command line options * `--qc-only`: ensures the tool will only perform the QCs, -other graphs or record analysis is turned off. This is the most efficient -QC mode. +other graphs or record analysis is turned off. +This command line option makes the quality checks and reporting faster. -* `--qc-config`: pass a configuration file for QC reporting and calculations management (see down below) +* `--qc-cfg`: QC analysis and reporting fine tuning (see down below). ## Basic QC (No NAV) @@ -82,83 +83,85 @@ Run this configuration for the most basic QC: ```bash rinex-cli \ - -P gps,glo \ + -P GPS,GLO \ --qc-only \ - --qc-conf rinex-cli/config/basic.json \ + --qc-cfg rinex-cli/config/gnss_snr30db.json \ --fp test_resources/CRNX/V3/ESBC00DNK_R_20201770000_01D_30S_MO.crx.gz ``` -`--qc-conf` is independent to `--qc` activation. If you activate QC without any configuration, the default configuration is used. -basic.json specifies that we want to report on a constellation basis. -If you compare this report to the previous one +## QC Configuration file + +The QC configuration file allows precise control of how analysis are performed. +We provide example configurations in `rinex-cli/config`. + +All parameters are optional, in a configuration file description. That means, +the default value is used if one parameter is omitted. For example, the classification +method is "mandatory" to a QC analysis. If not defined, the GNSS classification method is used. + + +### QC Classification method + +- "GNSS": the provided context is splitted into all constellations that remain +after pre processing. You can one sub report per constellation. +It is most convenient when one or two constellations are left out, for example with `-P GLO,GPS`. +The report is designed in alphabetical order. + +- "Sv" : split by Satellite vehicle. It is most convenient when a few vehicles are isolated, +for example either by retaining a single constellation : `-P GPS`, or by retaining +specific vehicles: `-P G08,G15,C19,C31,R08`. -1. you get one table per constellation. We only retained GPS and Glonass in this example, -therefore the report is made of two tables. -2. All statistical analysis are made on constellations separately and independently -3. A "25%" statistical window is specified. "window" accepts either a Duration -or a percentage of the file. Here we use the latter, and 25% means we will have 4 statistical -analysis performed over the course of 6 hours, because this file is 24h long. -4. You also have less information in this basic configuration, because most calculations are turned off +- "Physics": split by Observable. -Try this configuration now: +### Sv classification case + +The Sv classification case is interesting, because the Sampling/Gap anlalysis emphasizes +when isolated vehicles go out of sight. + +For example, compare the report generated by this "default" run : ```bash -rinex-cli \ - -F gps,glo \ - --qc-only \ - --qc-conf rinex-cli/config/basic_manual_gap.json \ - --fp test_resources/CRNX/V3/ESBC00DNK_R_20201770000_01D_30S_MO.crx.gz +./target/release/rinex-cli \ + --fp test_resources/CRNX/V3/MOJN00DNK_R_20201770000_01D_30S_MO.crx.gz \ + -P G08,G15,G16,R23,C19,C09 \ + --qc ``` -A "10%" statistical slot is specified, you get more statistical analysis because the time slot -now spans approximately 2 hours. -Also "manual_gap" is specified and set to "10 seconds". That means that a duration -of 10 seconds is now considered as an abnormal gap, in the data gap analysis. -In default configuration, there is no manual gap. That means an abnormal gap -is any abnormal duration above the dominant epoch interval (sample rate). +To this one : + +```bash +./target/release/rinex-cli \ + --fp test_resources/CRNX/V3/MOJN00DNK_R_20201770000_01D_30S_MO.crx.gz \ + -P G08,G15,G16,R23,C19,C09 \ + --qc --qc-cfg rinex-cli/config/sv_manual_gap.json +``` -When no configuration files are given, the default configuration is used +### SNR parametrization -```json TODO -``` -## Advanced configurations +### CS analysis -Now let's move on to more "advanced" configuration, in which basically all -calculations are active and customized +TODO -```bash -rinex-cli \ - -P gps,glo \ - --qc-separate \ - --qc-conf rinex-cli/config/basic_manual_gap.json \ - --fp test_resources/CRNX/V3/ESBC00DNK_R_20201770000_01D_30S_MO.crx.gz -``` +## QC with Navigation Data -## Basic QC (with NAV) +Providing Navigation Data is the only way for complete QCs. +With such context, we can now use the "elevation mask" in the configuration file: -Let's go back to our basic demo and provide Navigation context: +Compare these two runs. In the second one, Sv loss of sight is larger because +the go out of sight more rapidly, due to the stringent elevation criteria : ```bash rinex-cli \ -P gps,glo \ - --qc-separate \ - --qc-conf rinex-cli/config/basic.json \ + -P G08,G15,G16,R23,C19,C09 \ + --qc --qc-cfg rinex-cli/config/sv_manual_gap_ev35.json --fp test_resources/CRNX/V3/ESBC00DNK_R_20201770000_01D_30S_MO.crx.gz \ --nav test_resources/OBS/V3/ESBC00DNK_R_20201770000_01D_MN.rnx.gz ``` - -Navigation context is fully taken into account in advanced calculations -```bash -rinex-cli \ - -F gps,glo \ - --qc-separate \ - --qc-conf rinex-cli/config/advanced_study.json \ - --fp test_resources/CRNX/V3/ESBC00DNK_R_20201770000_01D_30S_MO.crx.gz \ - --nav test_resources/OBS/V3/ESBC00DNK_R_20201770000_01D_MN.rnx.gz -`` +## QC on other RINEX +TODO diff --git a/rinex-cli/src/cli.rs b/rinex-cli/src/cli.rs index f6e557d4c..97fe0de78 100644 --- a/rinex-cli/src/cli.rs +++ b/rinex-cli/src/cli.rs @@ -298,6 +298,7 @@ Refer to README")) if let Ok(content) = std::fs::read_to_string(path) { let opts = serde_json::from_str(&content); if let Ok(opts) = opts { + info!("qc parameter file \"{}\"", path); opts } else { error!("failed to parse parameter file \"{}\"", path); diff --git a/rinex-cli/src/main.rs b/rinex-cli/src/main.rs index 5c2bea3fd..671546bea 100644 --- a/rinex-cli/src/main.rs +++ b/rinex-cli/src/main.rs @@ -270,18 +270,18 @@ pub fn main() -> Result<(), rinex::Error> { * Code Multipath analysis */ if cli.multipath() { - let data = ctx - .primary_data() - .observation_phase_align_origin() - .observation_phase_carrier_cycles() - .mp(); - plot::plot_gnss_dcb( - &mut plot_ctx, - "Code Multipath Biases", - "Meters of delay", - &data, - ); - info!("--mp analysis"); + //let data = ctx + // .primary_data() + // .observation_phase_align_origin() + // .observation_phase_carrier_cycles() + // .mp(); + //plot::plot_gnss_dcb( + // &mut plot_ctx, + // "Code Multipath Biases", + // "Meters of delay", + // &data, + //); + warn!("--mp analysis not available yet"); } /* * [GF] recombination visualization requested @@ -471,8 +471,6 @@ pub fn main() -> Result<(), rinex::Error> { */ if qc { info!("entering qc mode"); - let mut qc_opts = cli.qc_config(); - /* * QC Config / versus command line * let the possibility to define some parameters diff --git a/rinex/Cargo.toml b/rinex/Cargo.toml index 93ad091d4..7b3aa21bc 100644 --- a/rinex/Cargo.toml +++ b/rinex/Cargo.toml @@ -15,11 +15,11 @@ rust-version = "1.64" [features] default = [] # no features by default sbas = ["geo", "wkt"] -obs = ["statrs"] +obs = [] meteo = [] nav = [] processing = [] -qc = ["horrorshow", "processing", "sp3"] # rinex Quality Check (mainly OBS RINEX) +qc = ["horrorshow", "processing", "sp3", "statrs"] # rinex Quality Check (mainly OBS RINEX) [package.metadata.docs.rs] all-features = true diff --git a/rinex/src/algorithm/averaging.rs b/rinex/src/algorithm/averaging.rs index eb7c163cb..e59d0ece0 100644 --- a/rinex/src/algorithm/averaging.rs +++ b/rinex/src/algorithm/averaging.rs @@ -1,58 +1,43 @@ -use crate::prelude::*; +use hifitime::{Duration, Epoch}; -pub enum AverageType { - MovingAverage, - Exponential, - Cummulative, +fn moving_average( + data: Vec<(Epoch, T)>, + window: Duration, +) -> Vec<(Epoch, T)> { + let mut acc = T::default(); + let mut prev_epoch: Option = None; + let mut ret: Vec<(Epoch, T)> = Vec::new(); + for (epoch, value) in data {} + ret } -pub struct Averager { - buffer: Vec, - next_epoch: Option, - window: Duration, - avgtype: AverageType, +#[derive(Debug, Clone, Copy)] +pub enum Averager { + MovingAverage(Duration), } -impl Averager { - pub fn new(window: Duration) -> Self { - Self { - buffer: Vec::new(), - next_epoch: None, - window, - } +impl Default for Averager { + fn default() -> Self { + Self::MovingAverage(Duration::from_seconds(600.0_f64)) } - pub fn average(&mut self, data: (f64, Epoch)) -> Option { - match self.avgtype { - AverageType::MovingAverage => self.moving_average(data), - AverageType::Exponential => self.moving_average(data), - AverageType::Cummulative => self.moving_average(data), - } +} + +impl Averager { + pub fn mov(window: Duration) -> Self { + Self::MovingAverage(window) } - pub fn moving_average(&mut self, data: (f64, Epoch)) -> Option { - self.buffer.push(data.0); - if self.next_epoch.is_none() { - self.next_epoch = Some(data.1 + self.window); - } - if let Some(next_epoch) = self.next_epoch { - if data.1 >= next_epoch { - self.next_epoch = Some(data.1 + self.window); - let mut avg = 0.0_f64; - for b in &self.buffer { - avg += b; - } - let ret = avg / self.buffer.len() as f64; - self.buffer.clear(); - return Some(ret); - } + pub fn eval(&self, input: Vec<(Epoch, T)>) -> Vec<(Epoch, T)> { + match self { + Self::MovingAverage(dt) => moving_average(input, *dt), } - None } } #[cfg(test)] mod test { + use super::*; #[test] fn test_moving_average() { - + let mov = Averager::mov(Duration::from_seconds(10.0_f64)); } } diff --git a/rinex/src/algorithm/derivative.rs b/rinex/src/algorithm/derivative.rs new file mode 100644 index 000000000..5de98965d --- /dev/null +++ b/rinex/src/algorithm/derivative.rs @@ -0,0 +1,38 @@ +use hifitime::Epoch; + +pub struct Derivative { + order: usize, +} + +/* + * Derivative of an number of array sorted by chronological Epoch + */ +pub(crate) fn derivative(input: Vec<(Epoch, f64)>, buf: &mut Vec<(Epoch, f64)>) { + let mut prev: Option<(Epoch, f64)> = None; + for (e, value) in input { + if let Some((prev_e, prev_v)) = prev { + let dt = e - prev_e; + let dy = (value - prev_v) / dt.to_seconds(); + buf.push((e, dy)); + } + prev = Some((e, value)); + } +} + +/* + * Derivative^2 of an number of array sorted by chronological Epoch + */ +impl Derivative { + pub fn new(order: usize) -> Self { + Self { order } + } + pub fn eval(&self, input: Vec<(Epoch, f64)>) -> Vec<(Epoch, f64)> { + let mut buf: Vec<(Epoch, f64)> = Vec::with_capacity(input.len()); + derivative(input, &mut buf); + //for i in 1..self.order { + // derivative(&ret, &mut ret); + //} + //ret + buf + } +} diff --git a/rinex/src/algorithm/mod.rs b/rinex/src/algorithm/mod.rs index 1b77a046e..5ee754691 100644 --- a/rinex/src/algorithm/mod.rs +++ b/rinex/src/algorithm/mod.rs @@ -1,3 +1,5 @@ +//mod averaging; +mod derivative; mod filters; mod target; @@ -7,3 +9,6 @@ pub use filters::{ Decimate, DecimationFilter, DecimationType, Filter, InterpFilter, InterpMethod, Interpolate, Mask, MaskFilter, MaskOperand, Preprocessing, Smooth, SmoothingFilter, SmoothingType, }; + +//pub use averaging::Averager; +pub use derivative::Derivative; diff --git a/rinex/src/hardware.rs b/rinex/src/hardware.rs index 78ed15509..cc57709ab 100644 --- a/rinex/src/hardware.rs +++ b/rinex/src/hardware.rs @@ -106,7 +106,7 @@ impl HtmlReport for Antenna { table(class="table is-bordered") { tr { th { - : "Antenna model" + : "Model" } th { : "SN#" @@ -167,7 +167,7 @@ impl HtmlReport for Rcvr { } fn to_inline_html(&self) -> Box { box_html! { - table(class="table is-bordered") { + table(class="table is-bordered; style=\"margin-bottom: 20px\"") { tr { th { : "Model" diff --git a/rinex/src/header.rs b/rinex/src/header.rs index afc574841..7a046b3a6 100644 --- a/rinex/src/header.rs +++ b/rinex/src/header.rs @@ -1833,14 +1833,32 @@ impl HtmlReport for Header { } fn to_inline_html(&self) -> Box { box_html! { - @ if let Some(antenna) = &self.rcvr_antenna { - div(id="antenna") { - : antenna.to_inline_html() + tr { + th { + : "Antenna" + } + @ if let Some(antenna) = &self.rcvr_antenna { + td { + : antenna.to_inline_html() + } + } else { + td { + : "No information" + } } } - @ if let Some(rcvr) = &self.rcvr { - div(id="receiver") { - : rcvr.to_inline_html() + tr { + th { + : "Receiver" + } + @ if let Some(rcvr) = &self.rcvr { + td { + : rcvr.to_inline_html() + } + } else { + td { + : "No information" + } } } } diff --git a/rinex/src/lib.rs b/rinex/src/lib.rs index 59867e605..701a10d99 100644 --- a/rinex/src/lib.rs +++ b/rinex/src/lib.rs @@ -984,36 +984,6 @@ impl Rinex { results } */ - /* - /// Extracts Raw Carrier Phase observations, - /// from this Observation record, on an epoch basis an per space vehicle. - /// Does not produce anything if self is not an Observation RINEX. - pub fn observation_phase(&self) -> BTreeMap<(Epoch, EpochFlag), HashMap>> { - let mut ret: BTreeMap>> = BTreeMap::new(); - if let Some(r) = self.record.as_obs() { - for ((e, _), (_, sv)) in record.iter() { - let mut map: BTreeMap> = BTreeMap::new(); - for (sv, obs) in sv.iter() { - let mut v: Vec<(String, f64)> = Vec::new(); - for (observable, data) in obs.iter() { - if observable.is_phase_observation() { - v.push((code.clone(), data.obs)); - } - } - if v.len() > 0 { - // did come with at least 1 Phase obs - map.insert(*sv, v); - } - } - if map.len() > 0 { - // did produce something - results.insert(*e, map); - } - } - } - ret - } - */ /* /// Extracts Carrier phases without Ionospheric path delay contributions, /// by extracting [carrier_phases] and using the differential (dual frequency) compensation. @@ -1802,14 +1772,52 @@ impl Rinex { } } +#[cfg(feature = "obs")] +use crate::observation::Snr; + /* * OBS RINEX specific methods: only available on crate feature. * Either specific Iterators, or meaningful data we can extract. */ - #[cfg(feature = "obs")] #[cfg_attr(docrs, doc(cfg(feature = "obs")))] impl Rinex { + /// Returns a Unique Iterator over identified [`Carrier`]s + pub fn carrier(&self) -> Box + '_> { + Box::new(self.observation().flat_map(|(_, (_, sv))| { + sv.iter().flat_map(|(sv, observations)| { + observations + .keys() + .filter_map(|observable| { + if let Ok(carrier) = observable.carrier(sv.constellation) { + Some(carrier) + } else { + None + } + }) + .fold(vec![], |mut list, item| { + if !list.contains(&item) { + list.push(item); + } + list + }) + .into_iter() + }) + })) + } + pub fn code(&self) -> Box + '_> { + Box::new( + self.observation() + .flat_map(|(_, (_, sv))| { + sv.iter().flat_map(|(_, observations)| { + observations + .keys() + .filter_map(|observable| observable.code()) + }) + }) + .unique(), + ) + } /// Returns ([`Epoch`] [`EpochFlag`]) iterator, where each {`EpochFlag`] /// validates or invalidates related [`Epoch`] /// ``` @@ -2023,6 +2031,21 @@ impl Rinex { }) })) } + /// Returns an Iterator over signal SNR indications. + /// All observation that did not come with such indication are filtered out. + pub fn snr(&self) -> Box + '_> { + Box::new(self.observation().flat_map(|(e, (_, vehicles))| { + vehicles.iter().flat_map(|(sv, observations)| { + observations.iter().filter_map(|(obs, obsdata)| { + if let Some(snr) = obsdata.snr { + Some((*e, *sv, obs, snr)) + } else { + None + } + }) + }) + })) + } } #[cfg(feature = "nav")] @@ -2815,82 +2838,6 @@ impl Decimate for Rinex { } } -#[cfg(feature = "obs")] -use crate::observation::Observation; - -#[cfg(feature = "obs")] -#[cfg_attr(docrs, doc(cfg(feature = "obs")))] -impl Observation for Rinex { - fn min(&self) -> (Option, HashMap>) { - if let Some(r) = self.record.as_obs() { - r.min() - } else { - (None, HashMap::new()) - } - } - fn min_observable(&self) -> HashMap { - if let Some(r) = self.record.as_obs() { - r.min_observable() - } else if let Some(r) = self.record.as_meteo() { - r.min_observable() - } else { - HashMap::new() - } - } - fn max(&self) -> (Option, HashMap>) { - if let Some(r) = self.record.as_obs() { - r.max() - } else { - (None, HashMap::new()) - } - } - fn max_observable(&self) -> HashMap { - if let Some(r) = self.record.as_obs() { - r.max_observable() - } else if let Some(r) = self.record.as_meteo() { - r.max_observable() - } else { - HashMap::new() - } - } - fn mean(&self) -> (Option, HashMap>) { - if let Some(r) = self.record.as_obs() { - r.mean() - } else { - (None, HashMap::new()) - } - } - fn mean_observable(&self) -> HashMap { - if let Some(r) = self.record.as_obs() { - r.mean_observable() - } else if let Some(r) = self.record.as_meteo() { - r.mean_observable() - } else { - HashMap::new() - } - } - fn std_dev(&self) -> (Option, HashMap>) { - if let Some(r) = self.record.as_obs() { - r.std_dev() - } else { - (None, HashMap::new()) - } - } - fn std_dev_observable(&self) -> HashMap { - HashMap::new() - } - fn std_var(&self) -> (Option, HashMap>) { - if let Some(r) = self.record.as_obs() { - r.std_var() - } else { - (None, HashMap::new()) - } - } - fn std_var_observable(&self) -> HashMap { - HashMap::new() - } -} - #[cfg(feature = "obs")] use observation::Dcb; @@ -2905,19 +2852,88 @@ impl Dcb for Rinex { } } -#[cfg(feature = "obs")] -use observation::Mp; - -#[cfg(feature = "obs")] -impl Mp for Rinex { - fn mp(&self) -> HashMap>> { - if let Some(r) = self.record.as_obs() { - r.dcb() - } else { - panic!("wrong rinex type"); - } - } -} +//#[cfg(feature = "obs")] +//use observation::Mp; +// +//#[cfg(feature = "obs")] +//impl Mp for Rinex { +// fn mp(&self) -> HashMap>> { +// /* +// * Determine mean value of all observed Phase and Pseudo Range +// * observations, for all Sv +// */ +// let mut mean: HashMap> = HashMap::new(); +// /* +// * Associate a Phase code to all PR codes +// */ +// let mut associations: HashMap = HashMap::new(); +// +// let pr_codes: Vec = self.observable() +// .filter_map(|obs| +// if obs.is_pseudorange_observable() { +// Some(obs.clone()) +// } else { +// None +// }) +// .collect(); +// +// for observable in self.observable() { +// if !observable.is_phase_observable() { +// if !observable.is_pseudorange_observable() { +// continue; +// } +// } +// // code associations (for future combianations) +// if observable.is_phase_observable() { +// for pr_code in &pr_codes { +// +// } +// } +// +// for sv in self.sv() { +// /* +// * average phase values +// */ +// let phases: Vec = self.carrier_phase() +// .filter_map(|(_, svnn, obs, ph)| { +// if sv == svnn && observable == obs { +// Some(ph) +// } else { +// None +// } +// }) +// .collect(); +// if let Some(data) = mean.get_mut(&sv) { +// data.insert(observable.clone(), phases.mean()); +// } else { +// let mut map: HashMap = HashMap::new(); +// map.insert(observable.clone(), phases.mean()); +// mean.insert(sv, map); +// } +// /* +// * average PR values +// */ +// let pr: Vec = self.pseudo_range() +// .filter_map(|(_, svnn, obs, pr)| { +// if sv == svnn && observable == obs { +// Some(pr) +// } else { +// None +// } +// }) +// .collect(); +// if let Some(data) = mean.get_mut(&sv) { +// data.insert(observable.clone(), pr.mean()); +// } else { +// let mut map: HashMap = HashMap::new(); +// map.insert(observable.clone(), pr.mean()); +// mean.insert(sv, map); +// } +// } +// } +// HashMap::new() +// } +//} #[cfg(feature = "obs")] use observation::Combine; diff --git a/rinex/src/meteo/record.rs b/rinex/src/meteo/record.rs index 2df3cd840..7e15a611f 100644 --- a/rinex/src/meteo/record.rs +++ b/rinex/src/meteo/record.rs @@ -436,69 +436,3 @@ impl Interpolate for Record { unimplemented!("meteo:record:interpolate_mut()") } } - -#[cfg(feature = "obs")] -use crate::observation::{Observation, StatisticalOps}; - -#[cfg(feature = "obs")] -use statrs::statistics::Statistics; - -#[cfg(feature = "obs")] -fn statistical_estimate(rec: &Record, ops: StatisticalOps) -> HashMap { - let mut ret: HashMap = HashMap::new(); - let mut dataset: HashMap> = HashMap::new(); - for (_t, observables) in rec { - for (observable, point) in observables { - if let Some(data) = dataset.get_mut(&observable) { - data.push(*point); - } else { - dataset.insert(observable.clone(), vec![*point]); - } - } - } - for (observable, data) in dataset { - let stats = match ops { - StatisticalOps::Min => data.min(), - StatisticalOps::Max => data.max(), - StatisticalOps::Mean => data.mean(), - StatisticalOps::StdDev => data.std_dev(), - StatisticalOps::StdVar => data.variance(), - }; - ret.insert(observable.clone(), stats); - } - ret -} - -#[cfg(feature = "obs")] -impl Observation for Record { - fn min(&self) -> (Option, HashMap>) { - unimplemented!("only on OBS RINEX!"); - } - fn max(&self) -> (Option, HashMap>) { - unimplemented!("only on OBS RINEX!"); - } - fn mean(&self) -> (Option, HashMap>) { - unimplemented!("only on OBS RINEX!"); - } - fn std_var(&self) -> (Option, HashMap>) { - unimplemented!("only on OBS RINEX!"); - } - fn std_dev(&self) -> (Option, HashMap>) { - unimplemented!("only on OBS RINEX!"); - } - fn min_observable(&self) -> HashMap { - statistical_estimate(self, StatisticalOps::Min) - } - fn max_observable(&self) -> HashMap { - statistical_estimate(self, StatisticalOps::Max) - } - fn std_dev_observable(&self) -> HashMap { - statistical_estimate(self, StatisticalOps::StdDev) - } - fn std_var_observable(&self) -> HashMap { - statistical_estimate(self, StatisticalOps::StdVar) - } - fn mean_observable(&self) -> HashMap { - statistical_estimate(self, StatisticalOps::Mean) - } -} diff --git a/rinex/src/observation/mod.rs b/rinex/src/observation/mod.rs index d311885a1..3725d2acb 100644 --- a/rinex/src/observation/mod.rs +++ b/rinex/src/observation/mod.rs @@ -174,122 +174,9 @@ impl HeaderFields { } } -#[cfg(feature = "obs")] -#[derive(Debug, Clone, Copy)] -pub(crate) enum StatisticalOps { - Min, - Max, - Mean, - StdDev, - StdVar, -} - #[cfg(feature = "obs")] use std::collections::BTreeMap; -/// OBS RINEX specific analysis trait. -/// Include this trait to unlock Observation analysis, mainly statistical analysis. -#[cfg(feature = "obs")] -#[cfg_attr(docrs, doc(cfg(feature = "obs")))] -pub trait Observation { - /// Returns minimum value observed, throughout all epochs, sorted by Observable. - /// This also applies to clock receiver estimate, - /// when requested on OBS RINEX files, not METEO files. - /// ``` - /// use rinex::*; // prelude + macros - /// use rinex::prelude::*; - /// use std::str::FromStr; // observable! - /// use rinex::observation::Observation; // .min_observable() - /// - /// // OBS RINEX example - /// let rinex = Rinex::from_file("../test_resources/OBS/V3/DUTH0630.22O") - /// .unwrap(); - /// let min_values = rinex.min_observable(); - /// for (observable, min_value) in min_values { - /// if observable == observable!("S1C") { - /// // minimum signal strength for carrier 1 - /// assert_eq!(min_value, 37.75); // L1 carrier min (worst) RSSI - /// } - /// } - /// - /// // METEO RINEX example - /// let rinex = Rinex::from_file("../test_resources/MET/V2/clar0020.00m") - /// .unwrap(); - /// let min_values = rinex.min_observable(); - /// for (observable, min_value) in min_values { - /// if observable == Observable::Temperature { - /// assert_eq!(min_value, 8.4); // min value encountered on that day - /// } - /// } - /// ``` - fn min_observable(&self) -> HashMap; - - /// Returns maximal value observed, throughout all epochs, sorted by Observable. - /// See [Self::min_observable()] for API example. - fn max_observable(&self) -> HashMap; - - /// Returns mean observation, throughout all epochs, sorted by Observable. - /// See [Self::min_observable()] for API example. - fn mean_observable(&self) -> HashMap; - - /// Returns standard deviation for all Observables. - /// See [Self::min_observable()] for API example. - fn std_dev_observable(&self) -> HashMap; - - /// Returns standard variance for all Observables. - /// See [Self::min_observable()] for API example. - fn std_var_observable(&self) -> HashMap; - - /// Returns minimum value observed throughout all epochs sorted by - /// Satellite vehicle and Observable. This does not apply to METEO - /// RINEX files. - /// ``` - /// use rinex::*; - /// use rinex::prelude::*; // basics - /// use std::str::FromStr; // sv!, observable! - /// use rinex::observation::Observation; // .min() - /// - /// // OBS RINEX example - /// let rinex = Rinex::from_file("../test_resources/OBS/V3/DUTH0630.22O") - /// .unwrap(); - /// - /// let (min_clock, min_values) = rinex.min(); - /// assert!(min_clock.is_none()); // we don't have an example file with such information yet - /// - /// for (sv, observables) in min_values { - /// if sv == sv!("G08") { - /// for (observable, min_value) in observables { - /// if observable == observable!("S1C") { - /// // minimum signal strength for carrier 1 - /// // for that particular vehicle - /// } - /// } - /// } - /// } - /// ``` - fn min(&self) -> (Option, HashMap>); - - /// Returns maximal value observed throughout all epochs sorted by - /// Satellite vehicle and Observable. This does not apply to METEO - /// RINEX files. See [Self::min()] for API example. - fn max(&self) -> (Option, HashMap>); - - /// Returns mean value observed throughout all epochs sorted by - /// Satellite vehicle and Observable. This does not apply to METEO - /// RINEX files. See [Self::min()] for API example. - fn mean(&self) -> (Option, HashMap>); - - /// Returns observations deviation throughout all epochs sorted by - /// Satellite vehicle and Observable. This does not apply to METEO - /// RINEX files. See [Self::min()] for API example. - fn std_dev(&self) -> (Option, HashMap>); - - /// Returns observations variance throughout all epochs sorted by - /// Satellite vehicle and Observable. This does not apply to METEO - /// RINEX files. See [Self::min()] for API example. - fn std_var(&self) -> (Option, HashMap>); -} - /// GNSS signal recombination trait. /// Import this to recombine OBS RINEX with usual recombination methods. /// This only applies to OBS RINEX records. diff --git a/rinex/src/observation/record.rs b/rinex/src/observation/record.rs index 09c200484..4e7c88aab 100644 --- a/rinex/src/observation/record.rs +++ b/rinex/src/observation/record.rs @@ -1282,167 +1282,6 @@ impl Decimate for Record { } } -#[cfg(feature = "obs")] -use statrs::statistics::Statistics; - -#[cfg(feature = "obs")] -use crate::observation::{Observation, StatisticalOps}; - -#[cfg(feature = "obs")] -/* - * Evaluate a specific statistical estimate on this record data set - */ -fn statistical_estimate( - rec: &Record, - ops: StatisticalOps, -) -> (Option, HashMap>) { - let mut ret: (Option, HashMap>) = (None, HashMap::new()); - - // Vectorize clock offset (in time), so we can invoke statrs() on it - let data: Vec = rec - .iter() - .filter_map(|(_, (clk, _))| { - if let Some(clk) = clk { - Some(*clk) - } else { - None - } - }) - .collect(); - - if data.len() > 0 { - // data exists - match ops { - StatisticalOps::Min => ret.0 = Some(data.min()), - StatisticalOps::Max => ret.0 = Some(data.max()), - StatisticalOps::Mean => ret.0 = Some(data.mean()), - StatisticalOps::StdDev => ret.0 = Some(data.std_dev()), - StatisticalOps::StdVar => ret.0 = Some(data.variance()), - } - } - - // Vectorize data (in time) so we can invoke statrs() on it - let mut data: HashMap>> = HashMap::new(); - - for (_epoch, (_clk, svnn)) in rec { - for (sv, observations) in svnn { - for (observable, obs_data) in observations { - if let Some(data) = data.get_mut(&sv) { - if let Some(data) = data.get_mut(&observable) { - data.push(obs_data.obs); // append - } else { - data.insert(observable.clone(), vec![obs_data.obs]); - } - } else { - let mut map: HashMap> = HashMap::new(); - map.insert(observable.clone(), vec![obs_data.obs]); - data.insert(sv.clone(), map); - } - } - } - } - - // build returned data - for (sv, observables) in data { - for (observable, data) in observables { - // invoke statrs - let stats = match ops { - StatisticalOps::Min => data.min(), - StatisticalOps::Max => data.max(), - StatisticalOps::Mean => data.mean(), - StatisticalOps::StdDev => data.std_dev(), - StatisticalOps::StdVar => data.variance(), - }; - if let Some(data) = ret.1.get_mut(&sv) { - data.insert(observable.clone(), stats); - } else { - let mut map: HashMap = HashMap::new(); - map.insert(observable.clone(), stats); - ret.1.insert(sv.clone(), map); - } - } - } - - ret -} - -#[cfg(feature = "obs")] -/* - * Evaluate a specific statistical estimate on this record data set - */ -fn statistical_observable_estimate(rec: &Record, ops: StatisticalOps) -> HashMap { - let mut ret: HashMap = HashMap::new(); - // For all statistical estimates but StdDev/StdVar - // We can compute across all vehicles then shrink 1D - match ops { - StatisticalOps::StdVar | StatisticalOps::StdDev => {}, - _ => { - // compute across all observables, then per Sv - let (_, data) = statistical_estimate(rec, ops); - let mut data_set: HashMap> = HashMap::new(); - for (_sv, observables) in data { - for (observable, value) in observables { - if let Some(data) = data_set.get_mut(&observable) { - data.push(value); - } else { - data_set.insert(observable.clone(), vec![value]); - } - } - } - for (observable, values) in data_set { - // invoke statrs - match ops { - StatisticalOps::Min => { - ret.insert(observable.clone(), values.min()); - }, - StatisticalOps::Max => { - ret.insert(observable.clone(), values.max()); - }, - StatisticalOps::Mean => { - ret.insert(observable.clone(), values.mean()); - }, - _ => {}, // does not apply - }; - } - }, - } - ret -} - -#[cfg(feature = "obs")] -impl Observation for Record { - fn min(&self) -> (Option, HashMap>) { - statistical_estimate(self, StatisticalOps::Min) - } - fn max(&self) -> (Option, HashMap>) { - statistical_estimate(self, StatisticalOps::Max) - } - fn mean(&self) -> (Option, HashMap>) { - statistical_estimate(self, StatisticalOps::Mean) - } - fn std_var(&self) -> (Option, HashMap>) { - statistical_estimate(self, StatisticalOps::StdVar) - } - fn std_dev(&self) -> (Option, HashMap>) { - statistical_estimate(self, StatisticalOps::StdDev) - } - fn min_observable(&self) -> HashMap { - statistical_observable_estimate(self, StatisticalOps::Min) - } - fn max_observable(&self) -> HashMap { - statistical_observable_estimate(self, StatisticalOps::Max) - } - fn std_dev_observable(&self) -> HashMap { - statistical_observable_estimate(self, StatisticalOps::StdDev) - } - fn std_var_observable(&self) -> HashMap { - statistical_observable_estimate(self, StatisticalOps::StdVar) - } - fn mean_observable(&self) -> HashMap { - statistical_observable_estimate(self, StatisticalOps::Mean) - } -} - #[cfg(feature = "obs")] use crate::observation::Combine; @@ -1834,7 +1673,7 @@ impl Combine for Record { #[cfg(feature = "obs")] use crate::{ carrier, - observation::{Dcb, Mp}, + observation::Dcb, //Mp}, }; #[cfg(feature = "obs")] @@ -2000,148 +1839,6 @@ impl IonoDelay for Record { } } -#[cfg(feature = "obs")] -impl Mp for Record { - fn mp(&self) -> HashMap>> { - let mut ret: HashMap>> = - HashMap::new(); - /* - * Determine mean value of all datasets - */ - let (_, mean) = self.mean(); // drop mean{clk} - //println!("MEAN VALUES {:?}", mean); //DEBUG - /* - * Run algorithm - */ - let mut associated: HashMap = HashMap::new(); // Ph code to associate to this Mpx - // for operation consistency - for (epoch, (_, vehicles)) in self { - for (sv, observations) in vehicles { - let _mean_sv = mean.get(&sv).unwrap(); - for (lhs_observable, lhs_data) in observations { - if lhs_observable.is_pseudorange_observable() { - let pr_i = lhs_data.obs; // - mean_sv.get(lhs_code).unwrap().1; - let lhs_code = lhs_observable.to_string(); - let mp_code = &lhs_code[2..]; //TODO will not work on RINEX2 - let lhs_carrier = &lhs_code[1..2]; - let mut ph_i: Option = None; - let mut ph_j: Option = None; - /* - * This will restrict combinations to - * 1 => 2 - * and M => 1 - */ - let rhs_carrier = match lhs_carrier { - "1" => "2", - _ => "1", - }; - /* - * locate related L_i PH code - */ - for (observable, data) in observations { - let ph_code = format!("L{}", mp_code); - let code = observable.to_string(); - if code.eq(&ph_code) { - ph_i = Some(data.obs); // - mean_sv.get(code).unwrap().1); - break; // DONE - } - } - /* - * locate another L_j PH code - */ - if let Some(to_locate) = associated.get(mp_code) { - /* - * We already have an association, keep it consistent throughout - * operations - */ - for (observable, data) in observations { - let code = observable.to_string(); - let carrier_code = &code[1..2]; - if carrier_code == rhs_carrier { - // correct carrier signal - if code.eq(to_locate) { - // match - ph_j = Some(data.obs); // - mean_sv.get(code).unwrap().1); - break; // DONE - } - } - } - } else { - // first: prefer the same code against rhs carrier - let to_locate = format!("L{}{}", rhs_carrier, &mp_code[1..]); - for (observable, data) in observations { - let code = observable.to_string(); - let carrier_code = &code[1..2]; - if carrier_code == rhs_carrier { - // correct carrier - if code.eq(&to_locate) { - // match - ph_j = Some(data.obs); // - mean_sv.get(code).unwrap().1); - associated.insert(mp_code.to_string(), code.clone()); - break; // DONE - } - } - } - if ph_j.is_none() { - /* - * Same code against different carrier does not exist - * try to grab another PH code, against rhs carrier - */ - for (observable, data) in observations { - let code = observable.to_string(); - let carrier_code = &code[1..2]; - if carrier_code == rhs_carrier { - if observable.is_phase_observable() { - ph_j = Some(data.obs); // - mean_sv.get(code).unwrap().1); - associated.insert(mp_code.to_string(), code.clone()); - break; // DONE - } - } - } - } - } - if ph_i.is_none() || ph_j.is_none() { - break; // incomplete associations, do not proceed further - } - let ph_i = ph_i.unwrap(); - let ph_j = ph_j.unwrap(); - let lhs_carrier = lhs_observable.carrier(sv.constellation).unwrap(); - let rhs_carrier = lhs_observable //rhs_observable TODO - .carrier(sv.constellation) - .unwrap(); - /*let gamma = (lhs_carrier.frequency() / rhs_carrier.frequency()).powf(2.0); - let alpha = (gamma +1.0_f64) / (gamma - 1.0_f64); - let beta = 2.0_f64 / (gamma - 1.0_f64); - let mp = pr_i - alpha * ph_i + beta * ph_j;*/ - - let alpha = 2.0_f64 * rhs_carrier.frequency().powf(2.0) - / (lhs_carrier.frequency().powf(2.0) - - rhs_carrier.frequency().powf(2.0)); - let mp = pr_i - ph_i - alpha * (ph_i - ph_j); - if let Some(data) = ret.get_mut(mp_code) { - if let Some(data) = data.get_mut(&sv) { - data.insert(*epoch, mp); - } else { - let mut bmap: BTreeMap<(Epoch, EpochFlag), f64> = BTreeMap::new(); - bmap.insert(*epoch, mp); - data.insert(*sv, bmap); - } - } else { - let mut bmap: BTreeMap<(Epoch, EpochFlag), f64> = BTreeMap::new(); - let mut map: BTreeMap> = - BTreeMap::new(); - bmap.insert(*epoch, mp); - map.insert(*sv, bmap); - ret.insert(mp_code.to_string(), map); - } - } - } - } - } - ret - } -} - #[cfg(test)] mod test { use super::*; diff --git a/rinex/src/qc/analysis/mod.rs b/rinex/src/qc/analysis/mod.rs index 504629ad8..3d1b862be 100644 --- a/rinex/src/qc/analysis/mod.rs +++ b/rinex/src/qc/analysis/mod.rs @@ -1,6 +1,6 @@ use super::{pretty_array, HtmlReport, QcOpts}; use crate::prelude::*; -use horrorshow::{helper::doctype, RenderBox}; +use horrorshow::{helper::doctype, RenderBox}; //table_lengthy_td mod sv; @@ -66,32 +66,33 @@ impl HtmlReport for QcAnalysis { } fn to_inline_html(&self) -> Box { box_html! { - div(id="analysis") { + div(id="analysis; style=\"display: flex; flex-direction: column; gap: 30px\"") { div(id="sampling") { - h3(class="title") { - : "Sampling" - } - table(class="table is-bordered") { + table(class="table is-bordered; style=\"margin-bottom: 30px\"") { + thead { + th { + : "Sampling" + } + } tbody { : self.sampling.to_inline_html() } } } div(id="sv") { - h3(class="title") { - : "Sv" - } - table(class="table is-bordered") { + table(class="table is-bordered; style=\"margin-bottom: 30px\"") { + thead { + th { + : "Sv" + } + } tbody { : self.sv.to_inline_html() } } } div(id="observations") { - h3(class="title") { - : "Observations" - } - table(class="table is-bordered") { + table(class="table is-bordered; style=\"margin-bottom: 30px\"") { tbody { : self.observ.to_inline_html() } diff --git a/rinex/src/qc/analysis/obs.rs b/rinex/src/qc/analysis/obs.rs index b4255842b..c5d5f3cbd 100644 --- a/rinex/src/qc/analysis/obs.rs +++ b/rinex/src/qc/analysis/obs.rs @@ -1,15 +1,16 @@ use crate::{carrier, observation::Snr, prelude::*, Carrier}; +use itertools::Itertools; +use std::str::FromStr; use super::{pretty_array, QcOpts}; use std::collections::HashMap; +use statrs::statistics::Statistics; + #[cfg(feature = "processing")] use crate::preprocessing::*; // include preprocessing toolkit when feasible, // for complete analysis -#[cfg(feature = "obs")] -use crate::observation::Observation; // having this feature unlocks full OBS RINEX analysis - /* * GNSS signal special formatting */ @@ -31,106 +32,135 @@ fn report_signals(list: &Vec) -> String { } /* - * Epoch anomalies formatter + * Report RX Clock drift analysis */ -fn report_anomalies(anomalies: &Vec<(Epoch, EpochFlag)>) -> Box { - if anomalies.len() == 0 { - box_html! { +fn report_clock_drift(data: &Vec<(Epoch, f64)>) -> Box { + box_html! { + @ if data.is_empty() { table(class="table is-bordered") { - th { - : "Anomalies" - } - td { - : "None" + tr { + th { + : "Unfeasible" + } + td { + : "Missing Data" + } } } - } - } else { - box_html! { + } else { table(class="table is-bordered") { - thead { + tr { th { - : "Anomalies" + : "Epoch" } th { - : "Power failure" + : "Mean Clock drift [s/s]" } - th { - : "Antenna movement detected" + } + @ for (epoch, drift) in data { + tr { + td { + : epoch.to_string() + } + td { + : format!("{:e}", drift) + } } + } + } + } + } +} + +/* + * Epoch anomalies formatter + */ +fn report_anomalies<'a>( + cs: &'a Vec, + power: &'a Vec, + other: &'a Vec<(Epoch, EpochFlag)>, +) -> Box { + box_html! { + tr { + th { + : "Power Failures" + } + @ if power.is_empty() { + td { + : "None" + } + } else { + td { + : format!("{:?}", power) + } + tr { th { - : "Kinematic" + : "Longest" } - th { - : "External event" + td { + //: power.iter().max_by(|(_, d1), (_, d2)| d1.cmp(d2)).unwrap().to_string() + : "TODO" } - th { - : "Cycle Slips" + td { + : "Average Duration" + } + td { + : "TODO" } } - tbody { - @ for (epoch, _flag) in anomalies { - tr { - td { - : epoch.to_string() - } - /*@match flag { - EpochFlag::PowerFailure => { - td { - : "x" - } - }, - EpochFlag::AntennaBeingMoved => { - td { - : "" - } - td { - : "x" - } - }, - EpochFlag::NewSiteOccupation => { - td { - : "" - } - td { - : "" - } - td { - : "x" - } - }, - EpochFlag::ExternalEvent => { - td { - : "" - } - td { - : "" - } - td { - : "" - } - td { - : "x" - } - }, - EpochFlag::CycleSlip => { - td { - : "" - } - td { - : "" - } - td { - : "" - } - td { - : "" - } - td { - : "x" - } - }, - }*/ + } + } + tr { + th { + : "Cycle slip(s)" + } + @ if cs.is_empty() { + td { + : "None" + } + } else { + td { + : format!("{:?}", cs) + } + } + } + tr { + th { + : "Other anomalies" + } + @ if other.is_empty() { + td { + : "None" + } + } else { + td { + : "Epoch" + } + td { + : "Event" + } + @ for (e, event) in other { + td { + : "" + } + td { + : format!("{}", e) + } + @ if *event == EpochFlag::AntennaBeingMoved { + td { + : "Antenna Being Moved" + } + } else if *event == EpochFlag::NewSiteOccupation { + td { + : "New Site Occupation" + } + } else if *event == EpochFlag::ExternalEvent { + td { + : "External Event" + } + } else { + td { + : "Other" } } } @@ -152,10 +182,9 @@ fn report_epoch_completion( ) -> Box { box_html! { table(class="table is-bordered") { - tbody { tr { th { - : "Total Epochs" + : "Total#" } td { : total.to_string() @@ -170,8 +199,19 @@ fn report_epoch_completion( } } tr { - th { - : "Complete" + td { + b { + : "Complete" + } + p { + : "Epochs with at least Phase + PR" + } + p { + : "in dual frequency, with" + } + p { + : "both SNR and elev above masks" + } } @ for (signal, count) in complete { td { @@ -179,12 +219,63 @@ fn report_epoch_completion( : format!("L1/{}", signal) } p { - : count.to_string() - } - b { - : format!("{}%", count * 100/total) + : format!("{} ({}%)", count, count * 100 / total) } - //: format!("L1/{}: {} ({}%)", signal, count, count *100/total) + } + } + } + } + } +} + +/* + * SNR analysis report + */ +fn report_snr_statistics( + snr_stats: &HashMap, +) -> Box { + box_html! { + tr { + td { + : "" + } + @ for (observable, _) in snr_stats { + @ if observable.is_phase_observable() || observable.is_pseudorange_observable() { + td { + : observable.to_string() + } + } + } + } + tr { + th { + : "Best" + } + @ for (observable, (_, max)) in snr_stats { + @ if observable.is_phase_observable() || observable.is_pseudorange_observable() { + td { + b { + : format!("{:e}", Snr::from_str(&format!("{}", max.1)).unwrap()) + } + p { + : format!("@{}", max.0) + } + } + } + } + } + tr { + th { + : "Worst" + } + @ for (observable, (min, _)) in snr_stats { + @ if observable.is_phase_observable() || observable.is_pseudorange_observable() { + td { + b { + : format!("{:e}", Snr::from_str(&format!("{}", min.1)).unwrap()) + } + p { + : format!("@{}", min.0) } } } @@ -196,9 +287,7 @@ fn report_epoch_completion( /* * Reports statistical analysis results for SSx observations */ -fn report_ssi_statistics( - ssi_stats: &HashMap, -) -> Box { +fn report_ssi_statistics(ssi_stats: &HashMap) -> Box { box_html! { table(class="table is-bordered") { thead { @@ -218,7 +307,7 @@ fn report_ssi_statistics( th { : "Mean" } - @ for (_, (mean, _, _)) in ssi_stats { + @ for (_, (mean, _)) in ssi_stats { td { : format!("{:.3} dB", mean) } @@ -228,19 +317,9 @@ fn report_ssi_statistics( th { : "Deviation" // (σ)" } - @ for (_, (_, sigma, _)) in ssi_stats { - td { - : format!("{:.3} dB", sigma) - } - } - } - tr { - th { - : "Skewness" //(μ / σγ)" - } - @ for (_, (_, _, sk)) in ssi_stats { + @ for (_, (_, std)) in ssi_stats { td { - : format!("{:.3}", sk) + : format!("{:.3} dB", std) } } } @@ -261,8 +340,12 @@ pub struct QcObsAnalysis { codes: Vec, /// true if doppler observations were identified has_doppler: bool, - /// Abnormal events, by chronological epochs - anomalies: Vec<(Epoch, EpochFlag)>, + /// CS anomalies + cs_anomalies: Vec, + /// Epochs where power failures took place, and their duration + power_failures: Vec, + /// Other abnormal events, by chronological epochs + other_anomalies: Vec<(Epoch, EpochFlag)>, /// Total number of epochs identified total_epochs: usize, /// Epochs with at least 1 observation @@ -270,87 +353,73 @@ pub struct QcObsAnalysis { /// Complete epochs, with respect to given signal complete_epochs: Vec<(Carrier, usize)>, #[cfg(feature = "obs")] - /// Min. Max. SNR (sv @ epoch) - min_max_snr: ((Sv, Epoch, Snr), (Sv, Epoch, Snr)), + /// Min. & Max. SNR (signal @ epoch) + snr_stats: HashMap, + #[cfg(feature = "obs")] + /// SSI statistical analysis (mean, stddev) + ssi_stats: HashMap, #[cfg(feature = "obs")] - /// SSi statistical analysis (mean, stddev, skew) - ssi_stats: HashMap, - #[cfg(feature = "processing")] - /// Receiver clock drift analysis - clock_drift: Option, + /// RX clock drift + clock_drift: Vec<(Epoch, f64)>, } impl QcObsAnalysis { pub fn new(rnx: &Rinex, opts: &QcOpts) -> Self { - let sv: Vec<_> = rnx.sv().collect(); - let obs = rnx.header.obs.as_ref().unwrap(); - let mut observables = obs.codes.clone(); - // TODO: improve this .get(0) - let observables = observables.get_mut(&sv[0].constellation).unwrap(); - let mut signals: Vec = Vec::new(); - let mut codes: Vec = Vec::new(); - let mut anomalies: Vec<(Epoch, EpochFlag)> = Vec::new(); - let mut total_epochs: usize = 0; - let mut epoch_with_obs: Vec = Vec::new(); - let mut complete_epochs: HashMap = HashMap::new(); - let mut min_max_snr = ( - (Sv::default(), Epoch::default(), Snr::DbHz54), - (Sv::default(), Epoch::default(), Snr::DbHz0), - ); - let mut ssi_stats: HashMap = HashMap::new(); + let doppler_obs = rnx.observable().filter(|obs| obs.is_doppler_observable()); - let mut clock_drift: Option = None; + let mut observables: Vec = rnx.observable().map(|obs| obs.to_string()).collect(); - if let Some(rec) = rnx.record.as_obs() { - let mask: Filter = Filter::from(MaskFilter { - operand: MaskOperand::Equals, - item: TargetItem::ClockItem, - }); - let clk_data = rec.filter(mask); - //TODO: run analysis on smoothed derivative - //let der = clk_data.derivative(); - //let mov = der.moving_average(Duration::from_seconds(600.0), None); - } + let mut signals: Vec<_> = rnx.carrier().unique().collect(); + let mut codes: Vec<_> = rnx.code().map(|c| c.to_string()).collect(); - if let Some(r) = rnx.record.as_obs() { - total_epochs = r.len(); - for ((epoch, flag), (clk, svs)) in r { - if let Some(_clk) = clk {} + let cs_anomalies: Vec<_> = rnx + .epoch_anomalies() + .filter_map(|(e, flag)| { + if flag == EpochFlag::CycleSlip { + Some(e) + } else { + None + } + }) + .collect(); - if !flag.is_ok() { - anomalies.push((*epoch, *flag)); + let power_failures: Vec<_> = rnx + .epoch_anomalies() + .filter_map(|(e, flag)| { + if flag == EpochFlag::PowerFailure { + Some(e) + } else { + None } - for (sv, observables) in svs { - if observables.len() > 0 && !epoch_with_obs.contains(&epoch) { - epoch_with_obs.push(*epoch); - } + }) + .collect(); - for (observable, observation) in observables { - let code = observable.code().unwrap(); - let carrier = observable.carrier(sv.constellation).unwrap(); - if !signals.contains(&carrier) { - signals.push(carrier); - } - if !codes.contains(&code) { - codes.push(code); - } + let other_anomalies: Vec<_> = rnx + .epoch_anomalies() + .filter_map(|(e, flag)| { + if flag != EpochFlag::PowerFailure && flag != EpochFlag::CycleSlip { + Some((e, flag)) + } else { + None + } + }) + .collect(); - if let Some(snr) = observation.snr { - if snr < min_max_snr.0 .2 { - min_max_snr.0 .0 = *sv; - min_max_snr.0 .1 = *epoch; - min_max_snr.0 .2 = snr; - } - if snr > min_max_snr.1 .2 { - min_max_snr.1 .0 = *sv; - min_max_snr.1 .1 = *epoch; - min_max_snr.1 .2 = snr; - } + let mut total_epochs = rnx.epoch().count(); + let mut epoch_with_obs: Vec = Vec::new(); + let mut complete_epochs: HashMap = HashMap::new(); + + if let Some(r) = rnx.record.as_obs() { + total_epochs = r.len(); + for ((epoch, _flag), (_clk, svs)) in r { + for (_sv, observables) in svs { + if !observables.is_empty() { + if !epoch_with_obs.contains(&epoch) { + epoch_with_obs.push(*epoch); } } } } - /* * Now that signals have been determined, * determine observation completion @@ -376,7 +445,6 @@ impl QcObsAnalysis { continue; // phase should have SNR information attached to it } } - /* * Signal condition */ @@ -423,52 +491,64 @@ impl QcObsAnalysis { } } } - - /* - * SSI statistical analysis - */ - let mut mean_ssi: HashMap<_, _> = r.mean_observable(); - mean_ssi.retain(|obs, _| obs.is_ssi_observable()); - for (obs, mean) in mean_ssi { - ssi_stats.insert(obs.clone(), (mean, 0.0_f64, 0.0_f64)); - } - - let mut stddev_ssi: HashMap<_, _> = r.mean_observable(); // TODO r.stddev_observable(); - stddev_ssi.retain(|obs, _| obs.is_ssi_observable()); - for (obs, stddev) in stddev_ssi { - if let Some((_, dev, _)) = ssi_stats.get_mut(&obs) { - *dev = stddev; - } + } + // append ssi: drop vehicle differentiation + let mut ssi: HashMap> = HashMap::new(); + for (_, _, obs, value) in rnx.ssi() { + if let Some(values) = ssi.get_mut(&obs) { + values.push(value); + } else { + ssi.insert(obs.clone(), vec![value]); } - - let mut skew_ssi: HashMap<_, _> = r.mean_observable(); // TODO r.skewness_observable(); - skew_ssi.retain(|obs, _| obs.is_ssi_observable()); - for (obs, skew) in skew_ssi { - if let Some((_, _, sk)) = ssi_stats.get_mut(&obs) { - *sk = skew; - } + } + /* + * SSI statistical analysis: {mean, stddev,} + * per signal: we do not differentiate vehicles + */ + let ssi_stats: HashMap = ssi + .iter() + .map(|(obs, values)| (obs.clone(), (values.mean(), values.std_dev()))) + .collect(); + // append snr: drop vehicle differentiation + let mut snr: HashMap> = HashMap::new(); + for ((e, _), _, obs, snr_value) in rnx.snr() { + let snr_f64: f64 = (snr_value as u8).into(); + if let Some(values) = snr.get_mut(&obs) { + values.push((e, snr_f64)); + } else { + snr.insert(obs.clone(), vec![(e, snr_f64)]); } } - + /* + * SNR analysis: {min, max} + * per signal: we do not differentiate vehicles + */ + let mut snr_stats: HashMap = HashMap::new(); + for (obs, data) in snr { + let values: Vec = data.iter().map(|(_e, value)| *value).collect(); + let min = values.clone().min(); + println!("MIN: {}", min); + let epoch_min = data.iter().find(|(_e, value)| *value == min).unwrap().0; + let max = values.clone().max(); + println!("MAX: {}", max); + let epoch_max = data.iter().find(|(_e, value)| *value == max).unwrap().0; + snr_stats.insert(obs, ((epoch_min, min), (epoch_max, max))); + } + /* + * sort, prior reporting + */ codes.sort(); - signals.sort(); observables.sort(); + signals.sort(); Self { - observables: { observables.iter().map(|v| v.to_string()).collect() }, - has_doppler: { - let mut ret = false; - for obs in observables.iter() { - if obs.is_doppler_observable() { - ret = true; - break; - } - } - ret - }, codes, signals, - anomalies, + observables, + has_doppler: doppler_obs.count() > 0, + cs_anomalies, + power_failures, + other_anomalies, total_epochs, total_with_obs: epoch_with_obs.len(), complete_epochs: { @@ -477,9 +557,20 @@ impl QcObsAnalysis { ret.sort(); ret }, - min_max_snr, + snr_stats, ssi_stats, - clock_drift, + clock_drift: { + let rx_clock: Vec<_> = rnx + .recvr_clock() + .map(|((e, _flag), value)| (e, value)) + .collect(); + let der = Derivative::new(1); + let rx_clock_drift: Vec<(Epoch, f64)> = der.eval(rx_clock); + //TODO + //let mov = Averager::mov(opts.clock_drift_window); + //mov.eval(rx_clock_drift) + rx_clock_drift + }, } } } @@ -532,56 +623,63 @@ impl HtmlReport for QcObsAnalysis { } } tr { - table { - : report_anomalies(&self.anomalies) + table(class="table is-bordered") { + thead { + th { + : "Anomalies" + } + } + tbody { + : report_anomalies(&self.cs_anomalies, &self.power_failures, &self.other_anomalies) + } } } tr { - table { - : report_epoch_completion(self.total_epochs, self.total_with_obs, &self.complete_epochs) - } - } - table(class="table is-bordered") { - thead { - tr { + table(class="table is-bordered") { + thead { th { - : "Worst SNR" + : "Epochs" } + } + tbody { + : report_epoch_completion(self.total_epochs, self.total_with_obs, &self.complete_epochs) + } + } + } + tr { + table(class="table is-bordered") { + thead { th { - : "Best SNR" + : "SNR" } } + tbody { + : report_snr_statistics(&self.snr_stats) + } } - tbody { - tr { - td { - p { - : self.min_max_snr.0.0.to_string() - } - b { - : format!("{:e}", self.min_max_snr.0.2) - } - p { - : format!("@{}", self.min_max_snr.0.1) - } - } - td { - p { - : self.min_max_snr.1.0.to_string() - } - b { - : format!("{:e}", self.min_max_snr.1.2) - } - p { - : format!("@{}", self.min_max_snr.1.1) - } + } + tr { + table(class="table is-bordered") { + thead { + th { + : "SSI" } } + tbody { + : report_ssi_statistics(&self.ssi_stats) + } } } tr { - table { - : report_ssi_statistics(&self.ssi_stats) + table(class="table is-bordered") { + thead { + th { + : "(RX) Clock Drift" + } + } + tbody { + : report_clock_drift(&self.clock_drift) + } } } } diff --git a/rinex/src/qc/analysis/sampling.rs b/rinex/src/qc/analysis/sampling.rs index 3266f6962..4c175113e 100644 --- a/rinex/src/qc/analysis/sampling.rs +++ b/rinex/src/qc/analysis/sampling.rs @@ -90,15 +90,29 @@ impl HtmlReport for QcSamplingAnalysis { } tr { th { - : "Sample rate" + : "Sample rate (Header)" } @ if let Some(rate) = self.sample_rate { td { : format!("{} ({:.3} Hz)", rate, 1.0 / rate.to_unit(Unit::Second)) } } else { + th { + : "Unspecified" + } + } + } + tr { + th { + : "Dominant Sample rate" + } + @ if let Some(rate) = self.dominant_sample_rate { td { - : "Unknown" + : format!("{} ({:.3} Hz)", rate, 1.0 / rate.to_unit(Unit::Second)) + } + } else { + th { + : "Undetermined" } } } @@ -106,21 +120,19 @@ impl HtmlReport for QcSamplingAnalysis { th { : "Gap analysis" } - } - @ if self.gaps.len() == 0 { - tr { + + @ if self.gaps.is_empty() { th { - : "No Data Gaps detected" + : "No gaps detected" } - } - } else { - @ for (epoch, dt) in &self.gaps { + } else { tr { td { - : format!("Start : {}", epoch) - } - td { - : format!("Duration: {}", dt) + @ for (epoch, dt) in &self.gaps { + p { + : format!("Start : {}, Duration: {}", epoch, dt) + } + } } } } diff --git a/rinex/src/qc/analysis/sv.rs b/rinex/src/qc/analysis/sv.rs index 01996cbc5..dcce4fa05 100644 --- a/rinex/src/qc/analysis/sv.rs +++ b/rinex/src/qc/analysis/sv.rs @@ -1,5 +1,5 @@ use super::{pretty_array, QcOpts}; -use crate::prelude::*; +use crate::prelude::*; //table_lengthy_td #[derive(Debug, Clone)] pub struct QcSvAnalysis { diff --git a/rinex/src/qc/mod.rs b/rinex/src/qc/mod.rs index bbe5c9c0f..3648c6595 100644 --- a/rinex/src/qc/mod.rs +++ b/rinex/src/qc/mod.rs @@ -13,6 +13,36 @@ use analysis::QcAnalysis; #[cfg(feature = "processing")] use crate::preprocessing::*; +/* + * Methods used when reporting lenghty vectors or data subsets in a table. + * Makes tables cleaner and nicer by wrapping string content, into several paragraphs. + * +pub(crate) fn table_lengthy_td( + list: &Vec, + max_items: usize, +) -> Box { + let mut content = String::with_capacity(64 * max_items); + let mut paragraphs: Vec = Vec::new(); + + for i in 0..list.len() { + content.push_str(&format!("{}, ", list[i])); + if i.rem_euclid(max_items) == 0 { + paragraphs.push(content.clone()); + content.clear(); + } else if i == list.len() - 1 { + paragraphs.push(content.clone()); + } + } + box_html! { + @ for paragraph in paragraphs { + p { + : paragraph.to_string() + } + } + } +} +*/ + /* * Array (CSV) pretty formatter */ @@ -139,7 +169,7 @@ impl QcReport { h2(class="title") { : "RINEX Quality Check summary" } - table(class="table is-bordered") { + table(class="table is-bordered; style=\"margin-bottom: 20px\"") { tbody { tr { th { @@ -153,30 +183,36 @@ impl QcReport { } }//div=header div(id="context") { - h3(class="title") { - : "Context" - } - table(class="table is-bordered") { + table(class="table is-bordered; style=\"margin-bottom: 20px\"") { + thead { + th { + : "Context" + } + } tbody { : context.to_inline_html() } } }//div=context div(id="parameters") { - h3(class="title") { - : "Parameters" - } - table(class="table is-bordered") { + table(class="table is-bordered; style=\"margin-bottom: 20px\"") { + thead { + th { + : "Parameters" + } + } tbody { : opts.to_inline_html() } } } //div=parameters div(id="header") { - h3(class="title") { - : "File Header" - } - table(class="table is-bordered") { + table(class="table is-bordered; style=\"margin-bottom: 20px\"") { + thead { + th { + : "File Header" + } + } tbody { : context.primary_data().header.to_inline_html() } @@ -186,10 +222,30 @@ impl QcReport { * Report all analysis that were performed */ div(id="analysis") { - @ for analysis in &analysis { - table(class="table is-bordered") { + /* + * Report all analysis + * and emphasize how they were sorted (self.opts.classfication) + */ + @ for i in 0..analysis.len() { + table(class="table is-bordered; style=\"margin-bottom: 20px\"") { + thead { + @ if opts.classification == QcClassification::GNSS { + th { + : format!("{} analysis", context.primary_data().constellation().nth(i).unwrap()) + } + } else if opts.classification == QcClassification::Sv { + th { + : format!("{} analysis", context.primary_data().sv().nth(i).unwrap()) + } + + } else if opts.classification == QcClassification::Physics { + th { + : format!("{} analysis", context.primary_data().observable().nth(i).unwrap()) + } + } + } tbody { - : analysis.to_inline_html() + : analysis[i].to_inline_html() } } } diff --git a/rinex/src/qc/opts.rs b/rinex/src/qc/opts.rs index 61ce82b90..5bb3acb2e 100644 --- a/rinex/src/qc/opts.rs +++ b/rinex/src/qc/opts.rs @@ -119,7 +119,7 @@ impl Default for ProcessingOpts { } /// Qc Report classification method -#[derive(Default, Debug, Clone)] +#[derive(Default, Debug, Clone, PartialEq)] #[cfg_attr(feature = "serde", derive(Serialize, Deserialize))] /// Classify the QC report by desired data set pub enum QcClassification { @@ -161,6 +161,9 @@ pub struct QcOpts { pub gap_tolerance: Option, /// Manually defined Ground position (ECEF) pub ground_position: Option, + /// Window duration to be used, during RX clock drift analysis + #[serde(default = "default_drift_window")] + pub clock_drift_window: Duration, } impl QcOpts { @@ -190,6 +193,10 @@ impl QcOpts { } } +fn default_drift_window() -> Duration { + Duration::from_seconds(3600.0) +} + impl Default for QcOpts { fn default() -> Self { Self { @@ -198,6 +205,7 @@ impl Default for QcOpts { min_snr_db: 20.0, // dB elev_mask: None, classification: QcClassification::default(), + clock_drift_window: default_drift_window(), } } } @@ -243,15 +251,23 @@ impl HtmlReport for QcOpts { : "Data gap" } @ if let Some(tol) = self.gap_tolerance { - th { - : format!("Tolerance: ({})", tol) + td { + : format!("{} tolerance", tol) } } else { - th { + td { : "No tolerance" } } } + tr { + th { + : "Clock Drift Window" + } + td { + : self.clock_drift_window.to_string() + } + } } } } diff --git a/rinex/src/qc/sampling.rs b/rinex/src/qc/sampling.rs deleted file mode 100644 index 7ad496aff..000000000 --- a/rinex/src/qc/sampling.rs +++ /dev/null @@ -1,128 +0,0 @@ -use crate::prelude::*; -use horrorshow::RenderBox; - -/// Sampling QC report -#[derive(Debug, Clone, PartialEq)] -#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))] -pub struct QcReport { - /// First epoch identified - pub first_epoch: Epoch, - /// Last epoch identified - pub last_epoch: Epoch, - /// Time line - pub time_line: Duration, - /// Dominant sample rate - pub sample_rate: Duration, - /// Unusual data gaps - pub gaps: Vec<(Epoch, Duration)>, -} - -impl QcReport { - pub fn new(rnx: &Rinex) -> Self { - let first_epoch = rnx - .first_epoch() - .expect("Sampling QC expects a RINEX indexed by epochs"); - let last_epoch = rnx - .last_epoch() - .expect("Sampling QC expects a RINEX indexed by epochs"); - let sample_rate = rnx - .sampling_interval() - .expect("failed to determine sample rate"); - Self { - first_epoch, - last_epoch, - sample_rate, - time_line: last_epoch - first_epoch, - gaps: rnx.data_gaps(), - } - } - fn gap_analysis(gaps: &Vec<(Epoch, Duration)>) -> Box { - box_html! { - @ if gaps.len() == 0 { - tr { - th { - b { - : "Gap analysis: " - } - } - td { - : "Data missing" - } - } - } else { - tr { - td { - b { - : "Gap Aanalysis:" - } - } - } - tr { - td { - : "" - } - td { - : "Epoch (start)" - } - td { - : "Duration" - } - } - @ for (epoch, duration) in gaps { - tr { - td { - : "" - } - td { - : epoch.to_string() - } - td { - : duration.to_string() - } - } - } - } - } - } - pub fn to_inline_html(&self) -> Box { - box_html! { - table(id="sampling") { - tr { - td { - : "First Epoch:" - } - td { - : self.first_epoch.to_string() - } - } - tr { - td { - : "Last Epoch:" - } - td { - : self.last_epoch.to_string() - } - } - tr { - td { - : "Time line:" - } - td { - : self.time_line.to_string() - } - } - tr { - td { - : "Sample rate:" - } - td { - : self.sample_rate.to_string() - } - } - table(id="gap-analysis") { - : Self::gap_analysis(&self.gaps) - } - } - } - } -} diff --git a/rinex/src/tests/mod.rs b/rinex/src/tests/mod.rs index a92ab8e64..dacb79326 100644 --- a/rinex/src/tests/mod.rs +++ b/rinex/src/tests/mod.rs @@ -17,4 +17,3 @@ mod parsing; mod production; mod sampling; mod smoothing; -mod stats; diff --git a/rinex/src/tests/obs.rs b/rinex/src/tests/obs.rs index 4ab469df9..6af66d117 100644 --- a/rinex/src/tests/obs.rs +++ b/rinex/src/tests/obs.rs @@ -1184,7 +1184,7 @@ mod test { let mut irnss_sv: Vec = rnx .sv() .filter_map(|sv| { - if (sv.constellation == Constellation::IRNSS) { + if sv.constellation == Constellation::IRNSS { Some(sv) } else { None diff --git a/rinex/src/tests/stats.rs b/rinex/src/tests/stats.rs deleted file mode 100644 index 5de0dfdbc..000000000 --- a/rinex/src/tests/stats.rs +++ /dev/null @@ -1,347 +0,0 @@ -#[cfg(test)] -mod test { - use crate::*; - use crate::{observation::*, prelude::*}; - use std::str::FromStr; - fn run_test(rinex: &Rinex, ops: &str, expected: Vec<(Sv, Vec<(Observable, f64)>)>) { - let (clk_stats, stats) = match ops { - "min" => rinex.min(), - "max" => rinex.max(), - "mean" => rinex.mean(), - _ => panic!("unknown test"), - }; - assert!(clk_stats.is_none()); // always none on our test pool - - for (sv, fields) in expected { - let sv_data = stats.get(&sv); - assert!( - sv_data.is_some(), - "{}", - format!("\"{}\" stats missing for vehicle \"{}\"", ops, sv) - ); - let sv_data = sv_data.unwrap(); - for (observable, expected) in fields { - let data = sv_data.get(&observable); - assert!( - data.is_some(), - "{}", - format!( - "\"{}\" stats missing for vehicle \"{}\" - \"{}\"", - ops, sv, observable - ) - ); - let data = data.unwrap(); - let e = *data - expected; - assert!( - e.abs() < 1E-5, - "{}", - format!( - "\"{}\" test failed for \"{}\" - \"{}\", expecting {} got {}", - ops, sv, observable, expected, *data - ) - ); - } - } - } - fn run_obsv_test(rinex: &Rinex, ops: &str, expected: Vec<(Observable, f64)>) { - let stats = match ops { - "min" => rinex.min_observable(), - "max" => rinex.max_observable(), - "mean" => rinex.mean_observable(), - _ => panic!("unknown test"), - }; - - for (observable, expected) in expected { - let results = stats.get(&observable); - assert!( - results.is_some(), - "{}", - format!( - "\"{}\" observable stats results missing for observable \"{}\"", - ops, observable - ) - ); - let results = results.unwrap(); - let err = *results - expected; - assert!( - err.abs() < 1E-5, - "{}", - format!( - "\"{}\" observable test failed for - \"{}\", expecting {} - got {}", - ops, observable, expected, results - ) - ); - } - } - #[test] - fn stats() { - let rinex = Rinex::from_file("../test_resources/OBS/V3/DUTH0630.22O").unwrap(); - - let expected: Vec<(Sv, Vec<(Observable, f64)>)> = vec![ - ( - sv!("g01"), - vec![ - ( - observable!("c1c"), - (20243517.560 + 20805393.080 + 21653418.260) / 3.0, - ), - (observable!("d1c"), (-1242.766 - 2193.055 - 2985.516) / 3.0), - (observable!("s1c"), (51.25 + 50.75 + 49.5) / 3.0), - ], - ), - ( - sv!("g03"), - vec![ - ( - observable!("c1c"), - (20619020.680 + 20425456.580 + 20410261.460) / 3.0, - ), - (observable!("d1c"), (852.785 + 328.797 - 244.031) / 3.0), - (observable!("s1c"), (50.75 + 51.0 + 51.0) / 3.0), - ], - ), - ]; - run_test(&rinex, "mean", expected); - - let expected: Vec<(Sv, Vec<(Observable, f64)>)> = vec![ - ( - sv!("g01"), - vec![ - (observable!("c1c"), 20243517.560), - (observable!("d1c"), -2985.516), - (observable!("s1c"), 49.5), - ], - ), - ( - sv!("g03"), - vec![ - (observable!("c1c"), 20410261.460), - (observable!("d1c"), -244.031), - (observable!("s1c"), 50.75), - ], - ), - ]; - run_test(&rinex, "min", expected); - - let expected: Vec<(Sv, Vec<(Observable, f64)>)> = vec![ - ( - sv!("g01"), - vec![ - (observable!("c1c"), 21653418.260), - (observable!("d1c"), -1242.766), - (observable!("s1c"), 51.25), - ], - ), - ( - sv!("g03"), - vec![ - (observable!("c1c"), 20619020.680), - (observable!("d1c"), 852.785), - (observable!("s1c"), 51.0), - ], - ), - ]; - run_test(&rinex, "max", expected); - - /* ********************************************* - * STATS FOR ALL SV / ALL OBS : meteo compatible - **********************************************/ - /*let expected : Vec<(Observable, f64)> = vec![ - (observable!("c1c"), (20243517.560 + 20619020.680 + 21542633.500 + 24438727.980 + 22978068.560 + 23460759.840 + 21923317.180 + 23434790.440 + 22401985.340 + 24991723.280 + 19727826.340 + 23171275.620 + 20662538.580 + 23450513.820 + 23044984.180 + 22909354.040 + 20116780.920 + 19708379.260 + 20805393.080 + 20425456.580 + 20887001.400 + 23371156.300 + 23031543.660 + 23117350.280 + 22726604.680 + 24425563.640 + 22689941.780 + 19677287.000 + 22265147.080 + 21462395.740 + 23740237.340 + 22432243.520 + 21750541.080 + 21199384.320 + 19680274.400 + 21653418.260 + 20410261.460 + 20488105.720 + 23647940.540 + 22436978.380 + 23392660.200 + 23154069.760 + 23689895.760 + 25161827.280 + 23333751.720 + 19831816.600 + 21490078.220 + 22328018.960 + 22235350.560 + 20915624.780 + 22543866.020 + 20147683.700) / 52.0), - ]; - run_obsv_test(&rinex, "mean", expected);*/ - - /* - * Test on METEO RINEX - */ - let rinex = Rinex::from_file("../test_resources/MET/V2/clar0020.00m").unwrap(); - let expected: Vec<(Observable, f64)> = vec![( - observable!("PR"), - (970.5 - + 970.4 - + 970.4 - + 970.3 - + 970.2 - + 970.3 - + 970.1 - + 970.3 - + 970.2 - + 970.2 - + 972.1 - + 972.3 - + 972.4 - + 972.6 - + 972.8 - + 973.1 - + 973.3 - + 973.4 - + 973.6 - + 973.6 - + 973.7 - + 973.7 - + 973.5 - + 973.6 - + 973.6 - + 973.6 - + 973.6 - + 973.5 - + 973.5 - + 973.4 - + 973.3 - + 973.2 - + 972.9 - + 972.8 - + 972.7 - + 972.5 - + 972.4 - + 972.4 - + 972.2 - + 972.2 - + 972.2 - + 972.3 - + 972.2 - + 972.2 - + 972.1 - + 972.2 - + 972.1 - + 972.2 - + 972.2 - + 972.2 - + 972.1 - + 972.2 - + 972.2 - + 972.3 - + 972.3 - + 972.5 - + 972.5) - / 57.0, - )]; - run_obsv_test(&rinex, "mean", expected); - - let expected: Vec<(Observable, f64)> = vec![( - observable!("TD"), - (10.7 - + 10.6 - + 10.4 - + 10.3 - + 10.1 - + 9.9 - + 9.8 - + 9.5 - + 9.6 - + 9.6 - + 8.4 - + 10.1 - + 10.1 - + 10.3 - + 10.2 - + 10.3 - + 10.6 - + 10.8 - + 11.1 - + 11.2 - + 11.4 - + 11.5 - + 11.8 - + 11.8 - + 11.8 - + 12.7 - + 12.7 - + 13.2 - + 13.5 - + 13.0 - + 13.5 - + 13.9 - + 14.0 - + 14.7 - + 15.0 - + 14.8 - + 15.6 - + 15.1 - + 14.7 - + 14.8 - + 14.9 - + 14.9 - + 15.9 - + 15.6 - + 16.2 - + 15.6 - + 15.5 - + 15.9 - + 15.8 - + 15.7 - + 15.7 - + 15.6 - + 15.4 - + 14.9 - + 14.6 - + 14.3 - + 14.2) - / 57.0, - )]; - run_obsv_test(&rinex, "mean", expected); - - let expected: Vec<(Observable, f64)> = vec![( - observable!("HR"), - (71.4 - + 72.2 - + 72.9 - + 72.9 - + 75.0 - + 76.7 - + 78.1 - + 79.7 - + 80.5 - + 79.1 - + 70.7 - + 65.3 - + 66.4 - + 68.0 - + 68.4 - + 68.6 - + 66.8 - + 65.4 - + 63.3 - + 62.7 - + 60.0 - + 57.2 - + 56.0 - + 54.9 - + 52.2 - + 49.3 - + 46.3 - + 43.8 - + 44.8 - + 42.7 - + 41.2 - + 40.2 - + 37.9 - + 36.4 - + 37.3 - + 36.1 - + 32.7 - + 31.6 - + 36.1 - + 35.7 - + 34.0 - + 32.6 - + 29.5 - + 28.0 - + 26.7 - + 25.5 - + 25.0 - + 30.7 - + 25.0 - + 23.5 - + 27.7 - + 28.2 - + 28.4 - + 30.9 - + 33.3 - + 33.3 - + 33.2) - / 57.0, - )]; - run_obsv_test(&rinex, "mean", expected); - } -}