From 00061fb176d1fdc9befb05b2dfc455827d18c1f4 Mon Sep 17 00:00:00 2001 From: "Guillaume W. Bres" Date: Sun, 17 Aug 2025 19:58:36 +0200 Subject: [PATCH 01/25] * fix docs * export BoundingRect trait Signed-off-by: Guillaume W. Bres --- src/lib.rs | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 2336d34..aefe89b 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -94,8 +94,8 @@ pub mod prelude { // pub re-export pub use geo::{ - algorithm::contains::Contains, coord, GeodesicArea, Geometry, LineString, Point, Polygon, - Rect, + algorithm::contains::Contains, coord, BoundingRect, GeodesicArea, Geometry, LineString, + Point, Polygon, Rect, }; pub use gnss::prelude::{Constellation, SV}; pub use hifitime::{Duration, Epoch, TimeScale, TimeSeries, Unit}; @@ -512,7 +512,7 @@ impl IONEX { false } - /// Returns map borders as [Retc]angle. This uses + /// Returns map borders as [Rect]angle. This uses /// the [Header] description and assumes all maps are within these borders. pub fn map_borders_degrees(&self) -> Rect { Rect::new( @@ -549,6 +549,7 @@ impl IONEX { let mut ionex = IONEX::default(); let bounding_rect = region.bounding_rect()?; + let (min_long, min_lat) = (bounding_rect.min().x, bounding_rect.min().y); let (max_long, max_lat) = (bounding_rect.max().x, bounding_rect.max().y); From 9af4319c0b33229a0d0ee0a73a19e00beaf7f278 Mon Sep 17 00:00:00 2001 From: "Guillaume W. Bres" Date: Sat, 1 Nov 2025 10:55:19 +0100 Subject: [PATCH 02/25] update deps Signed-off-by: Guillaume W. Bres --- README.md | 23 ++++++++++++++++++++++- 1 file changed, 22 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index dcdb53b..ec8b024 100644 --- a/README.md +++ b/README.md @@ -20,7 +20,7 @@ quantization spec that is constant over entire fileset and is described in the h ## Advantages - Fast and powerful parser -- Open source +- Open sources: read and access all the code - Seamless Gzip decompression (on `flate2` feature) - Full 2D support - TEC Root Mean Square is supported @@ -45,6 +45,27 @@ Contributions are welcomed: - follow our [Discussions on Github.com](https://github.com/nav-solutions/discussions) - join our [Discord channel](https://discord.gg/EqhEBXBmJh) +RINEX formats & applications +============================ + +| Type | Parser | Writer | Content | Record Indexing | Timescale | +|----------------------------|-------------------------------------------------------------------------|---------------------|-----------------------------------------------|----------------------------------------------------------------------------------| -----------| +| Navigation (NAV) | [Provided by the RINEX parser]((https://github.com/nav-solutions/rinex) | :heavy_check_mark: | Ephemerides, Ionosphere models | [NavKey](https://docs.rs/rinex/latest/rinex/navigation/struct.NavKey.html) | SV System time broadcasting this message | +| Observation (OBS) | [Provided by the RINEX parser]((https://github.com/nav-solutions/rinex) | :heavy_check_mark: | Phase, Pseudo Range, Doppler, SSI | [ObsKey](https://docs.rs/rinex/latest/rinex/observation/struct.ObsKey.html) | GNSS (any) | +| CRINEX (Compressed OBS) | [Provided by the RINEX parser]((https://github.com/nav-solutions/rinex) | :heavy_check_mark: | Phase, Pseudo Range, Doppler, SSI | [ObsKey](https://docs.rs/rinex/latest/rinex/observation/struct.ObsKey.html) | GNSS (any) | +| Meteorological data (MET) | [Provided by the RINEX parser]((https://github.com/nav-solutions/rinex) | :heavy_check_mark: | Meteo sensors data (Temperature, Moisture..) | [MeteoKey](https://docs.rs/rinex/latest/rinex/meteo/struct.MeteoKey.html) | UTC | +| Clocks (CLK) | [Provided by the RINEX parser]((https://github.com/nav-solutions/rinex) | :construction: | Precise temporal states | [ClockKey](https://docs.rs/rinex/latest/rinex/clock/record/struct.ClockKey.html) | GNSS (any) | +| Antenna (ATX) | [Provided by the RINEX parser]((https://github.com/nav-solutions/rinex) | :construction: | Precise RX/SV Antenna calibration | `antex::Antenna` | :heavy_minus_sign: | +| Ionosphere Maps (IONEX) | :heavy_check_mark: | :heavy_check_mark: | Ionosphere Electron density | [Record Key](https://docs.rs/ionex/latest/ionex/key/struct.Key.html) | UTC | +| DORIS RINEX | [Provided by the DORIS parser](https://github.com/nav-solutions/doris) | :heavy_check_mark: | Temperature, Moisture, Pseudo Range and Phase observations | [Record Key](https://docs.rs/doris-rs/latest/doris_rs/record/struct.Key.html) | TAI / "DORIS" timescale | + +Contributions +============= + +Contributions are welcomed, we still have a lot to accomplish, any help is always appreciated. +[We wrote these few lines](CONTRIBUTING.md) to help you understand the inner workings. +Join us on [Discord](https://discord.gg/EqhEBXBmJh) to discuss ongoing and future developments. + ## Getting started ```rust From 1a15db293c30cfc044b0fa95c6fe8b24657e5480 Mon Sep 17 00:00:00 2001 From: "Guillaume W. Bres" Date: Sat, 1 Nov 2025 10:56:16 +0100 Subject: [PATCH 03/25] update deps Signed-off-by: Guillaume W. Bres --- README.md | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/README.md b/README.md index ec8b024..3f29330 100644 --- a/README.md +++ b/README.md @@ -50,14 +50,14 @@ RINEX formats & applications | Type | Parser | Writer | Content | Record Indexing | Timescale | |----------------------------|-------------------------------------------------------------------------|---------------------|-----------------------------------------------|----------------------------------------------------------------------------------| -----------| -| Navigation (NAV) | [Provided by the RINEX parser]((https://github.com/nav-solutions/rinex) | :heavy_check_mark: | Ephemerides, Ionosphere models | [NavKey](https://docs.rs/rinex/latest/rinex/navigation/struct.NavKey.html) | SV System time broadcasting this message | -| Observation (OBS) | [Provided by the RINEX parser]((https://github.com/nav-solutions/rinex) | :heavy_check_mark: | Phase, Pseudo Range, Doppler, SSI | [ObsKey](https://docs.rs/rinex/latest/rinex/observation/struct.ObsKey.html) | GNSS (any) | -| CRINEX (Compressed OBS) | [Provided by the RINEX parser]((https://github.com/nav-solutions/rinex) | :heavy_check_mark: | Phase, Pseudo Range, Doppler, SSI | [ObsKey](https://docs.rs/rinex/latest/rinex/observation/struct.ObsKey.html) | GNSS (any) | -| Meteorological data (MET) | [Provided by the RINEX parser]((https://github.com/nav-solutions/rinex) | :heavy_check_mark: | Meteo sensors data (Temperature, Moisture..) | [MeteoKey](https://docs.rs/rinex/latest/rinex/meteo/struct.MeteoKey.html) | UTC | -| Clocks (CLK) | [Provided by the RINEX parser]((https://github.com/nav-solutions/rinex) | :construction: | Precise temporal states | [ClockKey](https://docs.rs/rinex/latest/rinex/clock/record/struct.ClockKey.html) | GNSS (any) | -| Antenna (ATX) | [Provided by the RINEX parser]((https://github.com/nav-solutions/rinex) | :construction: | Precise RX/SV Antenna calibration | `antex::Antenna` | :heavy_minus_sign: | -| Ionosphere Maps (IONEX) | :heavy_check_mark: | :heavy_check_mark: | Ionosphere Electron density | [Record Key](https://docs.rs/ionex/latest/ionex/key/struct.Key.html) | UTC | -| DORIS RINEX | [Provided by the DORIS parser](https://github.com/nav-solutions/doris) | :heavy_check_mark: | Temperature, Moisture, Pseudo Range and Phase observations | [Record Key](https://docs.rs/doris-rs/latest/doris_rs/record/struct.Key.html) | TAI / "DORIS" timescale | +| Navigation (NAV) | [Provided by the RINEX parser](https://github.com/nav-solutions/rinex) | :heavy_check_mark: | Ephemerides, Ionosphere models | [NavKey](https://docs.rs/rinex/latest/rinex/navigation/struct.NavKey.html) | SV System time broadcasting this message | +| Observation (OBS) | [Provided by the RINEX parser](https://github.com/nav-solutions/rinex) | :heavy_check_mark: | Phase, Pseudo Range, Doppler, SSI | [ObsKey](https://docs.rs/rinex/latest/rinex/observation/struct.ObsKey.html) | GNSS (any) | +| CRINEX (Compressed OBS) | [Provided by the RINEX parser](https://github.com/nav-solutions/rinex) | :heavy_check_mark: | Phase, Pseudo Range, Doppler, SSI | [ObsKey](https://docs.rs/rinex/latest/rinex/observation/struct.ObsKey.html) | GNSS (any) | +| Meteorological data (MET) | [Provided by the RINEX parser](https://github.com/nav-solutions/rinex) | :heavy_check_mark: | Meteo sensors data (Temperature, Moisture..) | [MeteoKey](https://docs.rs/rinex/latest/rinex/meteo/struct.MeteoKey.html) | UTC | +| Clocks (CLK) | [Provided by the RINEX parser](https://github.com/nav-solutions/rinex) | :construction: | Precise temporal states | [ClockKey](https://docs.rs/rinex/latest/rinex/clock/record/struct.ClockKey.html) | GNSS (any) | +| Antenna (ATX) | [Provided by the RINEX parser](https://github.com/nav-solutions/rinex) | :construction: | Precise RX/SV Antenna calibration | `antex::Antenna` | :heavy_minus_sign: | +| Ionosphere Maps (IONEX) | :heavy_check_mark: | :heavy_check_mark: | Ionosphere Electron density | [Record Key](https://docs.rs/ionex/latest/ionex/key/struct.Key.html) | UTC | +| DORIS RINEX | [Provided by the DORIS parser](https://github.com/nav-solutions/doris) | :heavy_check_mark: | Temperature, Moisture, Pseudo Range and Phase observations | [Record Key](https://docs.rs/doris-rs/latest/doris_rs/record/struct.Key.html) | TAI / "DORIS" timescale | Contributions ============= From f5a0ba85d0920727413167e908f017e85209554f Mon Sep 17 00:00:00 2001 From: "Guillaume W. Bres" Date: Sat, 1 Nov 2025 11:41:03 +0100 Subject: [PATCH 04/25] Improve specifications and features * removed num_integer::div_ceil Signed-off-by: Guillaume W. Bres --- Cargo.toml | 38 ++++++++++++++++++++++++-------------- src/lib.rs | 20 ++++++++++++++++---- 2 files changed, 40 insertions(+), 18 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 58b637f..46f6680 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -22,28 +22,37 @@ all-features = true rustdoc-args = ["--cfg", "docrs", "--generate-link-to-definition"] [features] -default = ["flate2"] # gzip files supported by default +# serde by default +# gzip files supported by default +default = ["serde", "flate2"] -[build-dependencies] -serde_json = { version = "1.0", features = ["preserve_order"] } -serde = { version = "1.0", default-features = false, features = ["derive"] } +serde = [ + "dep:serde", + "gnss-rs/serde", + "hifitime/serde", +] + +[dependencies.gnss-rs] +git = "https://github.com/nav-solutions/gnss" +rev = "dc4d4c2d413a3be90a3fa08a6ab29079eec13923" +features = ["std", "domes", "cospar", "sbas"] + +[dependencies.hifitime] +version = "4.1" +features = ["std"] [dependencies] -geo = "0.30" +geo = "0.31" thiserror = "2" -itertools = "0.14.0" -num-integer = "0.1.44" # TODO: see if we can get rid of div_ceil - -# activate parser logs +itertools = "0.14" log = { version = "0.4", optional = true } - -# support gzip files flate2 = { version = "1", optional = true } - -gnss-rs = { version = "2.4", features = ["serde"] } -hifitime = { version = "4.1", features = ["serde", "std"] } serde = { version = "1.0", optional = true, default-features = false, features = ["derive"] } +[build-dependencies] +serde_json = { version = "1.0", features = ["preserve_order"] } +serde = { version = "1.0", default-features = false, features = ["derive"] } + [dev-dependencies] log = "0.4" criterion = "0.7" @@ -56,3 +65,4 @@ harness = false [[bench]] name = "formatting" harness = false + diff --git a/src/lib.rs b/src/lib.rs index aefe89b..1edc1d7 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -104,13 +104,26 @@ pub mod prelude { /// IONEX comments are readable descriptions. pub type Comments = Vec; +/// Divides provided usize, returns ceil rounded value. +fn div_ceil(value: usize, divider: usize) -> usize { + let q = value.div_euclid(divider); + let r = value.rem_euclid(divider); + if r == 0 { + q + } else { + q + 1 + } +} + /// macro to format one header line or a comment pub(crate) fn fmt_ionex(content: &str, marker: &str) -> String { if content.len() < 60 { format!("{: Date: Sat, 1 Nov 2025 12:10:23 +0100 Subject: [PATCH 05/25] header: qc merge Signed-off-by: Guillaume W. Bres --- Cargo.toml | 12 +++++++ src/header/mod.rs | 3 ++ src/header/qc.rs | 88 +++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 103 insertions(+) create mode 100644 src/header/qc.rs diff --git a/Cargo.toml b/Cargo.toml index 46f6680..6feb7cd 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -32,11 +32,22 @@ serde = [ "hifitime/serde", ] +qc = [ + "dep:maud", + "dep:gnss-qc-traits", +] + [dependencies.gnss-rs] git = "https://github.com/nav-solutions/gnss" rev = "dc4d4c2d413a3be90a3fa08a6ab29079eec13923" features = ["std", "domes", "cospar", "sbas"] +[dependencies.gnss-qc-traits] +git = "https://github.com/nav-solutions/qc-traits" +rev = "63200d41bab2efa7d10082ae43a113d360722342" +optional = true +features = ["html"] + [dependencies.hifitime] version = "4.1" features = ["std"] @@ -47,6 +58,7 @@ thiserror = "2" itertools = "0.14" log = { version = "0.4", optional = true } flate2 = { version = "1", optional = true } +maud = { version = "0.26", optional = true } serde = { version = "1.0", optional = true, default-features = false, features = ["derive"] } [build-dependencies] diff --git a/src/header/mod.rs b/src/header/mod.rs index 3d5215f..53b1110 100644 --- a/src/header/mod.rs +++ b/src/header/mod.rs @@ -1,6 +1,9 @@ mod formatting; mod parsing; +#[cfg(feature = "serde")] +mod qc; + #[cfg(feature = "serde")] use serde::{Deserialize, Serialize}; diff --git a/src/header/qc.rs b/src/header/qc.rs new file mode 100644 index 0000000..e1bf38e --- /dev/null +++ b/src/header/qc.rs @@ -0,0 +1,88 @@ +use crate::prelude::Header; + +use gnss_qc_traits::{Merge, MergeError}; + +impl Merge for Header { + fn merge(&self, rhs: &Self) -> Result { + let mut s = self.clone(); + s.merge_mut(rhs)?; + Ok(s) + } + + fn merge_mut(&mut self, rhs: &Self) -> Result<(), MergeError> { + self.version = std::cmp::min(self.version, rhs.version); + + if self.program.is_none() { + if let Some(program) = &rhs.program { + self.program = Some(program.clone()); + } + } + + if self.run_by.is_none() { + if let Some(run_by) = &rhs.run_by { + self.run_by = Some(run_by.clone()); + } + } + + if self.date.is_none() { + if let Some(date) = &rhs.date { + self.date = Some(date.clone()); + } + } + + if self.license.is_none() { + if let Some(license) = &rhs.license { + self.license = Some(license.clone()); + } + } + + if self.doi.is_none() { + if let Some(doi) = &rhs.doi { + self.doi = Some(doi.clone()); + } + } + + match self.description { + Some(ref mut desc) => { + if let Some(rhs) = &rhs.description { + desc.push_str(rhs); + } + }, + None => { + if let Some(desc) = &rhs.description { + self.description = Some(desc.clone()); + } + }, + } + + self.epoch_of_first_map = std::cmp::min(self.epoch_of_first_map, rhs.epoch_of_first_map); + + self.epoch_of_last_map = std::cmp::max(self.epoch_of_last_map, rhs.epoch_of_last_map); + + if self.reference_system != rhs.reference_system { + // TODO: both must match + } + + if self.mapf != rhs.mapf { + // TODO: both must match + } + + if self.map_dimension != rhs.map_dimension { + // TODO: both must match + } + + self.sampling_period = std::cmp::min(self.sampling_period, rhs.sampling_period); + + if rhs.elevation_cutoff > self.elevation_cutoff { + self.elevation_cutoff = rhs.elevation_cutoff; + } + + // TODO: merge grid def + + for comment in rhs.comments.iter() { + self.comments.push(comment.clone()); + } + + Ok(()) + } +} From 63e8970b86cb43de95c5877a3e46e90c5997330a Mon Sep 17 00:00:00 2001 From: "Guillaume W. Bres" Date: Sat, 1 Nov 2025 12:14:18 +0100 Subject: [PATCH 06/25] record merge Signed-off-by: Guillaume W. Bres --- src/header/mod.rs | 2 +- src/record/mod.rs | 3 +++ src/record/qc.rs | 31 +++++++++++++++++++++++++++++++ src/tec.rs | 2 +- 4 files changed, 36 insertions(+), 2 deletions(-) create mode 100644 src/record/qc.rs diff --git a/src/header/mod.rs b/src/header/mod.rs index 53b1110..eaff34a 100644 --- a/src/header/mod.rs +++ b/src/header/mod.rs @@ -1,7 +1,7 @@ mod formatting; mod parsing; -#[cfg(feature = "serde")] +#[cfg(feature = "qc")] mod qc; #[cfg(feature = "serde")] diff --git a/src/record/mod.rs b/src/record/mod.rs index dc47da9..224cd2f 100644 --- a/src/record/mod.rs +++ b/src/record/mod.rs @@ -1,6 +1,9 @@ mod formatting; mod parsing; +#[cfg(feature = "qc")] +mod qc; + use std::collections::{btree_map::Iter, BTreeMap}; use itertools::Itertools; diff --git a/src/record/qc.rs b/src/record/qc.rs new file mode 100644 index 0000000..bf3471d --- /dev/null +++ b/src/record/qc.rs @@ -0,0 +1,31 @@ +use crate::prelude::Record; +use gnss_qc_traits::{Merge, MergeError}; + +impl Merge for Record { + fn merge(&self, rhs: &Self) -> Result { + let mut s = self.clone(); + s.merge_mut(rhs)?; + Ok(s) + } + + fn merge_mut(&mut self, rhs: &Self) -> Result<(), MergeError> { + for (rhs_k, rhs_v) in rhs.map.iter() { + if let Some(lhs_v) = self.map.get_mut(&rhs_k) { + if let Some(rhs_rms) = rhs_v.rms { + if lhs_v.rms.is_none() { + lhs_v.rms = Some(rhs_rms); + } + } + if let Some(rhs_height) = rhs_v.height { + if lhs_v.height.is_none() { + lhs_v.height = Some(rhs_height); + } + } + } else { + self.map.insert(*rhs_k, *rhs_v); + } + } + + Ok(()) + } +} diff --git a/src/tec.rs b/src/tec.rs index ce9f920..bf2ee60 100644 --- a/src/tec.rs +++ b/src/tec.rs @@ -14,7 +14,7 @@ pub struct TEC { pub(crate) rms: Option, /// Altitude offset for complex 3D height map - height: Option, + pub(crate) height: Option, } impl std::ops::Mul for TEC { From 9ced9685d9108f37e50c460df94d41595a403414 Mon Sep 17 00:00:00 2001 From: "Guillaume W. Bres" Date: Sat, 1 Nov 2025 12:27:01 +0100 Subject: [PATCH 07/25] file merging Signed-off-by: Guillaume W. Bres --- .gitignore | 9 ++---- src/lib.rs | 43 ++++++++++++++++++++++++++++ src/production.rs | 27 ++++++++++++++++++ src/tests/mod.rs | 1 + src/tests/qc.rs | 72 +++++++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 145 insertions(+), 7 deletions(-) create mode 100644 src/tests/qc.rs diff --git a/.gitignore b/.gitignore index 2a831cd..35518ea 100644 --- a/.gitignore +++ b/.gitignore @@ -25,10 +25,5 @@ target *.swo *.swp - -ckmg-v1.txt -jplg-v1.txt -custom.txt -region.txt -test.txt -test.txt.gz +*.txt +*.txt.gz diff --git a/src/lib.rs b/src/lib.rs index 1edc1d7..3281d11 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -844,6 +844,49 @@ impl IONEX { } } +/// Merge two [IONEX] structures into one. +/// This requires a few mandatory steps: +/// - reference systems must match +/// - maps dimension must match +/// - both must use the same mapping function +/// +/// Different sampling rate are supported, because the IONEX +/// description allows to describe that, but you will windup with +/// non constant sample rates. +#[cfg(feature = "qc")] +impl gnss_qc_traits::Merge for IONEX { + fn merge(&self, rhs: &Self) -> Result { + let mut s = self.clone(); + s.merge_mut(rhs)?; + Ok(s) + } + + fn merge_mut(&mut self, rhs: &Self) -> Result<(), gnss_qc_traits::MergeError> { + self.header.merge_mut(&rhs.header)?; + self.record.merge_mut(&rhs.record)?; + + match self.production { + Some(ref mut prods) => { + if let Some(rhs) = &rhs.production { + prods.merge_mut(rhs)?; + } + }, + None => { + if let Some(rhs) = &rhs.production { + self.production = Some(rhs.clone()); + } + }, + } + + // add new comments + for comment in rhs.comments.iter() { + self.comments.push(comment.clone()); + } + + Ok(()) + } +} + #[cfg(test)] mod test { use crate::{div_ceil, fmt_comment, prelude::*}; diff --git a/src/production.rs b/src/production.rs index 3848337..eeba49d 100644 --- a/src/production.rs +++ b/src/production.rs @@ -127,6 +127,33 @@ impl std::str::FromStr for ProductionAttributes { } } +#[cfg(feature = "qc")] +impl gnss_qc_traits::Merge for ProductionAttributes { + fn merge(&self, rhs: &Self) -> Result { + let mut s = self.clone(); + s.merge_mut(rhs)?; + Ok(s) + } + + fn merge_mut(&mut self, rhs: &Self) -> Result<(), gnss_qc_traits::MergeError> { + if rhs.year < self.year { + self.year = rhs.year; + self.doy = rhs.doy; + } + + if self.region == Region::Regional && rhs.region == Region::Global { + self.region = Region::Global; + } + + #[cfg(feature = "flate2")] + if self.gzip_compressed && !rhs.gzip_compressed { + self.gzip_compressed = false; + } + + Ok(()) + } +} + #[cfg(test)] mod test { use super::*; diff --git a/src/tests/mod.rs b/src/tests/mod.rs index 7baf321..becba47 100644 --- a/src/tests/mod.rs +++ b/src/tests/mod.rs @@ -3,6 +3,7 @@ pub mod toolkit; mod filename; mod parsing; +mod qc; mod v1; diff --git a/src/tests/qc.rs b/src/tests/qc.rs new file mode 100644 index 0000000..75f2d4d --- /dev/null +++ b/src/tests/qc.rs @@ -0,0 +1,72 @@ +use crate::{ + prelude::{coord, Duration, MappingFunction, Rect, Version, IONEX}, + tests::{ + init_logger, + toolkit::{generic_comparison, generic_test, TestPoint}, + }, +}; + +use gnss_qc_traits::Merge; + +use std::fs::File; +use std::io::BufWriter; + +#[test] +fn v1_self_merge() { + init_logger(); + + let ionex = IONEX::from_gzip_file("data/IONEX/V1/CKMG0020.22I.gz").unwrap_or_else(|e| { + panic!("Failed to parse CKMG0020: {}", e); + }); + + let merged = ionex.merge(&ionex).unwrap_or_else(|e| { + panic!("self::merge(self) should be feasible! {}", e); + }); + + // reciprocal + let fd = File::create("merged.txt").unwrap(); + let mut writer = BufWriter::new(fd); + + merged.format(&mut writer).unwrap_or_else(|e| { + panic!("failed to format V1 file: {}", e); + }); + + // parse back + let parsed = IONEX::from_file("merged.txt").unwrap_or_else(|e| { + panic!("failed to parse back CKMG V1: {}", e); + }); + + // full reciprocity + generic_comparison(&parsed, &ionex); +} + +#[test] +fn v1_files_merge() { + init_logger(); + + let file_a = IONEX::from_gzip_file("data/IONEX/V1/CKMG0090.21I.gz").unwrap_or_else(|e| { + panic!("Failed to parse V1 file: {}", e); + }); + let file_b = IONEX::from_gzip_file("data/IONEX/V1/CKMG0020.22I.gz").unwrap_or_else(|e| { + panic!("Failed to parse V1 file: {}", e); + }); + + let merged = file_a.merge(&file_b).unwrap_or_else(|e| { + panic!("failed to merge both files: {}", e); + }); + + // verifications + let fd = File::create("merged.txt").unwrap(); + let mut writer = BufWriter::new(fd); + + merged.format(&mut writer).unwrap_or_else(|e| { + panic!("failed to format V1 file: {}", e); + }); + + // parse back + let parsed = IONEX::from_file("ckmg-v1.txt").unwrap_or_else(|e| { + panic!("failed to parse back CKMG V1: {}", e); + }); + + // TODO +} From a1ff5371c8aca16093695cbd4aefc7c639fe4fb3 Mon Sep 17 00:00:00 2001 From: "Guillaume W. Bres" Date: Sat, 1 Nov 2025 12:28:04 +0100 Subject: [PATCH 08/25] file merging Signed-off-by: Guillaume W. Bres --- src/header/qc.rs | 4 +++- src/lib.rs | 4 +++- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/src/header/qc.rs b/src/header/qc.rs index e1bf38e..70056ce 100644 --- a/src/header/qc.rs +++ b/src/header/qc.rs @@ -80,7 +80,9 @@ impl Merge for Header { // TODO: merge grid def for comment in rhs.comments.iter() { - self.comments.push(comment.clone()); + if !self.comments.contains(&comment) { + self.comments.push(comment.clone()); + } } Ok(()) diff --git a/src/lib.rs b/src/lib.rs index 3281d11..858ea85 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -880,7 +880,9 @@ impl gnss_qc_traits::Merge for IONEX { // add new comments for comment in rhs.comments.iter() { - self.comments.push(comment.clone()); + if !self.comments.contains(&comment) { + self.comments.push(comment.clone()); + } } Ok(()) From 37b11f93b1b0788f271cd8f7a7189a13e0006ba5 Mon Sep 17 00:00:00 2001 From: "Guillaume W. Bres" Date: Sat, 1 Nov 2025 12:29:36 +0100 Subject: [PATCH 09/25] remove deprecated stuff Signed-off-by: Guillaume W. Bres --- src/rinex/feature.rs | 87 -------------------------------------------- src/rinex/mod.rs | 58 ----------------------------- 2 files changed, 145 deletions(-) delete mode 100644 src/rinex/feature.rs delete mode 100644 src/rinex/mod.rs diff --git a/src/rinex/feature.rs b/src/rinex/feature.rs deleted file mode 100644 index cf86641..0000000 --- a/src/rinex/feature.rs +++ /dev/null @@ -1,87 +0,0 @@ -use crate::{ - ionex::{IonexKey, TEC}, - prelude::{Epoch, Rinex}, -}; - -impl Rinex { - /// Returns Iterator over TEC expressed in TECu (10^-16 m-2) per lattitude, - /// longitude in decimal degrees, and altitude in km. - pub fn ionex_tecu_latlong_ddeg_alt_km_iter( - &self, - ) -> Box + '_> { - Box::new(self.ionex_tec_maps_iter().map(|(k, v)| { - ( - k.epoch, - k.coordinates.latitude_ddeg(), - k.coordinates.longitude_ddeg(), - k.coordinates.altitude_km(), - v.tecu(), - ) - })) - } - - /// Returns IONEX TEC map borders. - /// NB: this is only based on Header definitions. If provided - /// content did not follow those specs (incorrect), the returned value here will not - /// reflect actual content. - /// ## Output - /// - (lat_min, lat_max): southernmost and northernmost latitude, in decimal degrees - /// - (long_min, long_max): easternmost and westernmost longitude, in decimal degrees - pub fn tec_map_borders(&self) -> Option<((f64, f64), (f64, f64))> { - let ionex = self.header.ionex.as_ref()?; - Some(( - (ionex.grid.latitude.start, ionex.grid.latitude.end), - (ionex.grid.longitude.start, ionex.grid.longitude.end), - )) - } - - /// Designs an iterator over RMS TEC exclusively - pub fn ionex_rms_tec_maps_iter(&self) -> Box + '_> { - Box::new(self.ionex_tec_maps_iter().filter_map(|(k, tec)| { - if let Some(rms) = tec.rms_tec() { - Some((k, rms)) - } else { - None - } - })) - } - - /// Returns fixed altitude (in kilometers) if this is a 2D IONEX (only). - /// Note that this is only based on Header definitions. If provided - /// content does not follow those specs (incorrect data), the returned value here will not - /// reflect actual content. - pub fn ionex_2d_fixed_altitude_km(&self) -> Option { - if self.is_ionex_2d() { - let header = self.header.ionex.as_ref()?; - Some(header.grid.height.start) - } else { - None - } - } - - /// Returns altitude range of the IONEX TEC maps, expressed as {min, max} - /// both in kilometers. - /// - if this is a 2D IONEX: you will obtain (min, min) - /// - if this is a 3D IONEX: you will obtain (min, max) - /// Note that this is only based on Header definitions. If provided - /// content does not follow those specs (incorrect data), the returned value here will not - /// reflect actual content. - pub fn ionex_altitude_range_km(&self) -> Option<(f64, f64)> { - let header = self.header.ionex.as_ref()?; - Some((header.grid.height.start, header.grid.height.end)) - } - - /// Designs an TEC isosurface iterator starting at lowest altitude, - /// ending at highest altitude. See [Self::ionex_tec_maps_altitude_range] for - /// useful information. - pub fn ionex_tec_isosurface_iter(&self) -> Box + '_> { - Box::new([].into_iter()) - } - - /// Designs an RMS TEC isosurface iterator starting at lowest altitude, - /// ending at highest altitude. See [Self::ionex_tec_maps_altitude_range] for - /// useful information. - pub fn ionex_rms_tec_isosurface_iter(&self) -> Box + '_> { - Box::new([].into_iter()) - } -} diff --git a/src/rinex/mod.rs b/src/rinex/mod.rs deleted file mode 100644 index cfed3cd..0000000 --- a/src/rinex/mod.rs +++ /dev/null @@ -1,58 +0,0 @@ -#[cfg(feature = "ionex")] -#[cfg_attr(docsrs, doc(cfg(feature = "ionex")))] -mod feature; - -use crate::{ - ionex::{IonexKey, TEC}, - prelude::{Rinex, RinexType}, -}; - -use std::collections::btree_map::Keys; - -impl Rinex { - /// Returns true if this is a IONEX (special) [Rinex] - pub fn is_ionex(&self) -> bool { - self.header.rinex_type == RinexType::IonosphereMaps - } - - /// Returns true if this IONEX only contains a single isosurface - /// at fixed altitude. NB: this information only relies - /// on the [Header] section, not actual data analysis. - /// If [Record] content did not follow the specifications, this will be invalid. - pub fn is_ionex_2d(&self) -> bool { - if let Some(ionex) = &self.header.ionex { - ionex.map_dimension == 2 && self.is_ionex() - } else { - false - } - } - - /// Returns true if this IONEX contains several isosurface spanning [Self::ionex_altitude_range]. - /// NB: this information only relies - /// on the [Header] section, not actual data analysis. - /// If [Record] content did not follow the specifications, this will be invalid. - pub fn is_ionex_3d(&self) -> bool { - if let Some(ionex) = &self.header.ionex { - ionex.map_dimension == 3 && self.is_ionex() - } else { - false - } - } - - pub fn ionex_tec_maps_keys(&self) -> Keys<'_, IonexKey, TEC> { - if let Some(rec) = self.record.as_ionex() { - rec.keys() - } else { - panic!("invalid operation"); - } - } - - /// IONEX Total Electron Content Iterator. - pub fn ionex_tec_maps_iter(&self) -> Box + '_> { - if let Some(rec) = self.record.as_ionex() { - Box::new(rec.iter().map(|(k, v)| (*k, v))) - } else { - Box::new([].into_iter()) - } - } -} From ed49feb60842efc49a83eaab81bb0202e349d805 Mon Sep 17 00:00:00 2001 From: "Guillaume W. Bres" Date: Sun, 2 Nov 2025 10:11:10 +0100 Subject: [PATCH 10/25] update API Signed-off-by: Guillaume W. Bres --- README.md | 4 +- src/{production.rs => file_attributes.rs} | 43 ++++--- src/lib.rs | 132 ++++++++++++---------- src/tests/filename.rs | 2 +- 4 files changed, 98 insertions(+), 83 deletions(-) rename src/{production.rs => file_attributes.rs} (82%) diff --git a/README.md b/README.md index 3f29330..68f719e 100644 --- a/README.md +++ b/README.md @@ -96,14 +96,14 @@ assert_eq!(ionex.header.base_radius_km, 6371.0); assert!(ionex.is_2d()); // this file is named according to IGS standards -let descriptor = ionex.production.clone().unwrap(); +let descriptor = ionex.attributes.clone().unwrap(); // to obtain TEC values at any coordinates, you // should use the [MapCell] local region (rectangle quanta) // that offers many functions based off the Geo crate. // Convenient helper to follow standard conventions -let filename = ionex.standardized_filename(); +let filename = ionex.generate_standardized_filename(); // Dump to file let fd = File::create("custom.txt").unwrap(); diff --git a/src/production.rs b/src/file_attributes.rs similarity index 82% rename from src/production.rs rename to src/file_attributes.rs index eeba49d..a73d901 100644 --- a/src/production.rs +++ b/src/file_attributes.rs @@ -1,7 +1,6 @@ -//! File production infrastructure. use thiserror::Error; -/// File Production identification errors +/// File info parsing errors #[derive(Error, Debug)] pub enum Error { #[error("filename does not follow naming conventions")] @@ -10,19 +9,19 @@ pub enum Error { #[derive(Debug, Copy, Default, Clone, PartialEq)] pub enum Region { - /// Local IONEX (Regional maps) + /// Local/Regional IONEX map (specific ROI). Regional, - /// Global IONEX (Worldwide maps) + /// Worldwide IONEX map #[default] - Global, + WorldWide, } impl std::fmt::Display for Region { fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result { match self { Self::Regional => write!(f, "{}", 'R'), - Self::Global => write!(f, "{}", 'G'), + Self::WorldWide => write!(f, "{}", 'G'), } } } @@ -31,14 +30,14 @@ impl std::fmt::Display for Region { /// RINEX data that follows standard naming conventions, /// or attached to data parsed from such files. #[derive(Debug, Default, Clone, PartialEq)] -pub struct ProductionAttributes { - /// Production agency +pub struct FileAttributes { + /// File agency pub agency: String, /// Year of production pub year: u32, - /// Production Day of Year (DOY). + /// File Day of Year (DOY). /// We assume past J2000. pub doy: u32, @@ -51,7 +50,7 @@ pub struct ProductionAttributes { pub gzip_compressed: bool, } -impl std::fmt::Display for ProductionAttributes { +impl std::fmt::Display for FileAttributes { #[cfg(feature = "flate2")] fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result { let len = std::cmp::min(self.agency.len(), 3); @@ -86,7 +85,7 @@ impl std::fmt::Display for ProductionAttributes { } } -impl std::str::FromStr for ProductionAttributes { +impl std::str::FromStr for FileAttributes { type Err = Error; fn from_str(filename: &str) -> Result { @@ -107,7 +106,7 @@ impl std::str::FromStr for ProductionAttributes { .map_err(|_| Error::NonStandardFilename)?; let region = if filename[3..4].eq("G") { - Region::Global + Region::WorldWide } else { Region::Regional }; @@ -128,7 +127,7 @@ impl std::str::FromStr for ProductionAttributes { } #[cfg(feature = "qc")] -impl gnss_qc_traits::Merge for ProductionAttributes { +impl gnss_qc_traits::Merge for FileAttributes { fn merge(&self, rhs: &Self) -> Result { let mut s = self.clone(); s.merge_mut(rhs)?; @@ -141,8 +140,8 @@ impl gnss_qc_traits::Merge for ProductionAttributes { self.doy = rhs.doy; } - if self.region == Region::Regional && rhs.region == Region::Global { - self.region = Region::Global; + if self.region == Region::Regional && rhs.region == Region::WorldWide { + self.region = Region::WorldWide; } #[cfg(feature = "flate2")] @@ -162,15 +161,15 @@ mod test { #[test] fn standard_filenames() { for (filename, agency, year, doy, region) in [ - ("CKMG0020.22I", "CKM", 2022, 2, Region::Global), - ("CKMG0090.21I", "CKM", 2021, 9, Region::Global), - ("JPLG0010.17I", "JPL", 2017, 1, Region::Global), + ("CKMG0020.22I", "CKM", 2022, 2, Region::WorldWide), + ("CKMG0090.21I", "CKM", 2021, 9, Region::WorldWide), + ("JPLG0010.17I", "JPL", 2017, 1, Region::WorldWide), ("JPLR0010.17I", "JPL", 2017, 1, Region::Regional), ("JPLR0010.17I", "JPL", 2017, 1, Region::Regional), ] { println!("Testing IONEX filename \"{}\"", filename); - let attrs = ProductionAttributes::from_str(filename).unwrap(); + let attrs = FileAttributes::from_str(filename).unwrap(); assert_eq!(attrs.agency, agency); assert_eq!(attrs.year, year); @@ -186,12 +185,12 @@ mod test { #[test] fn gzip_filenames() { for (filename, agency, year, doy, region) in [ - ("CKMG0020.22I.gz", "CKM", 2022, 2, Region::Global), + ("CKMG0020.22I.gz", "CKM", 2022, 2, Region::WorldWide), ("CKMR0020.22I.gz", "CKM", 2022, 2, Region::Regional), - ("JPLG0010.17I.gz", "JPL", 2017, 1, Region::Global), + ("JPLG0010.17I.gz", "JPL", 2017, 1, Region::WorldWide), ("CKMR0020.22I.gz", "CKM", 2022, 2, Region::Regional), ] { - let attrs = ProductionAttributes::from_str(filename).unwrap(); + let attrs = FileAttributes::from_str(filename).unwrap(); assert_eq!(attrs.agency, agency); assert_eq!(attrs.year, year); diff --git a/src/lib.rs b/src/lib.rs index 858ea85..1aac8d2 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -23,12 +23,12 @@ extern crate gnss_rs as gnss; pub mod bias; pub mod error; +pub mod file_attributes; pub mod grid; pub mod header; pub mod key; pub mod linspace; pub mod mapf; -pub mod production; pub mod system; pub mod tec; pub mod version; @@ -64,9 +64,9 @@ use crate::{ cell::{MapCell, MapPoint}, coordinates::QuantizedCoordinates, error::{FormattingError, ParsingError}, + file_attributes::{FileAttributes, Region}, header::Header, key::Key, - production::{ProductionAttributes, Region}, quantized::Quantized, record::Record, tec::TEC, @@ -78,13 +78,13 @@ pub mod prelude { bias::BiasSource, cell::MapCell, error::{FormattingError, ParsingError}, + file_attributes::*, grid::Grid, header::Header, ionosphere::IonosphereParameters, key::Key, linspace::Linspace, mapf::MappingFunction, - production::*, record::Record, system::ReferenceSystem, tec::TEC, @@ -152,7 +152,7 @@ pub(crate) fn fmt_comment(content: &str) -> String { /// /// use ionex::prelude::*; /// -/// // Parse Global/worldwide map +/// // Worldwide (so called 'global') IONEX map /// let ionex = IONEX::from_gzip_file("data/IONEX/V1/CKMG0020.22I.gz") /// .unwrap(); /// @@ -174,14 +174,14 @@ pub(crate) fn fmt_comment(content: &str) -> String { /// assert!(ionex.is_2d()); /// /// // this file is named according to IGS standards -/// let descriptor = ionex.production.clone().unwrap(); +/// let descriptor = ionex.attributes.clone().unwrap(); /// /// // to obtain TEC values at any coordinates, you /// // should use the [MapCell] local region (rectangle quanta) /// // that offers many functions based off the Geo crate. /// /// // Convenient helper to follow standard conventions -/// let filename = ionex.standardized_filename(); +/// let filename = ionex.generate_standardized_filename(); /// /// // Dump to file /// let fd = File::create("custom.txt").unwrap(); @@ -209,9 +209,8 @@ pub struct IONEX { /// [Comments] found in file record pub comments: Comments, - /// [ProductionAttributes] resolved for file names that follow - /// according to the standards. - pub production: Option, + /// [FileAttributes] resolved for file names that follow the IGS conventions. + pub attributes: Option, } impl IONEX { @@ -220,7 +219,7 @@ impl IONEX { Self { header, record, - production: None, + attributes: None, comments: Default::default(), } } @@ -231,7 +230,7 @@ impl IONEX { header, record: self.record.clone(), comments: self.comments.clone(), - production: self.production.clone(), + attributes: self.attributes.clone(), } } @@ -246,7 +245,7 @@ impl IONEX { record, header: self.header.clone(), comments: self.comments.clone(), - production: self.production.clone(), + attributes: self.attributes.clone(), } } @@ -272,21 +271,21 @@ impl IONEX { /// Returns a file name that would describe [Self] according to the /// standards. - pub fn standardized_filename(&self) -> String { - let (agency, region, year, doy) = if let Some(production) = &self.production { + pub fn generate_standardized_filename(&self) -> String { + let (agency, region, year, doy) = if let Some(attributes) = &self.attributes { ( - production.agency.clone(), - production.region, - production.year - 2000, - production.doy, + attributes.agency.clone(), + attributes.region, + attributes.year - 2000, + attributes.doy, ) } else { ("XXX".to_string(), Region::default(), 0, 0) }; - let extension = if let Some(production) = &self.production { + let extension = if let Some(attributes) = &self.attributes { #[cfg(feature = "flate2")] - if production.gzip_compressed { + if attributes.gzip_compressed { ".gz" } else { "" @@ -298,12 +297,13 @@ impl IONEX { format!("{}{}{:03}0.{:02}I{}", agency, region, doy, year, extension) } - /// Guesses File [ProductionAttributes] from actual record content. - /// This is useful to generate a standardized file name, from good data content - /// parsed from files that do not follow the standard naming conventions. - /// The agency is only determined by the file name so you should provide one, - /// and it must be a at least 3 letter code. - pub fn guess_production_attributes(&self, agency: &str) -> Option { + /// Guesses [FileAttributes] from actual dataset. This is particularly useful + /// to generate a standardized file name, especially when arriving from data that + /// did not follow the conventions. + /// The name of the production agency (data provider) is only determined by the [FileAttributes] + /// itself, so we have no means to generate correct one, so we need you to define it right here. + /// The production agency should be at least a 3 letter code, for example: "IGS". + pub fn guess_file_attributes(&self, agency: &str) -> Option { if agency.len() < 3 { return None; } @@ -312,16 +312,16 @@ impl IONEX { let year = first_epoch.year(); let doy = first_epoch.day_of_year().round() as u32; - let region = Region::Global; // TODO: study the grid specs + let region = Region::WorldWide; // TODO: study the grid specs - Some(ProductionAttributes { + Some(FileAttributes { doy, region, year: year as u32, agency: agency.to_string(), #[cfg(feature = "flate2")] - gzip_compressed: if let Some(attributes) = &self.production { + gzip_compressed: if let Some(attributes) = &self.attributes { attributes.gzip_compressed } else { false @@ -344,7 +344,7 @@ impl IONEX { header, record, comments, - production: Default::default(), + attributes: Default::default(), }) } @@ -371,12 +371,12 @@ impl IONEX { /// See [Self::from_gzip_file] for seamless Gzip support. /// /// If file name follows standard naming conventions, then internal definitions - /// will truly be complete. Otherwise [ProductionAttributes] cannot be fully determined. + /// will truly be complete. Otherwise [FileAttributes] cannot be fully determined. /// If you want or need to you can either /// 1. define it yourself with further customization - /// 2. use the smart guesser (after parsing): [Self::guess_production_attributes] + /// 2. use the smart guesser (after parsing): [Self::guess_attributes_attributes] /// - /// This is typically needed in data production contexts. + /// This is typically needed in data attributes contexts. pub fn from_file>(path: P) -> Result { let path = path.as_ref(); @@ -384,7 +384,7 @@ impl IONEX { let file_attributes = match path.file_name() { Some(filename) => { let filename = filename.to_string_lossy().to_string(); - if let Ok(prod) = ProductionAttributes::from_str(&filename) { + if let Ok(prod) = FileAttributes::from_str(&filename) { Some(prod) } else { None @@ -397,7 +397,7 @@ impl IONEX { let mut reader = BufReader::new(fd); let mut ionex = Self::parse(&mut reader)?; - ionex.production = file_attributes; + ionex.attributes = file_attributes; Ok(ionex) } @@ -416,9 +416,9 @@ impl IONEX { /// assert!(ionex.to_file("test.txt").is_ok()); /// ``` /// - /// Other useful links are in data production contexts: + /// Other useful links are in data attributes contexts: /// * [Self::standard_filename] to generate a standardized filename - /// * [Self::guess_production_attributes] helps generate standardized filenames for + /// * [Self::guess_attributes_attributes] helps generate standardized filenames for /// files that do not follow naming conventions pub fn to_file>(&self, path: P) -> Result<(), FormattingError> { let fd = File::create(path)?; @@ -466,7 +466,7 @@ impl IONEX { let file_attributes = match path.file_name() { Some(filename) => { let filename = filename.to_string_lossy().to_string(); - if let Ok(prod) = ProductionAttributes::from_str(&filename) { + if let Ok(prod) = FileAttributes::from_str(&filename) { Some(prod) } else { None @@ -481,14 +481,15 @@ impl IONEX { let mut reader = BufReader::new(reader); let mut ionex = Self::parse(&mut reader)?; - ionex.production = file_attributes; + ionex.attributes = file_attributes; Ok(ionex) } /// Dumps and gzip encodes [IONEX] into writable local file, /// using efficient buffered formatting. - /// This is the mirror operation of [Self::from_gzip_file]. + /// + /// This operation is [Self::from_gzip_file] mirror operation. /// ``` /// // Read a IONEX and dump it without any modifications /// use ionex::prelude::*; @@ -499,10 +500,14 @@ impl IONEX { /// assert!(ionex.to_gzip_file("test.txt.gz").is_ok()); /// ``` /// - /// Other useful links are in data production contexts: - /// * [Self::standard_filename] to generate a standardized filename - /// * [Self::guess_production_attributes] helps generate standardized filenames for - /// files that do not follow naming conventions + /// Other useful methods are: + /// * [Self::generate_standardized_filename]: to generate a standardize file name + /// * [Self::guess_file_attributes] in close relation: helps generate standardized + /// filenames for files that did now follow the convention + /// * [gnss_qc_traits::Merge]: to combine two files to one another + /// * [Self::to_regional_roi]: to reduce a larger (possibly [Region::GlobalWorldWide] map) + /// to a given ROI + /// * [Self::to_worldwide_map]: to enlarge a possibly [Region::Regional] map to a [Region::GlobalWorldWide] IONEX. #[cfg(feature = "flate2")] #[cfg_attr(docsrs, doc(cfg(feature = "flate2")))] pub fn to_gzip_file>(&self, path: P) -> Result<(), FormattingError> { @@ -539,8 +544,8 @@ impl IONEX { pub fn to_worldwide_ionex(&self) -> IONEX { let mut ionex = self.clone(); - if let Some(production) = &mut ionex.production { - production.region = Region::Global; + if let Some(attributes) = &mut ionex.attributes { + attributes.region = Region::WorldWide; } // update grid specs, preserve accuracy @@ -554,23 +559,34 @@ impl IONEX { ionex } + /// Returns true if this [IONEX] is a Worldwide map + /// (as opposed to a local/regional ROI). + pub fn is_worldwide_map(&self) -> bool { + if let Some(attributes) = &self.attributes { + attributes.region == Region::WorldWide + } else { + let bounding_rect = self.map_borders_degrees(); + bounding_rect.width() == 360.0 && bounding_rect.height() == 175.0 + } + } + /// Converts this Global/Worldwide [IONEX] to Regional [IONEX] /// delimited by (possibly complex) bounding geometry, for which [BoundingRect] needs /// to return correctly. Each boundary coordinates must also be expressed in degrees. /// NB: work in progress, not validated yet. - pub fn to_regional_ionex(&self, region: Polygon) -> Option { + pub fn to_regional_roi(&self, roi: Polygon) -> Option { let mut ionex = IONEX::default(); - let bounding_rect = region.bounding_rect()?; + let bounding_rect = roi.bounding_rect()?; let (min_long, min_lat) = (bounding_rect.min().x, bounding_rect.min().y); let (max_long, max_lat) = (bounding_rect.max().x, bounding_rect.max().y); // copy attributes - ionex.production = self.production.clone(); + ionex.attributes = self.attributes.clone(); - if let Some(production) = &mut ionex.production { - production.region = Region::Regional; + if let Some(attributes) = &mut ionex.attributes { + attributes.region = Region::Regional; } // copy & rework header @@ -581,12 +597,12 @@ impl IONEX { ionex.header.grid.longitude.start = min_long; ionex.header.grid.longitude.end = max_long; - let region = Geometry::Polygon(region); + let roi = Geometry::Polygon(roi); // restrict area let cells = self .map_cell_iter() - .filter(|cell| cell.contains(®ion)) + .filter(|cell| cell.contains(&roi)) .collect::>(); ionex.record = Record::from_map_cells( @@ -865,15 +881,15 @@ impl gnss_qc_traits::Merge for IONEX { self.header.merge_mut(&rhs.header)?; self.record.merge_mut(&rhs.record)?; - match self.production { + match self.attributes { Some(ref mut prods) => { - if let Some(rhs) = &rhs.production { + if let Some(rhs) = &rhs.attributes { prods.merge_mut(rhs)?; } }, None => { - if let Some(rhs) = &rhs.production { - self.production = Some(rhs.clone()); + if let Some(rhs) = &rhs.attributes { + self.attributes = Some(rhs.clone()); } }, } @@ -940,7 +956,7 @@ mod test { let roi = Rect::new(coord!(x: -180.0, y: -85.0), coord!(x: 180.0, y: -82.5)); - let regional = ionex.to_regional_ionex(roi.into()).unwrap(); + let regional = ionex.to_regional_roi(roi.into()).unwrap(); // dump regional.to_file("region.txt").unwrap(); diff --git a/src/tests/filename.rs b/src/tests/filename.rs index 23d3dda..ccb37ab 100644 --- a/src/tests/filename.rs +++ b/src/tests/filename.rs @@ -5,6 +5,6 @@ fn filename_conventions() { for testfile in ["CKMG0020.22I.gz", "CKMG0080.09I.gz", "CKMG0090.21I.gz"] { let name = format!("data/IONEX/V1/{}", testfile); let ionex = IONEX::from_gzip_file(&name).unwrap(); - assert_eq!(ionex.standardized_filename(), testfile); + assert_eq!(ionex.generate_standardized_filename(), testfile); } } From a68bb254f6d267ba8fc87e4e8e574cab262a7823 Mon Sep 17 00:00:00 2001 From: "Guillaume W. Bres" Date: Sun, 2 Nov 2025 10:11:37 +0100 Subject: [PATCH 11/25] update API Signed-off-by: Guillaume W. Bres --- src/lib.rs | 1 - 1 file changed, 1 deletion(-) diff --git a/src/lib.rs b/src/lib.rs index 1aac8d2..d265e78 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -3,7 +3,6 @@ )] #![doc = include_str!("../README.md")] #![cfg_attr(docsrs, feature(doc_cfg))] -#![allow(clippy::type_complexity)] /* * IONEX is part of the nav-solutions framework. From 35633612ece2740ac2f47fca7e82e4f1615dddd9 Mon Sep 17 00:00:00 2001 From: "Guillaume W. Bres" Date: Sun, 2 Nov 2025 11:03:13 +0100 Subject: [PATCH 12/25] linspace stretching Signed-off-by: Guillaume W. Bres --- src/error.rs | 6 ++++ src/lib.rs | 57 +++++++++++++++++-------------- src/linspace.rs | 89 +++++++++++++++++++++++++++++++++++++++++++++++- src/tests/mod.rs | 1 + src/tests/roi.rs | 36 ++++++++++++++++++++ src/tests/v1.rs | 2 +- 6 files changed, 164 insertions(+), 27 deletions(-) create mode 100644 src/tests/roi.rs diff --git a/src/error.rs b/src/error.rs index 1f92308..1ab0faf 100644 --- a/src/error.rs +++ b/src/error.rs @@ -91,6 +91,12 @@ pub enum ParsingError { ExponentScaling, } +#[derive(Error, Debug)] +pub enum Error { + #[error("strech factor must be positive")] + NegativeStretchFactor, +} + /// Errors that may rise during Formatting process #[derive(Error, Debug)] pub enum FormattingError { diff --git a/src/lib.rs b/src/lib.rs index d265e78..bded9d6 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -76,7 +76,7 @@ pub mod prelude { pub use crate::{ bias::BiasSource, cell::MapCell, - error::{FormattingError, ParsingError}, + error::{FormattingError, ParsingError, Error}, file_attributes::*, grid::Grid, header::Header, @@ -529,17 +529,33 @@ impl IONEX { false } - /// Returns map borders as [Rect]angle. This uses - /// the [Header] description and assumes all maps are within these borders. - pub fn map_borders_degrees(&self) -> Rect { + /// Returns map borders as a [Rect]angle, with coordinates in decimal degrees. + /// This uses the [Header] description and assumes all maps are within these borders. + pub fn bounding_rect_degrees(&self) -> Rect { Rect::new( coord!( x: self.header.grid.longitude.start, y: self.header.grid.latitude.start ), coord!( x: self.header.grid.longitude.end, y: self.header.grid.latitude.end), ) } - /// Stretch this Regional [IONEX] to a Global/Worldwide [IONEX]. - /// NB: work in progress, not validated yet. + /// Returns true if this [IONEX] is a Worldwide map + /// (as opposed to a local/regional ROI). + pub fn is_worldwide_map(&self) -> bool { + if let Some(attributes) = &self.attributes { + attributes.region == Region::WorldWide + } else { + let bounding_rect = self.bounding_rect_degrees(); + bounding_rect.width() == 360.0 && bounding_rect.height() == 175.0 + } + } + + /// Returns true if this [IONEX] is not a Worldwide map but a regional/local ROI. + pub fn is_regional_map(&self) -> bool { + !self.is_worldwide_map() + } + + /// Stretch this [IONEX] definition so it becomes compatible + /// with the description of a Global/Worldwide [IONEX]. pub fn to_worldwide_ionex(&self) -> IONEX { let mut ionex = self.clone(); @@ -558,22 +574,10 @@ impl IONEX { ionex } - /// Returns true if this [IONEX] is a Worldwide map - /// (as opposed to a local/regional ROI). - pub fn is_worldwide_map(&self) -> bool { - if let Some(attributes) = &self.attributes { - attributes.region == Region::WorldWide - } else { - let bounding_rect = self.map_borders_degrees(); - bounding_rect.width() == 360.0 && bounding_rect.height() == 175.0 - } - } - - /// Converts this Global/Worldwide [IONEX] to Regional [IONEX] - /// delimited by (possibly complex) bounding geometry, for which [BoundingRect] needs - /// to return correctly. Each boundary coordinates must also be expressed in degrees. - /// NB: work in progress, not validated yet. - pub fn to_regional_roi(&self, roi: Polygon) -> Option { + /// Reduce this [IONEX] definition so it is reduced to a regional ROI, + /// described by a complex [Polygon] in decimal degrees. + /// [Polygon::bounding_rect] must be defined for this operation to work correctly. + pub fn to_regional_ionex(&self, roi: Polygon) -> Option { let mut ionex = IONEX::default(); let bounding_rect = roi.bounding_rect()?; @@ -588,9 +592,10 @@ impl IONEX { attributes.region = Region::Regional; } - // copy & rework header + // copy header ionex.header = self.header.clone(); + // rework header ionex.header.grid.latitude.start = max_lat; ionex.header.grid.latitude.end = min_lat; ionex.header.grid.longitude.start = min_long; @@ -773,7 +778,7 @@ impl IONEX { /// ]); /// /// // you can double check that - /// assert!(ionex.map_borders_degrees().contains(&roi_ddeg)); + /// assert!(ionex.bounding_rect_degrees().contains(&roi_ddeg)); /// /// // Epoch must exist in the record /// let noon = Epoch::from_str("2022-01-02T12:00:00 UTC") @@ -955,7 +960,8 @@ mod test { let roi = Rect::new(coord!(x: -180.0, y: -85.0), coord!(x: 180.0, y: -82.5)); - let regional = ionex.to_regional_roi(roi.into()).unwrap(); + let regional = ionex.to_regional_ionex(roi.into()) + .unwrap(); // dump regional.to_file("region.txt").unwrap(); @@ -965,4 +971,5 @@ mod test { panic!("Failed to parse region.txt: {}", e); }); } + } diff --git a/src/linspace.rs b/src/linspace.rs index 2a54218..96d73b8 100644 --- a/src/linspace.rs +++ b/src/linspace.rs @@ -1,6 +1,6 @@ use std::ops::Rem; -use crate::{error::ParsingError, quantized::Quantized}; +use crate::{error::{ParsingError, Error}, quantized::Quantized}; /// Quantized Linspace for iteration #[derive(Debug, Copy, Clone, Default, PartialEq, PartialOrd)] @@ -67,6 +67,51 @@ impl Linspace { (self.min(), self.max()) } + /// Stretch this mutable [Linspace] by a positive, possibly fractional number, + /// while preserving the initial grid quantization (point spacing). + /// This only modifies the [Linspace] dimensions. + /// To modify the sampling, use [Self::resample_mut]. + pub fn stretch_mut(&mut self, factor: f64) -> Result<(), Error> { + if factor.is_sign_negative() { + return Err(Error::NegativeStretchFactor); + } + self.start *= factor; + self.end *= factor; + Ok(()) + } + + /// Stretch this [Linspace] to a new [Linspace] definition. + /// The streetching factor must be a positive number. + /// this preserves the initial grid quantization (point spacing). + /// This only modifies the [Linspace] dimensions. + /// To modify the sampling, use [Self::resampled]. + pub fn stretched(&self, factor: f64) -> Result { + let mut s = self.clone(); + s.stretch_mut(factor)?; + Ok(s) + } + + /// Resample this mutable [Linspace] so the point spacing is modified. + /// This is a multiplicative stretching factor, which must be a positive number. + /// This does not modify the [Linspace] dimensions, use [Self::stretch_mut] to modify dimensions. + pub fn resample_mut(&mut self, factor: f64) -> Result<(), Error> { + if factor.is_sign_negative() { + return Err(Error::NegativeStretchFactor); + } + self.spacing *= factor; + Ok(()) + } + + /// Resample this [Linspace] to a new [Linspace] definition, modifying the point + /// spacing but preserving the initial dimensions. The resampling factor + /// must be positive, possibly fractional number. This does not modify the dimensions, + /// use [Self::stretched] to modify the dimensions. + pub fn resampled(&self, factor: f64) -> Result { + let mut s = self.clone(); + s.resample_mut(factor)?; + Ok(s) + } + /// Builds a new Linear space pub fn new(start: f64, end: f64, spacing: f64) -> Result { if start == end && spacing == 0.0 { @@ -189,4 +234,46 @@ mod test { assert_eq!(linspace.length(), 10); assert!(!linspace.is_single_point()); } + + #[test] + fn linspace_stretching() { + let mut linspace = Linspace::new(-180.0, 180.0, 5.0).unwrap(); + + linspace.stretch_mut(0.5).unwrap(); + assert_eq!(linspace.minmax(), (-90.0, 90.0)); + assert_eq!(linspace.spacing, 5.0, "linspace quantization not preserved!"); + + linspace.stretch_mut(0.75).unwrap(); + assert_eq!(linspace.minmax(), (-67.5, 67.5)); + assert_eq!(linspace.spacing, 5.0, "linspace quantization not preserved!"); + + linspace.stretch_mut(0.5).unwrap(); + assert_eq!(linspace.minmax(), (-33.75, 33.75)); + assert_eq!(linspace.spacing, 5.0, "linspace quantization not preserved!"); + + linspace.stretch_mut(2.0).unwrap(); + assert_eq!(linspace.minmax(), (-67.5, 67.5)); + assert_eq!(linspace.spacing, 5.0, "linspace quantization not preserved!"); + } + + #[test] + fn linspace_resampling() { + let mut linspace = Linspace::new(-180.0, 180.0, 5.0).unwrap(); + + linspace.resample_mut(0.5).unwrap(); + assert_eq!(linspace.spacing, 2.5); + assert_eq!(linspace.minmax(), (-180.0, 180.0), "dimensions not preserved!"); + + linspace.resample_mut(0.5).unwrap(); + assert_eq!(linspace.spacing, 1.25); + assert_eq!(linspace.minmax(), (-180.0, 180.0), "dimensions not preserved!"); + + linspace.resample_mut(2.0).unwrap(); + assert_eq!(linspace.spacing, 2.5); + assert_eq!(linspace.minmax(), (-180.0, 180.0), "dimensions not preserved!"); + + linspace.resample_mut(2.0).unwrap(); + assert_eq!(linspace.spacing, 5.0); + assert_eq!(linspace.minmax(), (-180.0, 180.0), "dimensions not preserved!"); + } } diff --git a/src/tests/mod.rs b/src/tests/mod.rs index becba47..0f1a3e6 100644 --- a/src/tests/mod.rs +++ b/src/tests/mod.rs @@ -4,6 +4,7 @@ pub mod toolkit; mod filename; mod parsing; mod qc; +mod roi; mod v1; diff --git a/src/tests/roi.rs b/src/tests/roi.rs new file mode 100644 index 0000000..3051c1c --- /dev/null +++ b/src/tests/roi.rs @@ -0,0 +1,36 @@ +use crate::{ + tests::init_logger, + prelude::{ + IONEX, + Rect, + coord, + }, +}; + +#[test] +fn worldwide_map() { + init_logger(); + + let mut ionex = IONEX::from_gzip_file("data/IONEX/V1/CKMG0020.22I.gz").unwrap_or_else(|e| { + panic!("Failed to parse CKMG0020: {}", e); + }); + + // returned from file name directly + assert!(ionex.is_worldwide_map(), "Worldwide map from standardized file name!"); + + // destroy and verify this still works + ionex.attributes = None; + assert!(ionex.is_worldwide_map(), "Guessed worldwide map !"); + + // reduce to ROI + let roi = Rect::new(coord!(x: -30.0, y: -30.0), coord!(x: 30.0, y: 30.0)); + + let reduced = ionex.to_regional_ionex(roi.into()) + .unwrap(); + + assert!(reduced.is_regional_map(), "Worldwide map reduced to Regional ROI"); + + // test new map + let bounding_rect = reduced.bounding_rect_degrees(); + assert_eq!(bounding_rect, roi); +} diff --git a/src/tests/v1.rs b/src/tests/v1.rs index da7c50f..c3e505b 100644 --- a/src/tests/v1.rs +++ b/src/tests/v1.rs @@ -220,7 +220,7 @@ fn parse_ckmg0020() { assert_eq!(ionex.header.exponent, -1); assert_eq!( - ionex.map_borders_degrees(), + ionex.bounding_rect_degrees(), Rect::new(coord!(x: -180.0, y: -87.5), coord!(x: 180.0, y: 87.5)) ); From 8676e58eb9f52e74da5183caa79477d2fe3d0f52 Mon Sep 17 00:00:00 2001 From: "Guillaume W. Bres" Date: Sun, 2 Nov 2025 11:08:22 +0100 Subject: [PATCH 13/25] linspace stretching Signed-off-by: Guillaume W. Bres --- src/error.rs | 4 ++-- src/linspace.rs | 22 +++++++++++----------- 2 files changed, 13 insertions(+), 13 deletions(-) diff --git a/src/error.rs b/src/error.rs index 1ab0faf..c9f97d3 100644 --- a/src/error.rs +++ b/src/error.rs @@ -93,8 +93,8 @@ pub enum ParsingError { #[derive(Error, Debug)] pub enum Error { - #[error("strech factor must be positive")] - NegativeStretchFactor, + #[error("strech factor must be positive finite number")] + InvalidStretchFactor, } /// Errors that may rise during Formatting process diff --git a/src/linspace.rs b/src/linspace.rs index 96d73b8..1b4f90c 100644 --- a/src/linspace.rs +++ b/src/linspace.rs @@ -67,13 +67,13 @@ impl Linspace { (self.min(), self.max()) } - /// Stretch this mutable [Linspace] by a positive, possibly fractional number, + /// Stretch this mutable [Linspace] by a finite positive number (possibly fractional), /// while preserving the initial grid quantization (point spacing). /// This only modifies the [Linspace] dimensions. /// To modify the sampling, use [Self::resample_mut]. pub fn stretch_mut(&mut self, factor: f64) -> Result<(), Error> { - if factor.is_sign_negative() { - return Err(Error::NegativeStretchFactor); + if !factor.is_normal() { + return Err(Error::InvalidStretchFactor); } self.start *= factor; self.end *= factor; @@ -81,10 +81,10 @@ impl Linspace { } /// Stretch this [Linspace] to a new [Linspace] definition. - /// The streetching factor must be a positive number. - /// this preserves the initial grid quantization (point spacing). - /// This only modifies the [Linspace] dimensions. - /// To modify the sampling, use [Self::resampled]. + /// The streetching factor must be a finite positive (possibly fractional) number. + /// This preserves the initial grid quantization (point spacing) and + /// only affects the [Linspace] dimensions. Use [Self::resampled] to resample + /// the [Linspace]. pub fn stretched(&self, factor: f64) -> Result { let mut s = self.clone(); s.stretch_mut(factor)?; @@ -92,11 +92,11 @@ impl Linspace { } /// Resample this mutable [Linspace] so the point spacing is modified. - /// This is a multiplicative stretching factor, which must be a positive number. + /// This is a multiplicative stretching factor, which must be a finite positive number (possibly fractional). /// This does not modify the [Linspace] dimensions, use [Self::stretch_mut] to modify dimensions. pub fn resample_mut(&mut self, factor: f64) -> Result<(), Error> { - if factor.is_sign_negative() { - return Err(Error::NegativeStretchFactor); + if !factor.is_normal() { + return Err(Error::InvalidStretchFactor); } self.spacing *= factor; Ok(()) @@ -104,7 +104,7 @@ impl Linspace { /// Resample this [Linspace] to a new [Linspace] definition, modifying the point /// spacing but preserving the initial dimensions. The resampling factor - /// must be positive, possibly fractional number. This does not modify the dimensions, + /// must be a finite positive, possibly fractional number. This does not modify the dimensions, /// use [Self::stretched] to modify the dimensions. pub fn resampled(&self, factor: f64) -> Result { let mut s = self.clone(); From 7fd7a3688fe49c6dfecc53520f449e654b2e3f77 Mon Sep 17 00:00:00 2001 From: "Guillaume W. Bres" Date: Sun, 2 Nov 2025 11:47:04 +0100 Subject: [PATCH 14/25] upscaling Signed-off-by: Guillaume W. Bres --- src/grid.rs | 19 +++++ src/lib.rs | 166 ++++++++++++++++++++++++++++++++-------- src/linspace.rs | 69 ++++++++++++----- src/tests/mod.rs | 1 + src/tests/roi.rs | 21 ++--- src/tests/stretching.rs | 130 +++++++++++++++++++++++++++++++ 6 files changed, 346 insertions(+), 60 deletions(-) create mode 100644 src/tests/stretching.rs diff --git a/src/grid.rs b/src/grid.rs index 1b72604..6335257 100644 --- a/src/grid.rs +++ b/src/grid.rs @@ -3,6 +3,25 @@ use crate::{error::ParsingError, linspace::Linspace}; #[cfg(feature = "serde")] use serde::{Deserialize, Serialize}; +#[derive(Debug, Copy, Clone, PartialEq, PartialOrd)] +#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))] +pub enum Axis { + /// Latitude axis (y) + Latitude, + + /// Longitude axis (x) + Longitude, + + /// Altitude axis (z), only meaningful in true 3D IONEX + Altitude, + + /// Latitude + Longitude axes plane + Planar, + + /// All (Latitude + Longitude plane, and z-axis when feasible) + All, +} + /// [Grid] used to describe latitude, longitude /// and altitude linar spaces, defining the entire map. #[derive(Debug, Copy, Clone, Default, PartialEq, PartialOrd)] diff --git a/src/lib.rs b/src/lib.rs index bded9d6..74c124a 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -62,8 +62,9 @@ use hifitime::prelude::Epoch; use crate::{ cell::{MapCell, MapPoint}, coordinates::QuantizedCoordinates, - error::{FormattingError, ParsingError}, + error::{Error, FormattingError, ParsingError}, file_attributes::{FileAttributes, Region}, + grid::{Axis, Grid}, header::Header, key::Key, quantized::Quantized, @@ -76,9 +77,9 @@ pub mod prelude { pub use crate::{ bias::BiasSource, cell::MapCell, - error::{FormattingError, ParsingError, Error}, + error::{Error, FormattingError, ParsingError}, file_attributes::*, - grid::Grid, + grid::{Axis, Grid}, header::Header, ionosphere::IonosphereParameters, key::Key, @@ -529,7 +530,7 @@ impl IONEX { false } - /// Returns map borders as a [Rect]angle, with coordinates in decimal degrees. + /// Returns map borders as a [Rect]angle, with coordinates in decimal degrees. /// This uses the [Header] description and assumes all maps are within these borders. pub fn bounding_rect_degrees(&self) -> Rect { Rect::new( @@ -554,7 +555,7 @@ impl IONEX { !self.is_worldwide_map() } - /// Stretch this [IONEX] definition so it becomes compatible + /// Stretch this [IONEX] definition so it becomes compatible /// with the description of a Global/Worldwide [IONEX]. pub fn to_worldwide_ionex(&self) -> IONEX { let mut ionex = self.clone(); @@ -595,6 +596,8 @@ impl IONEX { // copy header ionex.header = self.header.clone(); + // stretch factors + // rework header ionex.header.grid.latitude.start = max_lat; ionex.header.grid.latitude.end = min_lat; @@ -621,33 +624,136 @@ impl IONEX { Some(ionex) } - /// Stretch this [IONEX] into a new file [IONEX], modifying grid - /// spatial precision. - /// ## Inputs - /// - stretch > 1.0: increases precision. - /// The 2D planar interpolation is applied to define the TEC as originally - /// non existing coordinates. - /// - stretch < 1.0: reduces precision. - /// NB: work in progress, not validated yet. - pub fn grid_precision_stretch(&self, stretch_factor: f64) -> IONEX { + /// Modify the grid dimensions by a positive, possibly fractional number, + /// and interpolates the TEC values. + /// + /// This does not modify the grid quantization, only the dimensions. + /// + /// ## Input + /// - axis: [Axis] to be stretched + /// - factor: + /// - > 1.0: upscaling + /// - < 1.0: downscaling + /// - 0.0: invalid + /// + /// Example 1: spatial upscaling + /// ``` + /// use ionex::prelude::{ + /// IONEX, + /// Axis, + /// }; + /// + /// // grab a local dataset, upscaling hardly applies to worldwide maps.. + /// // this example is a data/IONEX/V1/CKMG0020.22I.gz that went through a x0.5 downscale. + /// let ionex = IONEX::from_gzip_file("data/IONEX/V1/CKMR0020.22I.gz").unwrap_or_else(|e| { + /// panic!("Failed to parse CKMR0020: {}", e); + /// }); + /// + /// // spatial dimensions upscale including interpolation: upscales to worldwide dimensions + /// let upscaled = ionex.spatially_stretched(Axis::Planar, 2.0) + /// .unwrap(); + /// + /// // testbench: compares these results to the original data + /// ``` + pub fn spatially_stretched(&self, axis: Axis, factor: f64) -> Result { let mut s = self.clone(); - s.grid_precision_stretch_mut(stretch_factor); - s + s.spatial_stretching_mut(axis, factor)?; + Ok(s) } - /// Stretch this mutable [IONEX], modifying the grid precision. - /// See [Self::grid_precision_stretch] for more information. - /// NB: work in progress, not validated yet. - pub fn grid_precision_stretch_mut(&mut self, stretch_factor: f64) { - // update grid - self.header.grid.latitude.start *= stretch_factor; - self.header.grid.latitude.end *= stretch_factor; - self.header.grid.longitude.start *= stretch_factor; - self.header.grid.longitude.end *= stretch_factor; + /// Modify the grid spacing (quantization) while preserving the dimensions, + /// and interpolates the TEC values. + /// + /// This only modify the grid quantization, the map dimensions are preserved. + /// + /// ## Input + /// - axis: [Axis] to be resampled + /// - factor: + /// - > 1.0: upscaling + /// - < 1.0: downscaling + /// - 0.0: invalid + pub fn spatially_resampled(&self, axis: Axis, factor: f64) -> Result { + let mut s = self.clone(); + s.spatial_resampling_mut(axis, factor)?; + Ok(s) + } - // // update map - // for epoch in self.record.epochs_iter() { - // } + /// Modify the grid dimensions by a positive, possibly fractional number, + /// and interpolates the TEC values. + /// + /// This does not modify the grid quantization, only the dimensions. + /// + /// ## Input + /// - axis: [Axis] to be stretched + /// - factor: + /// - > 1.0: upscaling + /// - < 1.0: downscaling + /// - 0.0: invalid + pub fn spatial_stretching_mut(&mut self, axis: Axis, factor: f64) -> Result<(), Error> { + // Stretch header specs + match axis { + Axis::Latitude => { + self.header.grid.latitude.stretch_mut(factor)?; + }, + Axis::Longitude => { + self.header.grid.longitude.stretch_mut(factor)?; + }, + Axis::Altitude => { + self.header.grid.altitude.stretch_mut(factor)?; + }, + Axis::Planar => { + self.header.grid.latitude.stretch_mut(factor)?; + self.header.grid.longitude.stretch_mut(factor)?; + }, + Axis::All => { + self.header.grid.latitude.stretch_mut(factor)?; + self.header.grid.longitude.stretch_mut(factor)?; + self.header.grid.altitude.stretch_mut(factor)?; + }, + } + + // TODO: data interpolation + + Ok(()) + } + + /// Modify the grid spacing (quantization) while preserving the dimensions, + /// and interpolates the TEC values. + /// + /// This only modify the grid quantization, the map dimensions are preserved. + /// + /// ## Input + /// - axis: [Axis] to be resampled + /// - factor: + /// - > 1.0: upscaling + /// - < 1.0: downscaling + /// - 0.0: invalid + pub fn spatial_resampling_mut(&mut self, axis: Axis, factor: f64) -> Result<(), Error> { + // Stretch header specs + match axis { + Axis::Latitude => { + self.header.grid.latitude.resample_mut(factor)?; + }, + Axis::Longitude => { + self.header.grid.longitude.resample_mut(factor)?; + }, + Axis::Altitude => { + self.header.grid.altitude.resample_mut(factor)?; + }, + Axis::Planar => { + self.header.grid.latitude.resample_mut(factor)?; + self.header.grid.longitude.resample_mut(factor)?; + }, + Axis::All => { + self.header.grid.latitude.resample_mut(factor)?; + self.header.grid.longitude.resample_mut(factor)?; + self.header.grid.altitude.resample_mut(factor)?; + }, + } + + // TODO: data interpolation + + Ok(()) } /// Designs a [MapCell] iterator (micro rectangle region made of 4 local points) @@ -960,8 +1066,7 @@ mod test { let roi = Rect::new(coord!(x: -180.0, y: -85.0), coord!(x: 180.0, y: -82.5)); - let regional = ionex.to_regional_ionex(roi.into()) - .unwrap(); + let regional = ionex.to_regional_ionex(roi.into()).unwrap(); // dump regional.to_file("region.txt").unwrap(); @@ -971,5 +1076,4 @@ mod test { panic!("Failed to parse region.txt: {}", e); }); } - } diff --git a/src/linspace.rs b/src/linspace.rs index 1b4f90c..ac81314 100644 --- a/src/linspace.rs +++ b/src/linspace.rs @@ -1,6 +1,9 @@ use std::ops::Rem; -use crate::{error::{ParsingError, Error}, quantized::Quantized}; +use crate::{ + error::{Error, ParsingError}, + quantized::Quantized, +}; /// Quantized Linspace for iteration #[derive(Debug, Copy, Clone, Default, PartialEq, PartialOrd)] @@ -80,9 +83,9 @@ impl Linspace { Ok(()) } - /// Stretch this [Linspace] to a new [Linspace] definition. + /// Stretch this [Linspace] to a new [Linspace] definition. /// The streetching factor must be a finite positive (possibly fractional) number. - /// This preserves the initial grid quantization (point spacing) and + /// This preserves the initial grid quantization (point spacing) and /// only affects the [Linspace] dimensions. Use [Self::resampled] to resample /// the [Linspace]. pub fn stretched(&self, factor: f64) -> Result { @@ -102,7 +105,7 @@ impl Linspace { Ok(()) } - /// Resample this [Linspace] to a new [Linspace] definition, modifying the point + /// Resample this [Linspace] to a new [Linspace] definition, modifying the point /// spacing but preserving the initial dimensions. The resampling factor /// must be a finite positive, possibly fractional number. This does not modify the dimensions, /// use [Self::stretched] to modify the dimensions. @@ -241,39 +244,67 @@ mod test { linspace.stretch_mut(0.5).unwrap(); assert_eq!(linspace.minmax(), (-90.0, 90.0)); - assert_eq!(linspace.spacing, 5.0, "linspace quantization not preserved!"); - + assert_eq!( + linspace.spacing, 5.0, + "linspace quantization not preserved!" + ); + linspace.stretch_mut(0.75).unwrap(); assert_eq!(linspace.minmax(), (-67.5, 67.5)); - assert_eq!(linspace.spacing, 5.0, "linspace quantization not preserved!"); - + assert_eq!( + linspace.spacing, 5.0, + "linspace quantization not preserved!" + ); + linspace.stretch_mut(0.5).unwrap(); assert_eq!(linspace.minmax(), (-33.75, 33.75)); - assert_eq!(linspace.spacing, 5.0, "linspace quantization not preserved!"); - + assert_eq!( + linspace.spacing, 5.0, + "linspace quantization not preserved!" + ); + linspace.stretch_mut(2.0).unwrap(); assert_eq!(linspace.minmax(), (-67.5, 67.5)); - assert_eq!(linspace.spacing, 5.0, "linspace quantization not preserved!"); + assert_eq!( + linspace.spacing, 5.0, + "linspace quantization not preserved!" + ); } - + #[test] fn linspace_resampling() { let mut linspace = Linspace::new(-180.0, 180.0, 5.0).unwrap(); linspace.resample_mut(0.5).unwrap(); assert_eq!(linspace.spacing, 2.5); - assert_eq!(linspace.minmax(), (-180.0, 180.0), "dimensions not preserved!"); - + assert_eq!( + linspace.minmax(), + (-180.0, 180.0), + "dimensions not preserved!" + ); + linspace.resample_mut(0.5).unwrap(); assert_eq!(linspace.spacing, 1.25); - assert_eq!(linspace.minmax(), (-180.0, 180.0), "dimensions not preserved!"); - + assert_eq!( + linspace.minmax(), + (-180.0, 180.0), + "dimensions not preserved!" + ); + linspace.resample_mut(2.0).unwrap(); assert_eq!(linspace.spacing, 2.5); - assert_eq!(linspace.minmax(), (-180.0, 180.0), "dimensions not preserved!"); - + assert_eq!( + linspace.minmax(), + (-180.0, 180.0), + "dimensions not preserved!" + ); + linspace.resample_mut(2.0).unwrap(); assert_eq!(linspace.spacing, 5.0); - assert_eq!(linspace.minmax(), (-180.0, 180.0), "dimensions not preserved!"); + assert_eq!( + linspace.minmax(), + (-180.0, 180.0), + "dimensions not preserved!" + ); } } diff --git a/src/tests/mod.rs b/src/tests/mod.rs index 0f1a3e6..b77be60 100644 --- a/src/tests/mod.rs +++ b/src/tests/mod.rs @@ -5,6 +5,7 @@ mod filename; mod parsing; mod qc; mod roi; +mod stretching; mod v1; diff --git a/src/tests/roi.rs b/src/tests/roi.rs index 3051c1c..a7a2b24 100644 --- a/src/tests/roi.rs +++ b/src/tests/roi.rs @@ -1,10 +1,6 @@ use crate::{ + prelude::{coord, Rect, IONEX}, tests::init_logger, - prelude::{ - IONEX, - Rect, - coord, - }, }; #[test] @@ -14,9 +10,12 @@ fn worldwide_map() { let mut ionex = IONEX::from_gzip_file("data/IONEX/V1/CKMG0020.22I.gz").unwrap_or_else(|e| { panic!("Failed to parse CKMG0020: {}", e); }); - + // returned from file name directly - assert!(ionex.is_worldwide_map(), "Worldwide map from standardized file name!"); + assert!( + ionex.is_worldwide_map(), + "Worldwide map from standardized file name!" + ); // destroy and verify this still works ionex.attributes = None; @@ -25,10 +24,12 @@ fn worldwide_map() { // reduce to ROI let roi = Rect::new(coord!(x: -30.0, y: -30.0), coord!(x: 30.0, y: 30.0)); - let reduced = ionex.to_regional_ionex(roi.into()) - .unwrap(); + let reduced = ionex.to_regional_ionex(roi.into()).unwrap(); - assert!(reduced.is_regional_map(), "Worldwide map reduced to Regional ROI"); + assert!( + reduced.is_regional_map(), + "Worldwide map reduced to Regional ROI" + ); // test new map let bounding_rect = reduced.bounding_rect_degrees(); diff --git a/src/tests/stretching.rs b/src/tests/stretching.rs new file mode 100644 index 0000000..879a690 --- /dev/null +++ b/src/tests/stretching.rs @@ -0,0 +1,130 @@ +use crate::{ + prelude::{coord, Axis, Duration, MappingFunction, Rect, Version, IONEX}, + tests::{ + init_logger, + toolkit::{generic_comparison, generic_test, TestPoint}, + }, +}; + +use std::fs::File; +use std::io::BufWriter; + +#[test] +fn unitary_spatial() { + init_logger(); + + let ionex = IONEX::from_gzip_file("data/IONEX/V1/CKMG0020.22I.gz").unwrap_or_else(|e| { + panic!("Failed to parse CKMG0020: {}", e); + }); + + // dummy case + let stretched = ionex + .spatially_stretched(Axis::All, 1.0) + .unwrap_or_else(|e| { + panic!("unexpected error: {}", e); + }); + + assert_eq!(stretched, ionex); +} + +#[test] +fn spatial_planar_downscaling() { + init_logger(); + + // 2D IONEX + let ionex = IONEX::from_gzip_file("data/IONEX/V1/CKMG0020.22I.gz").unwrap_or_else(|e| { + panic!("Failed to parse CKMG0020: {}", e); + }); + + for factor in [0.1, 0.2, 0.25, 0.33, 0.45, 0.5, 0.6, 0.7, 0.8, 0.9] { + let bounding_rect = ionex.bounding_rect_degrees(); + let dimensions_deg = (bounding_rect.width(), bounding_rect.height()); + + // Verify that z-axis does not modify anything (2D case) + let stretched = ionex + .spatially_stretched(Axis::Altitude, factor) + .unwrap_or_else(|e| { + panic!("unexpected error: {}", e); + }); + + assert_eq!(stretched.record, ionex.record, "z-axis on 2D IONEX"); + assert_eq!( + stretched.bounding_rect_degrees(), + bounding_rect, + "z-axis on 2D IONEX" + ); + + let stretched = ionex + .spatially_stretched(Axis::All, factor) + .unwrap_or_else(|e| { + panic!("unexpected error: {}", e); + }); + + let new_dimensions_deg = ( + stretched.bounding_rect_degrees().width(), + stretched.bounding_rect_degrees().height(), + ); + + assert_eq!( + new_dimensions_deg.0, + dimensions_deg.0 * factor, + "incorrect post-stretching dimensions" + ); + assert_eq!( + new_dimensions_deg.1, + dimensions_deg.1 * factor, + "incorrect post-stretching dimensions" + ); + } +} + +#[test] +fn spatial_planar_upscaling() { + init_logger(); + + // 2D IONEX + let ionex = IONEX::from_gzip_file("data/IONEX/V1/CKMG0020.22I.gz").unwrap_or_else(|e| { + panic!("Failed to parse CKMG0020: {}", e); + }); + + for factor in [1.1, 1.2, 1.25, 1.5, 1.75, 2.0, 3.0, 4.0] { + let bounding_rect = ionex.bounding_rect_degrees(); + let dimensions_deg = (bounding_rect.width(), bounding_rect.height()); + + // Verify that z-axis does not modify anything (2D case) + let stretched = ionex + .spatially_stretched(Axis::Altitude, factor) + .unwrap_or_else(|e| { + panic!("unexpected error: {}", e); + }); + + assert_eq!(stretched.record, ionex.record, "z-axis on 2D IONEX"); + assert_eq!( + stretched.bounding_rect_degrees(), + bounding_rect, + "z-axis on 2D IONEX" + ); + + let stretched = ionex + .spatially_stretched(Axis::All, factor) + .unwrap_or_else(|e| { + panic!("unexpected error: {}", e); + }); + + let new_dimensions_deg = ( + stretched.bounding_rect_degrees().width(), + stretched.bounding_rect_degrees().height(), + ); + + assert_eq!( + new_dimensions_deg.0, + dimensions_deg.0 * factor, + "incorrect post-stretching dimensions" + ); + assert_eq!( + new_dimensions_deg.1, + dimensions_deg.1 * factor, + "incorrect post-stretching dimensions" + ); + } +} From ae2a6983f9314887625180be4c0898e61d400339 Mon Sep 17 00:00:00 2001 From: "Guillaume W. Bres" Date: Sun, 2 Nov 2025 13:20:32 +0100 Subject: [PATCH 15/25] working on api Signed-off-by: Guillaume W. Bres --- src/cell.rs | 114 ++++++++++++++++++++++++++++++++++++++++- src/error.rs | 6 +++ src/file_attributes.rs | 35 +++++++++---- src/grid.rs | 22 ++++++-- src/lib.rs | 101 +++++++++++++++++++++++++++++++++--- 5 files changed, 252 insertions(+), 26 deletions(-) diff --git a/src/cell.rs b/src/cell.rs index 3c96a4a..8f4c8b9 100644 --- a/src/cell.rs +++ b/src/cell.rs @@ -1,6 +1,6 @@ use geo::{Contains, GeodesicArea, Geometry, Point, Rect}; -use crate::prelude::{Epoch, TEC}; +use crate::prelude::{Epoch, Error, TEC}; #[derive(Debug, Default, Clone, Copy, PartialEq)] pub struct MapPoint { @@ -11,7 +11,9 @@ pub struct MapPoint { pub tec: TEC, } -/// [MapCell] describing a region that we can then interpolate. +/// [MapCell] describing a 4 corner region that we can interpolate. +/// In the processing workflow, [MapCell]s are constructed from individual +/// quanta (smallest ROI) described in a IONEX map. #[derive(Debug, Default, Clone, Copy, PartialEq)] pub struct MapCell { /// Epoch of observation @@ -116,6 +118,26 @@ impl MapCell { } } + /// Returns true if both [MapCell]s describe the same spatial ROI + pub fn spatial_match(&self, rhs: &Self) -> bool { + if self.north_east.point == rhs.north_east.point { + if self.north_west.point == rhs.north_west.point { + if self.south_east.point == rhs.south_east.point { + if self.south_west.point == rhs.south_west.point { + return true; + } + } + } + } + + false + } + + /// Returns true if both [MapCell]s describe the same point in time + pub fn temporal_match(&self, rhs: &Self) -> bool { + self.epoch == rhs.epoch + } + /// Defines a unitary [MapCell] ((0,0), (0,1), (1,0), (1,1)) with associated TEC values, /// where /// - (x/long=0, y/lat=0) is the SW corner @@ -327,6 +349,94 @@ impl MapCell { TEC::from_tecu(tecu) } + // /// Interpolates two [MapCell]s that must describe the same area, + // /// but a different point in time. + // /// + // /// ## Input + // /// - epoch: [Epoch] of interpolation, must be temporally in between + // /// [Self] and rhs. + // pub fn temporally_interpolated(&self, epoch: Epoch, rhs: &Self) -> Result { + // if !self.spatial_match(rhs) { + // return Err(Error::SpatialMismatch); + // } + + // let (min_t, max_t) = + // ( + // std::cmp::min(self.epoch, rhs.epoch), + // std::cmp::max(self.epoch, rhs.epoch), + // ); + + // if epoch < min_t && epoch > max_t { + // return Err(Error::InvalidTemporalPoint); + // } + + // let (num_1, num_2, dt) = if self.epoch > rhs.epoch { + // ( + // (self.epoch - epoch).to_seconds(), + // (epoch - rhs.epoch).to_seconds(), + // (self.epoch - rhs.epoch).to_seconds(), + // ) + // } else { + // ( + // (rhs.epoch - epoch).to_seconds(), + // (epoch - self.epoch).to_seconds(), + // (rhs.epoch - self.epoch).to_seconds(), + // ) + // }; + + // let (ne_1, ne_2) = if self.epoch > rhs.epoch { + // (self.north_east.tec.tecu(), rhs.north_east.tec.tecu()) + // } else { + // (rhs.north_east.tec.tecu(), self.north_east.tec.tecu()) + // }; + // + // let (nw_1, nw_2) = if self.epoch > rhs.epoch { + // (self.north_west.tec.tecu(), rhs.north_west.tec.tecu()) + // } else { + // (rhs.north_west.tec.tecu(), self.north_west.tec.tecu()) + // }; + + // let (se_1, se_2) = if self.epoch > rhs.epoch { + // (self.south_east.tec.tecu(), rhs.south_east.tec.tecu()) + // } else { + // (rhs.south_east.tec.tecu(), self.south_east.tec.tecu()) + // }; + // + // let (sw_1, sw_2) = if self.epoch > rhs.epoch { + // (self.south_west.tec.tecu(), rhs.south_west.tec.tecu()) + // } else { + // (rhs.south_west.tec.tecu(), self.south_west.tec.tecu()) + // }; + + // Ok(Self { + // epoch, + // north_east: MapPoint { + // point: self.north_east.point, + // tec: TEC::from_tecu( + // num_1 * ne_1 /dt + num_2 * ne_2 /dt + // ), + // }, + // north_west: MapPoint { + // point: self.north_west.point, + // tec: TEC::from_tecu( + // num_1 * nw_1 /dt + num_2 * nw_2 /dt + // ), + // }, + // south_east: MapPoint { + // point: self.south_east.point, + // tec: TEC::from_tecu( + // num_1 * se_1 /dt + num_2 * se_2 /dt + // ), + // }, + // south_west: MapPoint { + // point: self.south_west.point, + // tec: TEC::from_tecu( + // num_1 * sw_1 /dt + num_2 * sw_2 /dt + // ), + // }, + // }) + // } + /// Spatial + Temporal Interpolation of [TEC] value using planery equation /// and rhs [MapCell], which should be closely sampled in time. /// [MapCell::contains] should be true both [MapCell]s and proposed geometry diff --git a/src/error.rs b/src/error.rs index c9f97d3..311beb7 100644 --- a/src/error.rs +++ b/src/error.rs @@ -95,6 +95,12 @@ pub enum ParsingError { pub enum Error { #[error("strech factor must be positive finite number")] InvalidStretchFactor, + + #[error("both regions do not describe the same spatial ROI")] + SpatialMismatch, + + #[error("invalid temporal interpolation instant")] + InvalidTemporalPoint, } /// Errors that may rise during Formatting process diff --git a/src/file_attributes.rs b/src/file_attributes.rs index a73d901..4082a3a 100644 --- a/src/file_attributes.rs +++ b/src/file_attributes.rs @@ -14,14 +14,14 @@ pub enum Region { /// Worldwide IONEX map #[default] - WorldWide, + Worldwide, } impl std::fmt::Display for Region { fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result { match self { Self::Regional => write!(f, "{}", 'R'), - Self::WorldWide => write!(f, "{}", 'G'), + Self::Worldwide => write!(f, "{}", 'G'), } } } @@ -29,7 +29,7 @@ impl std::fmt::Display for Region { /// File production attributes. Used when generating /// RINEX data that follows standard naming conventions, /// or attached to data parsed from such files. -#[derive(Debug, Default, Clone, PartialEq)] +#[derive(Debug, Clone, PartialEq)] pub struct FileAttributes { /// File agency pub agency: String, @@ -50,6 +50,19 @@ pub struct FileAttributes { pub gzip_compressed: bool, } +impl Default for FileAttributes { + fn default() -> Self { + Self { + doy: 1, + year: 2000, + agency: "XXX".to_string(), // valid + region: Default::default(), + #[cfg(feature = "flate2")] + gzip_compressed: Default::default(), + } + } +} + impl std::fmt::Display for FileAttributes { #[cfg(feature = "flate2")] fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result { @@ -106,7 +119,7 @@ impl std::str::FromStr for FileAttributes { .map_err(|_| Error::NonStandardFilename)?; let region = if filename[3..4].eq("G") { - Region::WorldWide + Region::Worldwide } else { Region::Regional }; @@ -140,8 +153,8 @@ impl gnss_qc_traits::Merge for FileAttributes { self.doy = rhs.doy; } - if self.region == Region::Regional && rhs.region == Region::WorldWide { - self.region = Region::WorldWide; + if self.region == Region::Regional && rhs.region == Region::Worldwide { + self.region = Region::Worldwide; } #[cfg(feature = "flate2")] @@ -161,9 +174,9 @@ mod test { #[test] fn standard_filenames() { for (filename, agency, year, doy, region) in [ - ("CKMG0020.22I", "CKM", 2022, 2, Region::WorldWide), - ("CKMG0090.21I", "CKM", 2021, 9, Region::WorldWide), - ("JPLG0010.17I", "JPL", 2017, 1, Region::WorldWide), + ("CKMG0020.22I", "CKM", 2022, 2, Region::Worldwide), + ("CKMG0090.21I", "CKM", 2021, 9, Region::Worldwide), + ("JPLG0010.17I", "JPL", 2017, 1, Region::Worldwide), ("JPLR0010.17I", "JPL", 2017, 1, Region::Regional), ("JPLR0010.17I", "JPL", 2017, 1, Region::Regional), ] { @@ -185,9 +198,9 @@ mod test { #[test] fn gzip_filenames() { for (filename, agency, year, doy, region) in [ - ("CKMG0020.22I.gz", "CKM", 2022, 2, Region::WorldWide), + ("CKMG0020.22I.gz", "CKM", 2022, 2, Region::Worldwide), ("CKMR0020.22I.gz", "CKM", 2022, 2, Region::Regional), - ("JPLG0010.17I.gz", "JPL", 2017, 1, Region::WorldWide), + ("JPLG0010.17I.gz", "JPL", 2017, 1, Region::Worldwide), ("CKMR0020.22I.gz", "CKM", 2022, 2, Region::Regional), ] { let attrs = FileAttributes::from_str(filename).unwrap(); diff --git a/src/grid.rs b/src/grid.rs index 6335257..f24061c 100644 --- a/src/grid.rs +++ b/src/grid.rs @@ -38,15 +38,27 @@ pub struct Grid { } impl Grid { - /// Returns true if self is compatible with a 3D TEC map. + /// Returns true when this [Grid] is not compatible with a 3D TEC map. + /// That means the altitude is a single point with null width. + pub fn is_2d_grid(&self) -> bool { + self.altitude.is_single_point() + } + + /// Returns true if this [Grid] is compatible with a 3D TEC map. pub fn is_3d_grid(&self) -> bool { !self.is_2d_grid() } - /// Returns true if self is not compatible with a 3D TEC map. - /// That means the altitude is a single point with null width. - pub fn is_2d_grid(&self) -> bool { - self.altitude.is_single_point() + /// Returns true if this [Grid] matches the description of a worldwide map. + pub fn is_worldwide(&self) -> bool { + let latitude_deg = self.latitude.minmax(); + let longitude_deg = self.longitude.minmax(); + latitude_deg == (-87.5, 87.5) && longitude_deg == (-180.0, 180.0) + } + + /// Returns true if this [Grid] does not match the description of a worldwide map. + pub fn is_regional(&self) -> bool { + !self.is_worldwide() } /// Defines a new [Grid] with updated latitude space diff --git a/src/lib.rs b/src/lib.rs index 74c124a..7d30483 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -312,7 +312,11 @@ impl IONEX { let year = first_epoch.year(); let doy = first_epoch.day_of_year().round() as u32; - let region = Region::WorldWide; // TODO: study the grid specs + let region = if self.header.grid.is_worldwide() { + Region::Worldwide + } else { + Region::Regional + }; Some(FileAttributes { doy, @@ -505,9 +509,9 @@ impl IONEX { /// * [Self::guess_file_attributes] in close relation: helps generate standardized /// filenames for files that did now follow the convention /// * [gnss_qc_traits::Merge]: to combine two files to one another - /// * [Self::to_regional_roi]: to reduce a larger (possibly [Region::GlobalWorldWide] map) + /// * [Self::to_regional_roi]: to reduce a larger (possibly [Region::GlobalWorldwide] map) /// to a given ROI - /// * [Self::to_worldwide_map]: to enlarge a possibly [Region::Regional] map to a [Region::GlobalWorldWide] IONEX. + /// * [Self::to_worldwide_map]: to enlarge a possibly [Region::Regional] map to a [Region::GlobalWorldwide] IONEX. #[cfg(feature = "flate2")] #[cfg_attr(docsrs, doc(cfg(feature = "flate2")))] pub fn to_gzip_file>(&self, path: P) -> Result<(), FormattingError> { @@ -543,7 +547,7 @@ impl IONEX { /// (as opposed to a local/regional ROI). pub fn is_worldwide_map(&self) -> bool { if let Some(attributes) = &self.attributes { - attributes.region == Region::WorldWide + attributes.region == Region::Worldwide } else { let bounding_rect = self.bounding_rect_degrees(); bounding_rect.width() == 360.0 && bounding_rect.height() == 175.0 @@ -561,7 +565,7 @@ impl IONEX { let mut ionex = self.clone(); if let Some(attributes) = &mut ionex.attributes { - attributes.region = Region::WorldWide; + attributes.region = Region::Worldwide; } // update grid specs, preserve accuracy @@ -690,7 +694,7 @@ impl IONEX { /// - < 1.0: downscaling /// - 0.0: invalid pub fn spatial_stretching_mut(&mut self, axis: Axis, factor: f64) -> Result<(), Error> { - // Stretch header specs + // Update the header specs match axis { Axis::Latitude => { self.header.grid.latitude.stretch_mut(factor)?; @@ -712,11 +716,38 @@ impl IONEX { }, } - // TODO: data interpolation + // update the attributes + match self.attributes { + Some(ref mut attributes) => { + if self.header.grid.is_worldwide() { + attributes.region = Region::Worldwide; + } else { + attributes.region = Region::Regional; + } + }, + None => { + let mut attributes = FileAttributes::default(); + if self.header.grid.is_worldwide() { + attributes.region = Region::Worldwide; + } else { + attributes.region = Region::Regional; + } + + self.attributes = Some(attributes); + }, + } + + // synchronous spatial interpolation + for epoch in self.epoch_iter() {} Ok(()) } + /// Returns an iterator over [Epoch]s in chronological order. + pub fn epoch_iter(&self) -> Box + '_> { + Box::new(self.record.map.keys().map(|k| k.epoch).unique().sorted()) + } + /// Modify the grid spacing (quantization) while preserving the dimensions, /// and interpolates the TEC values. /// @@ -751,11 +782,65 @@ impl IONEX { }, } - // TODO: data interpolation + // synchronous spatial interpolation + for epoch in self.epoch_iter() {} Ok(()) } + /// Upscale (upsample) or Downscale (downsample) this mutable [IONEX], + /// modifying the stretch on the temporal axis. + /// + /// ## Input + /// - factor: a positive finite number + /// - 0.0: is not valid + /// - >1.0: upscaling case. For example, 2.0 means x2 sample rate increase (+100%). + /// - and 1.5 means +50% sample rate increase. + /// - <1.0: downscaling case. For example, 0.5 means /2 sample rate decrease (-50%). + pub fn temporal_stretching_mut(&mut self, factor: f64) -> Result<(), Error> { + if !factor.is_normal() { + return Err(Error::InvalidStretchFactor); + } + + // update header + self.header.sampling_period = self.header.sampling_period / factor; + + // TODO: (t1, t2) iter + Ok(()) + } + + /// Stretchs this mutable [IONEX] both in temporal and spatial dimensions, + /// modifying both dimensions. + /// When factor > 1.0 this is an upscaling operation, when factor < 1.0, + /// this is a downscaling operation. + /// This uses both spatial and temporal interpolation to form valid data points. + pub fn temporal_spatial_stretching_mut( + &mut self, + temporal_factor: f64, + axis: Axis, + spatial_factor: f64, + ) -> Result<(), Error> { + self.spatial_stretching_mut(axis, spatial_factor)?; + self.temporal_stretching_mut(temporal_factor)?; + Ok(()) + } + + /// Stretchs this [IONEX] into a new [IONEX], both temporally and spatially stretched, + /// modifying both dimensions. + /// When factor > 1.0 this is an upscaling operation, when factor < 1.0, + /// this is a downscaling operation. + /// This uses both spatial and temporal interpolation to form valid data points. + pub fn temporally_spatially_stretched( + &self, + temporal_factor: f64, + axis: Axis, + spatial_factor: f64, + ) -> Result { + let mut s = self.clone(); + s.temporal_spatial_stretching_mut(temporal_factor, axis, spatial_factor)?; + Ok(s) + } + /// Designs a [MapCell] iterator (micro rectangle region made of 4 local points) /// that can then be interpolated. pub fn map_cell_iter(&self) -> Box + '_> { From 836aba56d384f445eb7d9d59306d4aaa39d26b06 Mon Sep 17 00:00:00 2001 From: "Guillaume W. Bres" Date: Sun, 2 Nov 2025 18:02:48 +0100 Subject: [PATCH 16/25] upgrading API Signed-off-by: Guillaume W. Bres --- src/cell.rs | 122 ++++++++++++++--- src/error.rs | 9 ++ src/file_attributes.rs | 14 ++ src/lib.rs | 297 ++++++++++++++++++++++++----------------- 4 files changed, 300 insertions(+), 142 deletions(-) diff --git a/src/cell.rs b/src/cell.rs index 8f4c8b9..883dc52 100644 --- a/src/cell.rs +++ b/src/cell.rs @@ -2,6 +2,22 @@ use geo::{Contains, GeodesicArea, Geometry, Point, Rect}; use crate::prelude::{Epoch, Error, TEC}; +#[derive(Debug, Copy, Clone, PartialEq, PartialOrd)] +#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))] +pub enum Cardinal { + /// NE [Cardinal] + NorthEast, + + /// NW [Cardinal] + NorthWest, + + /// SE [Cardinal] + SouthEast, + + /// SW [Cardinal] + SouthWest, +} + #[derive(Debug, Default, Clone, Copy, PartialEq)] pub struct MapPoint { /// [Point] @@ -138,6 +154,44 @@ impl MapCell { self.epoch == rhs.epoch } + /// Returns true if self is the northern neighbor of provided rhs [MapCell]. + pub fn is_northern_neighbor(&self, rhs: &Self) -> bool { + rhs.north_east.point == self.south_east.point + && rhs.north_west.point == self.south_west.point + } + + /// Returns true if self is the southern neighbor of provided rhs [MapCell]. + pub fn is_southern_neighbor(&self, rhs: &Self) -> bool { + rhs.south_east.point == self.north_east.point + && rhs.south_west.point == self.north_west.point + } + + /// Returns true if self is the easthern neighbor of provided rhs [MapCell]. + pub fn is_eastern_neighbor(&self, rhs: &Self) -> bool { + rhs.north_west.point == self.north_east.point + && rhs.south_west.point == self.south_east.point + } + + /// Returns true if self is the westhern neighbor of provided rhs [MapCell]. + pub fn is_western_neighbor(&self, rhs: &Self) -> bool { + rhs.north_east.point == self.north_west.point + && rhs.south_east.point == self.south_west.point + } + + /// Returns true if both cells are neighbors, meaning, they share two corners. + pub fn is_neighbor(&self, rhs: &Self) -> bool { + self.is_northern_neighbor(rhs) + || self.is_southern_neighbor(rhs) + || self.is_western_neighbor(rhs) + || self.is_eastern_neighbor(rhs) + } + + /// Returns true if this [MapCell] totally contains (wrapps) rhs cell. + /// Meaning, rhs is fully contained within Self. + pub fn contains_entirely(&self, rhs: &Self) -> bool { + self.contains(&Geometry::Rect(rhs.bounding_rect_degrees())) + } + /// Defines a unitary [MapCell] ((0,0), (0,1), (1,0), (1,1)) with associated TEC values, /// where /// - (x/long=0, y/lat=0) is the SW corner @@ -213,27 +267,28 @@ impl MapCell { /// Returns central [Point] of this [MapCell]. pub fn center(&self) -> Point { - geo::Point(self.borders().center()) + geo::Point(self.bounding_rect_degrees().center()) } - /// Returns borders of this [MapCell] expressed as a [Rect]angle. - pub fn borders(&self) -> Rect { + /// Returns borders of this [MapCell] expressed as a [Rect]angle, in decimal degrees. + pub fn bounding_rect_degrees(&self) -> Rect { Rect::new(self.south_west.point, self.north_east.point) } /// Returns geodesic perimeter (in meters) of this [MapCell]. pub fn geodesic_perimeter(&self) -> f64 { - self.borders().geodesic_perimeter() + self.bounding_rect_degrees().geodesic_perimeter() } /// Returns geodesic area (in squared meters) of this [MapCell]. pub fn geodesic_area(&self) -> f64 { - self.borders().geodesic_area_unsigned() + self.bounding_rect_degrees().geodesic_area_unsigned() } - /// Returns true if following [Geometry] is contained within this [MapCell]. + /// Returns true if following [Geometry], expressed in decimal degrees, + /// is contained within this [MapCell]. pub fn contains(&self, geometry: &Geometry) -> bool { - self.borders().contains(geometry) + self.bounding_rect_degrees().contains(geometry) } /// Copies and updates the Northeastern TEC component @@ -268,13 +323,13 @@ impl MapCell { /// Returns latitude span of this [MapCell] in degrees pub fn latitude_span_degrees(&self) -> f64 { - let borders = self.borders(); + let borders = self.bounding_rect_degrees(); borders.max().y - borders.min().y } /// Returns longitude span of this [MapCell] in degrees pub fn longitude_span_degrees(&self) -> f64 { - let borders = self.borders(); + let borders = self.bounding_rect_degrees(); borders.max().x - borders.min().x } @@ -331,7 +386,11 @@ impl MapCell { /// let tec = cell.spatial_interpolation(Point::new(0.01, 0.01)); /// assert_eq!(tec.tecu(), 0.9801); /// ``` - pub fn spatial_interpolation(&self, point: Point) -> TEC { + pub fn spatial_tec_interp(&self, point: Point) -> Result { + if !self.contains(&Geometry::Point(point)) { + return Err(Error::OutsideBoundaries); + } + let (latitude_span, longitude_span) = self.latitude_longitude_span_degrees(); let (p, q) = (point.y() / latitude_span, point.x() / longitude_span); @@ -346,7 +405,32 @@ impl MapCell { let tecu = (1.0 - p) * (1.0 - q) * e00 + p * (1.0 - q) * e10 + q * (1.0 - p) * e01 + p * q * e11; - TEC::from_tecu(tecu) + Ok(TEC::from_tecu(tecu)) + } + + /// Merges two neighboring [MapCell]s forming a new [MapCell] with best precision. + /// Both cells must be synchronous. + pub fn merge_neighbors(&self, rhs: &Self) -> Result { + if !self.temporal_match(rhs) { + return Err(Error::TemporalMismatch); + } + + // 1: determine the matching corners. The center + // point of the matching line is to become the new center. + // 2: form a new cell, of the boundary rect, interpolate + // both TEC at the new center + + if self.is_northern_neighbor(rhs) { + Ok(Self::default()) // TODO + } else if self.is_southern_neighbor(rhs) { + Ok(Self::default()) // TODO + } else if self.is_western_neighbor(rhs) { + Ok(Self::default()) // TODO + } else if self.is_eastern_neighbor(rhs) { + Ok(Self::default()) // TODO + } else { + Err(Error::SpatialMismatch) + } } // /// Interpolates two [MapCell]s that must describe the same area, @@ -476,16 +560,16 @@ impl MapCell { /// let tec = cell0.temporal_spatial_interpolation(t_ok, center, &cell1).unwrap(); /// assert_eq!(tec.tecu(), 1.0); /// ``` - pub fn temporal_spatial_interpolation( + pub fn temporal_spatial_tec_interp( &self, epoch: Epoch, - point: Point, + coordinates: Point, rhs: &Self, - ) -> Option { + ) -> Result { // interpolate at exact coordinates let (tecu_0, tecu_1) = ( - self.spatial_interpolation(point).tecu(), - rhs.spatial_interpolation(point).tecu(), + self.spatial_tec_interp(coordinates)?.tecu(), + rhs.spatial_tec_interp(coordinates)?.tecu(), ); if epoch >= self.epoch && epoch < rhs.epoch { @@ -495,7 +579,7 @@ impl MapCell { let tecu = (rhs.epoch - epoch).to_seconds() / dt * tecu_0 + (epoch - self.epoch).to_seconds() / dt * tecu_1; - Some(TEC::from_tecu(tecu)) + Ok(TEC::from_tecu(tecu)) } else if epoch >= rhs.epoch && epoch < self.epoch { // backwards let dt = (self.epoch - rhs.epoch).to_seconds(); @@ -503,9 +587,9 @@ impl MapCell { let tecu = (self.epoch - epoch).to_seconds() / dt * tecu_1 + (epoch - rhs.epoch).to_seconds() / dt * tecu_0; - Some(TEC::from_tecu(tecu)) + Ok(TEC::from_tecu(tecu)) } else { - None + Err(Error::TemporalMismatch) } } } diff --git a/src/error.rs b/src/error.rs index 311beb7..871b6c5 100644 --- a/src/error.rs +++ b/src/error.rs @@ -96,9 +96,18 @@ pub enum Error { #[error("strech factor must be positive finite number")] InvalidStretchFactor, + #[error("undefined ROI exterior boundaries")] + UndefinedBoundaries, + + #[error("provided coordinates are outside boundaries")] + OutsideBoundaries, + #[error("both regions do not describe the same spatial ROI")] SpatialMismatch, + #[error("both regions are not synchronous in time")] + TemporalMismatch, + #[error("invalid temporal interpolation instant")] InvalidTemporalPoint, } diff --git a/src/file_attributes.rs b/src/file_attributes.rs index 4082a3a..cca551a 100644 --- a/src/file_attributes.rs +++ b/src/file_attributes.rs @@ -50,6 +50,20 @@ pub struct FileAttributes { pub gzip_compressed: bool, } +impl FileAttributes { + /// Copies and returns a [FileAttributes] converted to global/worldwide map + pub fn globalized(mut self) -> Self { + self.region = Region::Worldwide; + self + } + + /// Copies and returns a [FileAttributes] converted to regional map + pub fn regionalized(mut self) -> Self { + self.region = Region::Regional; + self + } +} + impl Default for FileAttributes { fn default() -> Self { Self { diff --git a/src/lib.rs b/src/lib.rs index 7d30483..37a8ba7 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -60,7 +60,7 @@ use flate2::{read::GzDecoder, write::GzEncoder, Compression as GzCompression}; use hifitime::prelude::Epoch; use crate::{ - cell::{MapCell, MapPoint}, + cell::{Cardinal, MapCell, MapPoint}, coordinates::QuantizedCoordinates, error::{Error, FormattingError, ParsingError}, file_attributes::{FileAttributes, Region}, @@ -76,7 +76,7 @@ pub mod prelude { // export pub use crate::{ bias::BiasSource, - cell::MapCell, + cell::{Cardinal, MapCell}, error::{Error, FormattingError, ParsingError}, file_attributes::*, grid::{Axis, Grid}, @@ -581,39 +581,45 @@ impl IONEX { /// Reduce this [IONEX] definition so it is reduced to a regional ROI, /// described by a complex [Polygon] in decimal degrees. + /// The quantization (both spatial and temporal) is preserved, only the + /// spatial dimensions are modified by this operation. /// [Polygon::bounding_rect] must be defined for this operation to work correctly. - pub fn to_regional_ionex(&self, roi: Polygon) -> Option { + pub fn to_regional_ionex(&self, roi: Polygon) -> Result { let mut ionex = IONEX::default(); - let bounding_rect = roi.bounding_rect()?; + let bounding_rect = roi.bounding_rect().ok_or(Error::UndefinedBoundaries)?; let (min_long, min_lat) = (bounding_rect.min().x, bounding_rect.min().y); let (max_long, max_lat) = (bounding_rect.max().x, bounding_rect.max().y); - // copy attributes - ionex.attributes = self.attributes.clone(); - - if let Some(attributes) = &mut ionex.attributes { - attributes.region = Region::Regional; + // attributes + match &self.attributes { + Some(attributes) => { + ionex.attributes = Some(attributes.clone().regionalized()); + }, + None => { + let mut attributes = FileAttributes::default().regionalized(); + ionex.attributes = Some(attributes); + }, } // copy header ionex.header = self.header.clone(); - // stretch factors + // upscale so the grid boundaries become a perfect multiple + let long_upscaling = 0.0; + let lat_upscaling = 0.0; - // rework header - ionex.header.grid.latitude.start = max_lat; - ionex.header.grid.latitude.end = min_lat; - ionex.header.grid.longitude.start = min_long; - ionex.header.grid.longitude.end = max_long; + // spatial upscaling + ionex.spatial_stretching_mut(Axis::Latitude, lat_upscaling)?; + ionex.spatial_stretching_mut(Axis::Longitude, long_upscaling)?; - let roi = Geometry::Polygon(roi); + // simply restrict to ROI area + let geometry = Geometry::Polygon(roi); - // restrict area let cells = self .map_cell_iter() - .filter(|cell| cell.contains(&roi)) + .filter(|cell| cell.contains(&geometry)) .collect::>(); ionex.record = Record::from_map_cells( @@ -625,7 +631,7 @@ impl IONEX { &cells, ); - Some(ionex) + Ok(ionex) } /// Modify the grid dimensions by a positive, possibly fractional number, @@ -924,7 +930,52 @@ impl IONEX { ) } - /// Returns the [MapCell] that contains following [Geometry]. + /// Obtain a [MapCell] at specified point in time, wrapping this ROI. + /// The temporal instant must be within the current temporal axis. + /// The returned [MapCell] is a bounding rectangle, of 4 corners + /// correctly interpolated in space and time, than you can further manipulate. + /// Depending on the ROI dimensions, this is either: + /// - a unitary [MapCell]: the smallest region we can represent in this [IONEX]. + /// Prefer [Self::unitary_roi_at] in this case, to not bother designing the unit [Rect]. + /// - or the closest fitting [MapCell] (wrapping boundaries) that we could fit + pub fn roi_at(&self, epoch: Epoch, roi: Rect) -> Result { + if epoch < self.header.epoch_of_first_map { + return Err(Error::TemporalMismatch); + } + + if epoch > self.header.epoch_of_last_map { + return Err(Error::TemporalMismatch); + } + + // for each corner of the ROI, determine the best suited + // map cell for spatial interpolation + + // TODO + Err(Error::TemporalMismatch) + } + + /// Obtain the smallest [MapCell] (map quantization) at provided point in time + /// and coordinates. + /// + /// ## Input + /// - epoch: [Epoch] that must fit within the temporal axis. + /// If this instant aligns with the sampling rate, this process will not require temporal interpolation. + /// - coordinates: 2D [Point] that must lie within the map. + /// If the coordinates align with the grid, this process will not require spatial interpolation. + /// - card: [Cardinal] defining which corner these coordinates will represent in the returned [MapCell]. + /// For example: + /// - [Cardinal::NorthEast] means the northeastern component of fitted ROI [MapCell] will lie at these coordinates. + pub fn unitary_roi_at( + &self, + epoch: Epoch, + coordinates: Point, + card: Cardinal, + ) -> Result { + // TODO + Err(Error::TemporalMismatch) + } + + /// Obtain the best suited [MapCell] spatially wrapping this Geometry that contains following [Geometry]. /// /// ## Input /// - epoch: [Epoch] that must exist in this [IONEX] @@ -950,109 +1001,109 @@ impl IONEX { ) } - /// Interpolate TEC values for all discrete coordinates described by the following [LineString] - /// (in decimal degrees), at specific point in time that must exist within this record. - /// Otherwise, you should use [Self::temporal_spatial_area_interpolation] to also - /// use temporal interpolation from two existing data points. - /// ``` - /// use std::str::FromStr; - /// use ionex::prelude::{IONEX, Epoch, LineString, coord, Contains}; - /// - /// let ionex = IONEX::from_gzip_file("data/IONEX/V1/CKMG0020.22I.gz") - /// .unwrap(); - /// - /// // ROI (ddeg) must be within borders - /// let roi_ddeg = LineString::new(vec![ - /// coord! { x: -50.0, y: -23.0 }, - /// coord! { x: -50.25, y: -23.1 }, - /// coord! { x: -50.5, y: -23.2 }, - /// ]); - /// - /// // you can double check that - /// assert!(ionex.bounding_rect_degrees().contains(&roi_ddeg)); - /// - /// // Epoch must exist in the record - /// let noon = Epoch::from_str("2022-01-02T12:00:00 UTC") - /// .unwrap(); - /// - /// - /// ``` - pub fn spatial_area_interpolation( - &self, - area: &LineString, - epoch: Epoch, - ) -> BTreeMap { - let mut values = BTreeMap::new(); - - let fixed_altitude_km = self.header.grid.altitude.start; - - // for all requested coordinates - for point in area.points() { - let geo = Geometry::Point(point); - - let key = Key::from_decimal_degrees_km(epoch, point.y(), point.x(), fixed_altitude_km); - - for cell in self.synchronous_map_cell_iter(epoch) { - if cell.contains(&geo) { - let tec = cell.spatial_interpolation(point); - values.insert(key, tec); - } - } - } - - values - } - - /// Interpolate TEC values for all discrete coordinates described by the following [LineString] - /// (in decimal degrees), at specific point in time that does exist within this record. - /// We will interpolate the two neighbouring data points, which is not feasible at day boundaries. - pub fn temporal_spatial_area_interpolation( - &self, - area: &LineString, - epoch: Epoch, - ) -> BTreeMap { - let mut values = BTreeMap::new(); - - let (min_t, max_t) = ( - (epoch - self.header.sampling_period), - (epoch + self.header.sampling_period), - ); - - let (min_t, max_t) = ( - min_t.round(self.header.sampling_period), - max_t.round(self.header.sampling_period), - ); - - // for all requested coordinates - for point in area.points() { - let geo = Geometry::Point(point); - - for cell_t0 in self.synchronous_map_cell_iter(min_t) { - if cell_t0.contains(&geo) { - for cell_t1 in self.synchronous_map_cell_iter(max_t) { - if cell_t1.contains(&geo) { - if let Some(interpolated) = - cell_t0.temporal_spatial_interpolation(epoch, point, &cell_t1) - { - let key = Key { - epoch, - coordinates: QuantizedCoordinates::from_decimal_degrees( - point.y(), - point.x(), - self.header.grid.altitude.start, - ), - }; - - values.insert(key, interpolated); - } - } - } - } - } - } - - values - } + // /// Interpolate TEC values for all discrete coordinates described by the following [LineString] + // /// (in decimal degrees), at specific point in time that must exist within this record. + // /// Otherwise, you should use [Self::temporal_spatial_area_interpolation] to also + // /// use temporal interpolation from two existing data points. + // /// ``` + // /// use std::str::FromStr; + // /// use ionex::prelude::{IONEX, Epoch, LineString, coord, Contains}; + // /// + // /// let ionex = IONEX::from_gzip_file("data/IONEX/V1/CKMG0020.22I.gz") + // /// .unwrap(); + // /// + // /// // ROI (ddeg) must be within borders + // /// let roi_ddeg = LineString::new(vec![ + // /// coord! { x: -50.0, y: -23.0 }, + // /// coord! { x: -50.25, y: -23.1 }, + // /// coord! { x: -50.5, y: -23.2 }, + // /// ]); + // /// + // /// // you can double check that + // /// assert!(ionex.bounding_rect_degrees().contains(&roi_ddeg)); + // /// + // /// // Epoch must exist in the record + // /// let noon = Epoch::from_str("2022-01-02T12:00:00 UTC") + // /// .unwrap(); + // /// + // /// + // /// ``` + // pub fn spatial_area_interpolation( + // &self, + // area: &LineString, + // epoch: Epoch, + // ) -> BTreeMap { + // let mut values = BTreeMap::new(); + + // let fixed_altitude_km = self.header.grid.altitude.start; + + // // for all requested coordinates + // for point in area.points() { + // let geo = Geometry::Point(point); + + // let key = Key::from_decimal_degrees_km(epoch, point.y(), point.x(), fixed_altitude_km); + + // for cell in self.synchronous_map_cell_iter(epoch) { + // if cell.contains(&geo) { + // let tec = cell.spatial_interpolation(point); + // values.insert(key, tec); + // } + // } + // } + + // values + // } + + // /// Interpolate TEC values for all discrete coordinates described by the following [LineString] + // /// (in decimal degrees), at specific point in time that does exist within this record. + // /// We will interpolate the two neighbouring data points, which is not feasible at day boundaries. + // pub fn temporal_spatial_area_interpolation( + // &self, + // area: &LineString, + // epoch: Epoch, + // ) -> BTreeMap { + // let mut values = BTreeMap::new(); + + // let (min_t, max_t) = ( + // (epoch - self.header.sampling_period), + // (epoch + self.header.sampling_period), + // ); + + // let (min_t, max_t) = ( + // min_t.round(self.header.sampling_period), + // max_t.round(self.header.sampling_period), + // ); + + // // for all requested coordinates + // for point in area.points() { + // let geo = Geometry::Point(point); + + // for cell_t0 in self.synchronous_map_cell_iter(min_t) { + // if cell_t0.contains(&geo) { + // for cell_t1 in self.synchronous_map_cell_iter(max_t) { + // if cell_t1.contains(&geo) { + // if let Some(interpolated) = + // cell_t0.temporal_spatial_interpolation(epoch, point, &cell_t1) + // { + // let key = Key { + // epoch, + // coordinates: QuantizedCoordinates::from_decimal_degrees( + // point.y(), + // point.x(), + // self.header.grid.altitude.start, + // ), + // }; + + // values.insert(key, interpolated); + // } + // } + // } + // } + // } + // } + + // values + // } } /// Merge two [IONEX] structures into one. From 85e08e6f163f496d14d546ad82072b7a611d7e6f Mon Sep 17 00:00:00 2001 From: "Guillaume W. Bres" Date: Sun, 2 Nov 2025 18:05:02 +0100 Subject: [PATCH 17/25] upgrading API Signed-off-by: Guillaume W. Bres --- src/cell.rs | 42 +++++++++++++++++++++++++++--------------- 1 file changed, 27 insertions(+), 15 deletions(-) diff --git a/src/cell.rs b/src/cell.rs index 883dc52..806c3b3 100644 --- a/src/cell.rs +++ b/src/cell.rs @@ -350,7 +350,7 @@ impl MapCell { /// /// // central point /// let center = Point::new(0.5, 0.5); - /// let tec = cell.spatial_interpolation(center); + /// let tec = cell.spatial_tec_interp(center); /// assert_eq!(tec.tecu(), 1.0); /// ``` /// @@ -371,19 +371,19 @@ impl MapCell { /// let cell = MapCell::from_unitary_tec(t0, gradient.0, gradient.1, gradient.2, gradient.3); /// /// // central point - /// let tec = cell.spatial_interpolation(Point::new(0.5, 0.5)); + /// let tec = cell.spatial_tec_interp(Point::new(0.5, 0.5)); /// assert_eq!(tec.tecu(), 0.25); /// /// // SW boundary - /// let tec = cell.spatial_interpolation(Point::new(0.0, 0.0)); + /// let tec = cell.spatial_tec_interp(Point::new(0.0, 0.0)); /// assert_eq!(tec.tecu(), 1.0); /// /// // SWern point - /// let tec = cell.spatial_interpolation(Point::new(0.1, 0.1)); + /// let tec = cell.spatial_tec_interp(Point::new(0.1, 0.1)); /// assert_eq!(tec.tecu(), 0.81); /// /// // SWwern point - /// let tec = cell.spatial_interpolation(Point::new(0.01, 0.01)); + /// let tec = cell.spatial_tec_interp(Point::new(0.01, 0.01)); /// assert_eq!(tec.tecu(), 0.9801); /// ``` pub fn spatial_tec_interp(&self, point: Point) -> Result { @@ -545,19 +545,19 @@ impl MapCell { /// let cell1 = MapCell::from_unitary_tec(t1, one_tec, one_tec, one_tec, one_tec); /// /// // verify central point value - /// let central_tec0 = cell0.spatial_interpolation(center); + /// let central_tec0 = cell0.spatial_tec_interp(center); /// assert_eq!(central_tec0.tecu(), 1.0); /// /// // verify central point value - /// let central_tec1 = cell1.spatial_interpolation(center); + /// let central_tec1 = cell1.spatial_tec_interp(center); /// assert_eq!(central_tec1.tecu(), 1.0); /// /// // spatial + temporal interpolation /// // outside sampling interval - /// assert!(cell0.temporal_spatial_interpolation(t_nok, center, &cell1).is_none()); + /// assert!(cell0.temporal_spatial_tec_interp(t_nok, center, &cell1).is_none()); /// /// // spatial + temporal interpolation - /// let tec = cell0.temporal_spatial_interpolation(t_ok, center, &cell1).unwrap(); + /// let tec = cell0.temporal_spatial_tec_interp(t_ok, center, &cell1).unwrap(); /// assert_eq!(tec.tecu(), 1.0); /// ``` pub fn temporal_spatial_tec_interp( @@ -629,7 +629,10 @@ mod test { assert!(!cell.contains(&Geometry::Point(outside))); assert_eq!(center, cell.center()); - let interpolated = cell.spatial_interpolation(center); + let interpolated = cell.spatial_tec_interp(center).unwrap_or_else(|e| { + panic!("should have been feasible: {}", e); + }); + assert_eq!(interpolated.tecu(), 1.0); } @@ -657,7 +660,14 @@ mod test { (0.0, 0.0, 1.0), ] { let point = Point::new(x_deg, y_deg); - let interpolated = cell.spatial_interpolation(point).tecu(); + + let interpolated = cell + .spatial_tec_interp(point) + .unwrap_or_else(|e| { + panic!("should have been feasible! {}", e); + }) + .tecu(); + assert_eq!(interpolated, tecu, "failed at (x={}, y={})", x_deg, y_deg); } } @@ -700,14 +710,16 @@ mod test { assert!( cell0 - .temporal_spatial_interpolation(t_nok, center, &cell1) - .is_none(), + .temporal_spatial_tec_interp(t_nok, center, &cell1) + .is_err(), "interpolation is temporally incorrect!" ); let tec = cell0 - .temporal_spatial_interpolation(t_ok, center, &cell1) - .unwrap(); + .temporal_spatial_tec_interp(t_ok, center, &cell1) + .unwrap_or_else(|e| { + panic!("should have been feasible! {}", e); + }); assert_eq!(tec.tecu(), 1.0); } From 601ae7804d91b6514d1b7f863d8f4978d6d63160 Mon Sep 17 00:00:00 2001 From: "Guillaume W. Bres" Date: Sun, 2 Nov 2025 18:46:43 +0100 Subject: [PATCH 18/25] upgrading API Signed-off-by: Guillaume W. Bres --- src/cell.rs | 4 +- src/error.rs | 7 +++- src/header/qc.rs | 7 ++++ src/lib.rs | 99 +++++++++++++++++++++++++++++++++++++++++++----- src/tests/qc.rs | 41 +++++++++++++++++--- 5 files changed, 141 insertions(+), 17 deletions(-) diff --git a/src/cell.rs b/src/cell.rs index 806c3b3..1d4e6e7 100644 --- a/src/cell.rs +++ b/src/cell.rs @@ -271,6 +271,8 @@ impl MapCell { } /// Returns borders of this [MapCell] expressed as a [Rect]angle, in decimal degrees. + /// This is a direct conversion of this [MapCell] in terms of spatial dimensions, + /// discarding the associated TEC values. pub fn bounding_rect_degrees(&self) -> Rect { Rect::new(self.south_west.point, self.north_east.point) } @@ -388,7 +390,7 @@ impl MapCell { /// ``` pub fn spatial_tec_interp(&self, point: Point) -> Result { if !self.contains(&Geometry::Point(point)) { - return Err(Error::OutsideBoundaries); + return Err(Error::OutsideSpatialBoundaries); } let (latitude_span, longitude_span) = self.latitude_longitude_span_degrees(); diff --git a/src/error.rs b/src/error.rs index 871b6c5..ffe51bd 100644 --- a/src/error.rs +++ b/src/error.rs @@ -99,8 +99,11 @@ pub enum Error { #[error("undefined ROI exterior boundaries")] UndefinedBoundaries, - #[error("provided coordinates are outside boundaries")] - OutsideBoundaries, + #[error("coordinates are outside spatial boundaries")] + OutsideSpatialBoundaries, + + #[error("temporal coordinates outside this temporal axis")] + OutsideTemporalBoundaries, #[error("both regions do not describe the same spatial ROI")] SpatialMismatch, diff --git a/src/header/qc.rs b/src/header/qc.rs index 70056ce..c755780 100644 --- a/src/header/qc.rs +++ b/src/header/qc.rs @@ -85,6 +85,13 @@ impl Merge for Header { } } + // insert special comment + let merge_comment = "FILE MERGE".to_string(); + + if !self.comments.contains(&merge_comment) { + self.comments.push(merge_comment.clone()); + } + Ok(()) } } diff --git a/src/lib.rs b/src/lib.rs index 37a8ba7..c8c5e27 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -52,7 +52,7 @@ use std::{ use itertools::Itertools; -use geo::{coord, BoundingRect, Geometry, LineString, Point, Polygon, Rect}; +use geo::{coord, BoundingRect, Contains, Geometry, LineString, Point, Polygon, Rect}; #[cfg(feature = "flate2")] use flate2::{read::GzDecoder, write::GzEncoder, Compression as GzCompression}; @@ -939,16 +939,42 @@ impl IONEX { /// Prefer [Self::unitary_roi_at] in this case, to not bother designing the unit [Rect]. /// - or the closest fitting [MapCell] (wrapping boundaries) that we could fit pub fn roi_at(&self, epoch: Epoch, roi: Rect) -> Result { - if epoch < self.header.epoch_of_first_map { - return Err(Error::TemporalMismatch); + // determine whether this is within the temporal axis + if epoch < self.header.epoch_of_first_map || epoch > self.header.epoch_of_last_map { + return Err(Error::OutsideTemporalBoundaries); } - if epoch > self.header.epoch_of_last_map { - return Err(Error::TemporalMismatch); + // determine this is within the map + if !self.bounding_rect_degrees().contains(&Geometry::Rect(roi)) { + return Err(Error::OutsideSpatialBoundaries); } - // for each corner of the ROI, determine the best suited - // map cell for spatial interpolation + // determine whether we need temporal interpolation or not + let mut needs_temporal_interp = true; + let mut t = self.header.epoch_of_first_map; + + while t < self.header.epoch_of_last_map { + if t == epoch { + needs_temporal_interp = false; + break; + } + + t += self.header.sampling_period; + } + + // we don't need spatial interpolation if the ROI height and width + // are a perfect multiple of the grid dimensions + let (width, height) = (roi.width(), roi.height()); + + let mut needs_lat_interpolation = + height.rem_euclid(self.header.grid.latitude.spacing) == 0.0; + let mut needs_long_interpolation = + width.rem_euclid(self.header.grid.longitude.spacing) == 0.0; + + let needs_spatial_interpolation = needs_lat_interpolation || needs_long_interpolation; + if needs_spatial_interpolation { + panic!("not supported yet"); + } // TODO Err(Error::TemporalMismatch) @@ -971,8 +997,63 @@ impl IONEX { coordinates: Point, card: Cardinal, ) -> Result { - // TODO - Err(Error::TemporalMismatch) + // builds the unitary ROI (as Rect) matching those coordinates and cardinal + let (lat_pairs, long_pairs) = ( + self.header.grid.latitude.quantize().tuple_windows(), + self.header.grid.longitude.quantize().tuple_windows(), + ); + + let (lat_spacing, long_spacing) = ( + self.header.grid.latitude.spacing, + self.header.grid.longitude.spacing, + ); + + for (lat1, lat2) in lat_pairs { + for (long1, long2) in long_pairs.clone() { + let (lat1, lat2, long1, long2) = ( + lat1.real_value(), + lat2.real_value(), + long1.real_value(), + long2.real_value(), + ); + match card { + Cardinal::NorthEast => { + if lat1 == coordinates.y() && long1 == coordinates.x() { + let roi = + Rect::new(coord!(x: long1, y: lat1), coord!(x: long2, y: lat2)); + + return self.roi_at(epoch, roi); + } + }, + Cardinal::SouthEast => { + if lat1 == coordinates.y() && long1 == coordinates.x() { + let roi = + Rect::new(coord!(x: long1, y: lat1), coord!(x: long2, y: lat2)); + + return self.roi_at(epoch, roi); + } + }, + Cardinal::SouthWest => { + if lat1 == coordinates.y() && long1 == coordinates.x() { + let roi = + Rect::new(coord!(x: long1, y: lat1), coord!(x: long2, y: lat2)); + + return self.roi_at(epoch, roi); + } + }, + Cardinal::NorthWest => { + if lat1 == coordinates.y() && long1 == coordinates.x() { + let roi = + Rect::new(coord!(x: long1, y: lat1), coord!(x: long2, y: lat2)); + + return self.roi_at(epoch, roi); + } + }, + } + } + } + + Err(Error::OutsideSpatialBoundaries) } /// Obtain the best suited [MapCell] spatially wrapping this Geometry that contains following [Geometry]. diff --git a/src/tests/qc.rs b/src/tests/qc.rs index 75f2d4d..09522c7 100644 --- a/src/tests/qc.rs +++ b/src/tests/qc.rs @@ -23,6 +23,8 @@ fn v1_self_merge() { panic!("self::merge(self) should be feasible! {}", e); }); + assert!(merged.is_merged(), "merged file not declared as merged!"); + // reciprocal let fd = File::create("merged.txt").unwrap(); let mut writer = BufWriter::new(fd); @@ -32,11 +34,18 @@ fn v1_self_merge() { }); // parse back - let parsed = IONEX::from_file("merged.txt").unwrap_or_else(|e| { + let mut parsed = IONEX::from_file("merged.txt").unwrap_or_else(|e| { panic!("failed to parse back CKMG V1: {}", e); }); - // full reciprocity + assert!(parsed.is_merged(), "merged file not declared as merged!"); + + // remove special comment and run strict equality test + parsed + .header + .comments + .retain(|comment| !comment.eq("FILE MERGE")); + generic_comparison(&parsed, &ionex); } @@ -55,7 +64,19 @@ fn v1_files_merge() { panic!("failed to merge both files: {}", e); }); - // verifications + // testbench + assert!(merged.is_merged(), "merged file not declared as merged!"); + assert_eq!( + merged.header.epoch_of_first_map.to_string(), + "2021-01-09T00:00:00 UTC" + ); + assert_eq!( + merged.header.epoch_of_last_map.to_string(), + "2022-01-03T00:00:00 UTC" + ); + assert!(merged.is_worldwide_map()); + + // Dump to file let fd = File::create("merged.txt").unwrap(); let mut writer = BufWriter::new(fd); @@ -64,9 +85,19 @@ fn v1_files_merge() { }); // parse back - let parsed = IONEX::from_file("ckmg-v1.txt").unwrap_or_else(|e| { + let parsed = IONEX::from_file("merged.txt").unwrap_or_else(|e| { panic!("failed to parse back CKMG V1: {}", e); }); - // TODO + // testbench + assert!(parsed.is_merged(), "merged file not declared as merged!"); + assert_eq!( + parsed.header.epoch_of_first_map.to_string(), + "2021-01-09T00:00:00 UTC" + ); + assert_eq!( + parsed.header.epoch_of_last_map.to_string(), + "2022-01-03T00:00:00 UTC" + ); + assert!(parsed.is_worldwide_map()); } From f010550387f7cb362b5ee454f415a9e245854629 Mon Sep 17 00:00:00 2001 From: "Guillaume W. Bres" Date: Sun, 2 Nov 2025 19:53:36 +0100 Subject: [PATCH 19/25] upgrading API Signed-off-by: Guillaume W. Bres --- src/cell.rs | 150 +++++++++++++++++++++++++++++++++++++++++++++++++++- src/lib.rs | 66 ++++++++++++++++------- 2 files changed, 196 insertions(+), 20 deletions(-) diff --git a/src/cell.rs b/src/cell.rs index 1d4e6e7..b4337d8 100644 --- a/src/cell.rs +++ b/src/cell.rs @@ -1,6 +1,9 @@ use geo::{Contains, GeodesicArea, Geometry, Point, Rect}; -use crate::prelude::{Epoch, Error, TEC}; +use crate::{ + prelude::{Epoch, Error, TEC}, + rectangle_to_cardinals, +}; #[derive(Debug, Copy, Clone, PartialEq, PartialOrd)] #[cfg_attr(feature = "serde", derive(Serialize, Deserialize))] @@ -410,7 +413,150 @@ impl MapCell { Ok(TEC::from_tecu(tecu)) } - /// Merges two neighboring [MapCell]s forming a new [MapCell] with best precision. + /// Determines the northeastern cell amongst a grouping of 4 + pub fn northeasternmost_cell4(cell1: &Self, cell2: &Self, cell3: &Self, cell4: &Self) -> Self { + let min_x = cell1 + .point + .x() + .min(cell2.point.x()) + .min(cell3.point.x()) + .min(cell4.point.x()); + let min_y = cell1 + .point + .y() + .min(cell2.point.y()) + .min(cell3.point.y()) + .min(cell4.point.y()); + + if cell1.point.x() == min_x && cell1.point.y() == min_y { + *cell_1 + } else if cell2.point.x() == min_x && cell2.point.y() == min_y { + *cell_2 + } else if cell3.point.x() == min_x && cell3.point.y() == min_y { + *cell_3 + } else { + *cell_4 + } + } + + /// Determines the northwestern cell amongst a grouping of 4 + pub fn northeasternmost_cell4(cell1: &Self, cell2: &Self, cell3: &Self, cell4: &Self) -> Self { + let min_x = cell1 + .point + .x() + .min(cell2.point.x()) + .min(cell3.point.x()) + .min(cell4.point.x()); + let min_y = cell1 + .point + .y() + .min(cell2.point.y()) + .min(cell3.point.y()) + .min(cell4.point.y()); + + if cell1.point.x() == min_x && cell1.point.y() == min_y { + *cell_1 + } else if cell2.point.x() == min_x && cell2.point.y() == min_y { + *cell_2 + } else if cell3.point.x() == min_x && cell3.point.y() == min_y { + *cell_3 + } else { + *cell_4 + } + } + + /// Determines the southwestern cell amongst a grouping of 4 + pub fn southwesternmost_cell4(cell1: &Self, cell2: &Self, cell3: &Self, cell4: &Self) -> Self { + let min_x = cell1 + .point + .x() + .min(cell2.point.x()) + .min(cell3.point.x()) + .min(cell4.point.x()); + let min_y = cell1 + .point + .y() + .min(cell2.point.y()) + .min(cell3.point.y()) + .min(cell4.point.y()); + + if cell1.point.x() == min_x && cell1.point.y() == min_y { + *cell_1 + } else if cell2.point.x() == min_x && cell2.point.y() == min_y { + *cell_2 + } else if cell3.point.x() == min_x && cell3.point.y() == min_y { + *cell_3 + } else { + *cell_4 + } + } + + /// Determines the southwestern cell amongst a grouping of 4 + pub fn southwesternmost_cell4(cell1: &Self, cell2: &Self, cell3: &Self, cell4: &Self) -> Self { + let min_x = cell1 + .point + .x() + .min(cell2.point.x()) + .min(cell3.point.x()) + .min(cell4.point.x()); + let min_y = cell1 + .point + .y() + .min(cell2.point.y()) + .min(cell3.point.y()) + .min(cell4.point.y()); + + if cell1.point.x() == min_x && cell1.point.y() == min_y { + *cell_1 + } else if cell2.point.x() == min_x && cell2.point.y() == min_y { + *cell_2 + } else if cell3.point.x() == min_x && cell3.point.y() == min_y { + *cell_3 + } else { + *cell_4 + } + } + + /// Interpolate the TEC at 4 points described by provided [Rect]angle, returning a new [MapCell]. + /// We use 4 neighboring [MapCell]s to apply the interpolation equation with the highest precision + /// on each corner of the [Rect]angle. This is to be used when the ROI does not align with the grid space. + pub fn interpolate_at_grouping4( + &self, + roi: Rect, + neighbor_1: &Self, + neighbor_2: &Self, + neighbor_3: &Self, + ) -> Result { + let ( + (roi_ne_lat, roi_ne_long), + (roi_se_lat, roi_se_long), + (roi_sw_lat, roi_sw_long), + (roi_nw_lat, roi_nw_long), + ) = rectangle_to_cardinals(roi); + + // determines NE, SE, NW, NE cells conveniently + // so we tolerate a random order (but they need to be neighboring cells) + let ne_cell = Self::northeasternmost_cell4(self, neighbor_1, neighbor_2, neighbor_3); + let se_cell = Self::southasternmost_cell4(self, neighbor_1, neighbor_2, neighbor_3); + let nw_cell = Self::northwesternmost_cell4(self, neighbor_1, neighbor_2, neighbor_3); + let ne_cell = Self::northeasternmost_cell4(self, neighbor_1, neighbor_2, neighbor_3); + + // verifies they are all neighboring cells + if !nw_cell.is_western_neighbor(ne_cell) { + return Err(); + } + if !ne_cell.is_northern_neighbor(se_cell) { + return Err(); + } + if !se_cell.is_eastern_neighbor(sw_cell) { + return Err(); + } + if !nw_cell.is_northern_neighbor(sw_cell) { + return Err(); + } + } + + /// Merges two neighboring [MapCell]s forming a new (upscaled) [MapCell]. /// Both cells must be synchronous. pub fn merge_neighbors(&self, rhs: &Self) -> Result { if !self.temporal_match(rhs) { diff --git a/src/lib.rs b/src/lib.rs index c8c5e27..938b5a8 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -115,6 +115,20 @@ fn div_ceil(value: usize, divider: usize) -> usize { } } +/// Converts a geo [Rect]angle to NE, SE, SW, NW (latitude, longitude) quadruplets +pub(crate) fn rectangle_to_cardinals( + rect: Rect, +) -> ((f64, f64), (f64, f64), (f64, f64), (f64, f64)) { + let (min, width, height) = (rect.min(), rect.width(), rect.height()); + let (x0, y0) = (min.x, min.y); + ( + (x0, y0), + (x0 + width, y0), + (x0 + width, y0 + height), + (x0, y0 + height), + ) +} + /// macro to format one header line or a comment pub(crate) fn fmt_ionex(content: &str, marker: &str) -> String { if content.len() < 60 { @@ -1234,7 +1248,7 @@ impl gnss_qc_traits::Merge for IONEX { #[cfg(test)] mod test { - use crate::{div_ceil, fmt_comment, prelude::*}; + use crate::{div_ceil, fmt_comment, prelude::*, rectangle_to_cardinals}; #[test] fn fmt_comments_singleline() { @@ -1275,22 +1289,38 @@ mod test { } #[test] - #[ignore] - fn regional_ionex() { - let ionex = IONEX::from_gzip_file("data/IONEX/V1/CKMG0020.22I.gz").unwrap_or_else(|e| { - panic!("Failed to parse CKMG0020: {}", e); - }); - - let roi = Rect::new(coord!(x: -180.0, y: -85.0), coord!(x: 180.0, y: -82.5)); - - let regional = ionex.to_regional_ionex(roi.into()).unwrap(); - - // dump - regional.to_file("region.txt").unwrap(); - - // parse - let parsed = IONEX::from_file("region.txt").unwrap_or_else(|e| { - panic!("Failed to parse region.txt: {}", e); - }); + fn rect_to_cardinals() { + for (rect, ((lat11, long11), (lat12, long12), (lat21, long21), (lat22, long22))) in [ + ( + Rect::new(coord!(x: -30.0, y: -30.0), coord!(x: 30.0, y: 30.0)), + ( + (-30.0, -30.0), + (-30.0, -30.0), + (-30.0, -30.0), + (-30.0, -30.0), + ), + ), + ( + Rect::new(coord!(x: -40.0, y: -40.0), coord!(x: 40.0, y: 50.0)), + ( + (-30.0, -30.0), + (-30.0, -30.0), + (-30.0, -30.0), + (-30.0, -30.0), + ), + ), + ] { + assert_eq!( + rectangle_to_cardinals(rect), + ( + (lat11, long11), + (lat12, long12), + (lat21, long21), + (lat22, long22) + ), + "failed for {:?}", + rect + ); + } } } From 77484354b29d8ee426a4e04e566a0ef74cbcebd5 Mon Sep 17 00:00:00 2001 From: "Guillaume W. Bres" Date: Mon, 10 Nov 2025 17:51:47 +0100 Subject: [PATCH 20/25] up Signed-off-by: Guillaume W. Bres --- src/cell.rs | 506 +++++++++++++++++++++++++++++++--------------------- src/lib.rs | 39 ++-- 2 files changed, 334 insertions(+), 211 deletions(-) diff --git a/src/cell.rs b/src/cell.rs index b4337d8..d01bb86 100644 --- a/src/cell.rs +++ b/src/cell.rs @@ -1,3 +1,5 @@ +use std::collections::HashMap; + use geo::{Contains, GeodesicArea, Geometry, Point, Rect}; use crate::{ @@ -5,29 +7,41 @@ use crate::{ rectangle_to_cardinals, }; -#[derive(Debug, Copy, Clone, PartialEq, PartialOrd)] +#[derive(Debug, Copy, Clone, PartialEq, Eq, Hash, PartialOrd)] #[cfg_attr(feature = "serde", derive(Serialize, Deserialize))] pub enum Cardinal { /// NE [Cardinal] NorthEast, + /// N [Cardinal] + North, + /// NW [Cardinal] NorthWest, - /// SE [Cardinal] - SouthEast, + /// W [Cardinal] + West, /// SW [Cardinal] SouthWest, + + /// S [Cardinal] + South, + + /// SE [Cardinal] + SouthEast, + + /// E [Cardinal] + East, } #[derive(Debug, Default, Clone, Copy, PartialEq)] -pub struct MapPoint { - /// [Point] - pub point: Point, - +pub struct TecPoint { /// TEC pub tec: TEC, + + /// [Point] + pub point: Point, } /// [MapCell] describing a 4 corner region that we can interpolate. @@ -38,17 +52,88 @@ pub struct MapCell { /// Epoch of observation pub epoch: Epoch, - /// North East [MapPoint] - pub north_east: MapPoint, + /// North East [TecPoint] + pub north_east: TecPoint, + + /// North West [TecPoint] + pub north_west: TecPoint, + + /// South East [TecPoint] + pub south_east: TecPoint, + + /// South West [TecPoint] + pub south_west: TecPoint, +} + +/// A structure holding 9 synchronous neighboring [MapCells] +#[derive(Debug, Default, Clone, PartialEq)] +pub struct Cell9 { + /// The central [MapCell] + pub central: MapCell, + + /// 8 Neighboring [MapCells] + pub neighbors: HashMap, +} + +impl Cell9 { + /// Returns true if this [Cell9] definition is "complete" + pub(crate) fn is_complete(&self) -> bool { + self.neighbors.len() == 8 + } + + /// Builds a new updated [Cell9] with updated central cell. + pub fn with_central_cell(mut self, center: MapCell) -> Self { + self.central = center; + self + } - /// North West [MapPoint] - pub north_west: MapPoint, + /// Builds a new [Cell9] conveniently from 9 [MapCells], correctly + /// electing the central ROI, and the 8 neighbors. + /// Returns None if no central ROI exists or could not determine 8 neighbors. + pub fn from_slice(cells: [MapCell; 9]) -> Option { + for i in 0..9 { + // determine whether the ith cell is the potential center + let mut all_neighbors = true; + // must be neighbor with all other 8 ROIs + for j in 0..9 { + if j != i { + if !cells[i].is_neighbor(&cells[j]) { + all_neighbors = false; + break; + } + } + } + + if all_neighbors { + let mut ret = Self::default().with_central_cell(cells[i]); + + // i is the central ROI, order other ROIs correctly & exit + for j in 0..9 { + if j != i { + if cells[j].is_northwestern_neighbor(&cells[i]) { + ret.neighbors.insert(Cardinal::NorthWest, cells[j]); + } else if cells[j].is_northeastern_neighbor(&cells[i]) { + ret.neighbors.insert(Cardinal::NorthEast, cells[j]); + } else if cells[j].is_northern_neighbor(&cells[i]) { + ret.neighbors.insert(Cardinal::North, cells[j]); + } else if cells[j].is_southwestern_neighbor(&cells[i]) { + ret.neighbors.insert(Cardinal::SouthWest, cells[j]); + } else if cells[j].is_southeastern_neighbor(&cells[i]) { + ret.neighbors.insert(Cardinal::SouthEast, cells[j]); + } else if cells[j].is_southern_neighbor(&cells[i]) { + ret.neighbors.insert(Cardinal::South, cells[j]); + } + } + } - /// South East [MapPoint] - pub south_east: MapPoint, + if ret.is_complete() { + return Some(ret); + } + } + } - /// South West [MapPoint] - pub south_west: MapPoint, + None + } } impl MapCell { @@ -67,19 +152,19 @@ impl MapCell { ) -> Self { Self { epoch, - north_east: MapPoint { + north_east: TecPoint { point: Point::new(northeast_ddeg.0, northeast_ddeg.1), tec: northeast_tec, }, - north_west: MapPoint { + north_west: TecPoint { point: Point::new(northwest_ddeg.0, northwest_ddeg.1), tec: northwest_tec, }, - south_east: MapPoint { + south_east: TecPoint { point: Point::new(southeast_ddeg.0, southeast_ddeg.1), tec: southeast_tec, }, - south_west: MapPoint { + south_west: TecPoint { point: Point::new(southwest_ddeg.0, southwest_ddeg.1), tec: southwest_tec, }, @@ -101,32 +186,32 @@ impl MapCell { ) -> Self { Self { epoch, - north_east: MapPoint { + north_east: TecPoint { point: Point::new(northeast_rad.0.to_degrees(), northeast_rad.1.to_degrees()), tec: northeast_tec, }, - north_west: MapPoint { + north_west: TecPoint { point: Point::new(northwest_rad.0.to_degrees(), northwest_rad.1.to_degrees()), tec: northwest_tec, }, - south_east: MapPoint { + south_east: TecPoint { point: Point::new(southeast_rad.0.to_degrees(), southeast_rad.1.to_degrees()), tec: southeast_tec, }, - south_west: MapPoint { + south_west: TecPoint { point: Point::new(southwest_rad.0.to_degrees(), southwest_rad.1.to_degrees()), tec: southwest_tec, }, } } - /// Define a new [MapCell] from all 4 [MapPoint]s describing each corner at this [Epoch]. + /// Define a new [MapCell] from all 4 [TecPoint]s describing each corner at this [Epoch]. pub fn from_cardinal_points( epoch: Epoch, - north_east: MapPoint, - north_west: MapPoint, - south_east: MapPoint, - south_west: MapPoint, + north_east: TecPoint, + north_west: TecPoint, + south_east: TecPoint, + south_west: TecPoint, ) -> Self { Self { epoch, @@ -157,25 +242,41 @@ impl MapCell { self.epoch == rhs.epoch } - /// Returns true if self is the northern neighbor of provided rhs [MapCell]. + /// Returns true if self is the northern neighbor of provided (rhs) [MapCell]. pub fn is_northern_neighbor(&self, rhs: &Self) -> bool { rhs.north_east.point == self.south_east.point && rhs.north_west.point == self.south_west.point } - /// Returns true if self is the southern neighbor of provided rhs [MapCell]. + pub fn is_northwestern_neighbor(&self, rhs: &Self) -> bool { + false // TODO + } + + pub fn is_northeastern_neighbor(&self, rhs: &Self) -> bool { + false // TODO + } + + /// Returns true if self is the southern neighbor of provided (rhs) [MapCell]. pub fn is_southern_neighbor(&self, rhs: &Self) -> bool { rhs.south_east.point == self.north_east.point && rhs.south_west.point == self.north_west.point } - /// Returns true if self is the easthern neighbor of provided rhs [MapCell]. + pub fn is_southeastern_neighbor(&self, rhs: &Self) -> bool { + false // TODO + } + + pub fn is_southwestern_neighbor(&self, rhs: &Self) -> bool { + false // TODO + } + + /// Returns true if self is the easthern neighbor of provided (rhs) [MapCell]. pub fn is_eastern_neighbor(&self, rhs: &Self) -> bool { rhs.north_west.point == self.north_east.point && rhs.south_west.point == self.south_east.point } - /// Returns true if self is the westhern neighbor of provided rhs [MapCell]. + /// Returns true if self is the westhern neighbor of provided (rhs) [MapCell]. pub fn is_western_neighbor(&self, rhs: &Self) -> bool { rhs.north_east.point == self.north_west.point && rhs.south_east.point == self.south_west.point @@ -184,14 +285,19 @@ impl MapCell { /// Returns true if both cells are neighbors, meaning, they share two corners. pub fn is_neighbor(&self, rhs: &Self) -> bool { self.is_northern_neighbor(rhs) + || self.is_northwestern_neighbor(rhs) + || self.is_northeastern_neighbor(rhs) || self.is_southern_neighbor(rhs) + || self.is_southwestern_neighbor(rhs) + || self.is_southeastern_neighbor(rhs) || self.is_western_neighbor(rhs) || self.is_eastern_neighbor(rhs) } - /// Returns true if this [MapCell] totally contains (wrapps) rhs cell. - /// Meaning, rhs is fully contained within Self. - pub fn contains_entirely(&self, rhs: &Self) -> bool { + /// Returns true if this [MapCell] contains (wrapps) entirely the spatial ROI + /// described by the provided (rhs) [MapCell]. + /// Meaning, rhs is fully contained within self. + pub fn wrapps_entirely(&self, rhs: &Self) -> bool { self.contains(&Geometry::Rect(rhs.bounding_rect_degrees())) } @@ -413,173 +519,173 @@ impl MapCell { Ok(TEC::from_tecu(tecu)) } - /// Determines the northeastern cell amongst a grouping of 4 - pub fn northeasternmost_cell4(cell1: &Self, cell2: &Self, cell3: &Self, cell4: &Self) -> Self { - let min_x = cell1 - .point - .x() - .min(cell2.point.x()) - .min(cell3.point.x()) - .min(cell4.point.x()); - let min_y = cell1 - .point - .y() - .min(cell2.point.y()) - .min(cell3.point.y()) - .min(cell4.point.y()); - - if cell1.point.x() == min_x && cell1.point.y() == min_y { - *cell_1 - } else if cell2.point.x() == min_x && cell2.point.y() == min_y { - *cell_2 - } else if cell3.point.x() == min_x && cell3.point.y() == min_y { - *cell_3 - } else { - *cell_4 - } - } + // /// Determines the northeastern cell amongst a grouping of 4 + // pub fn northeasternmost_cell4(cell1: &Self, cell2: &Self, cell3: &Self, cell4: &Self) -> Self { + // let min_x = cell1 + // .point + // .x() + // .min(cell2.point.x()) + // .min(cell3.point.x()) + // .min(cell4.point.x()); + // let min_y = cell1 + // .point + // .y() + // .min(cell2.point.y()) + // .min(cell3.point.y()) + // .min(cell4.point.y()); + + // if cell1.point.x() == min_x && cell1.point.y() == min_y { + // *cell_1 + // } else if cell2.point.x() == min_x && cell2.point.y() == min_y { + // *cell_2 + // } else if cell3.point.x() == min_x && cell3.point.y() == min_y { + // *cell_3 + // } else { + // *cell_4 + // } + // } - /// Determines the northwestern cell amongst a grouping of 4 - pub fn northeasternmost_cell4(cell1: &Self, cell2: &Self, cell3: &Self, cell4: &Self) -> Self { - let min_x = cell1 - .point - .x() - .min(cell2.point.x()) - .min(cell3.point.x()) - .min(cell4.point.x()); - let min_y = cell1 - .point - .y() - .min(cell2.point.y()) - .min(cell3.point.y()) - .min(cell4.point.y()); - - if cell1.point.x() == min_x && cell1.point.y() == min_y { - *cell_1 - } else if cell2.point.x() == min_x && cell2.point.y() == min_y { - *cell_2 - } else if cell3.point.x() == min_x && cell3.point.y() == min_y { - *cell_3 - } else { - *cell_4 - } - } + // /// Determines the northwestern cell amongst a grouping of 4 + // pub fn northeasternmost_cell4(cell1: &Self, cell2: &Self, cell3: &Self, cell4: &Self) -> Self { + // let min_x = cell1 + // .point + // .x() + // .min(cell2.point.x()) + // .min(cell3.point.x()) + // .min(cell4.point.x()); + // let min_y = cell1 + // .point + // .y() + // .min(cell2.point.y()) + // .min(cell3.point.y()) + // .min(cell4.point.y()); + + // if cell1.point.x() == min_x && cell1.point.y() == min_y { + // *cell_1 + // } else if cell2.point.x() == min_x && cell2.point.y() == min_y { + // *cell_2 + // } else if cell3.point.x() == min_x && cell3.point.y() == min_y { + // *cell_3 + // } else { + // *cell_4 + // } + // } - /// Determines the southwestern cell amongst a grouping of 4 - pub fn southwesternmost_cell4(cell1: &Self, cell2: &Self, cell3: &Self, cell4: &Self) -> Self { - let min_x = cell1 - .point - .x() - .min(cell2.point.x()) - .min(cell3.point.x()) - .min(cell4.point.x()); - let min_y = cell1 - .point - .y() - .min(cell2.point.y()) - .min(cell3.point.y()) - .min(cell4.point.y()); - - if cell1.point.x() == min_x && cell1.point.y() == min_y { - *cell_1 - } else if cell2.point.x() == min_x && cell2.point.y() == min_y { - *cell_2 - } else if cell3.point.x() == min_x && cell3.point.y() == min_y { - *cell_3 - } else { - *cell_4 - } - } + // /// Determines the southwestern cell amongst a grouping of 4 + // pub fn southwesternmost_cell4(cell1: &Self, cell2: &Self, cell3: &Self, cell4: &Self) -> Self { + // let min_x = cell1 + // .point + // .x() + // .min(cell2.point.x()) + // .min(cell3.point.x()) + // .min(cell4.point.x()); + // let min_y = cell1 + // .point + // .y() + // .min(cell2.point.y()) + // .min(cell3.point.y()) + // .min(cell4.point.y()); + + // if cell1.point.x() == min_x && cell1.point.y() == min_y { + // *cell_1 + // } else if cell2.point.x() == min_x && cell2.point.y() == min_y { + // *cell_2 + // } else if cell3.point.x() == min_x && cell3.point.y() == min_y { + // *cell_3 + // } else { + // *cell_4 + // } + // } - /// Determines the southwestern cell amongst a grouping of 4 - pub fn southwesternmost_cell4(cell1: &Self, cell2: &Self, cell3: &Self, cell4: &Self) -> Self { - let min_x = cell1 - .point - .x() - .min(cell2.point.x()) - .min(cell3.point.x()) - .min(cell4.point.x()); - let min_y = cell1 - .point - .y() - .min(cell2.point.y()) - .min(cell3.point.y()) - .min(cell4.point.y()); - - if cell1.point.x() == min_x && cell1.point.y() == min_y { - *cell_1 - } else if cell2.point.x() == min_x && cell2.point.y() == min_y { - *cell_2 - } else if cell3.point.x() == min_x && cell3.point.y() == min_y { - *cell_3 - } else { - *cell_4 - } - } + // /// Determines the southwestern cell amongst a grouping of 4 + // pub fn southwesternmost_cell4(cell1: &Self, cell2: &Self, cell3: &Self, cell4: &Self) -> Self { + // let min_x = cell1 + // .point + // .x() + // .min(cell2.point.x()) + // .min(cell3.point.x()) + // .min(cell4.point.x()); + // let min_y = cell1 + // .point + // .y() + // .min(cell2.point.y()) + // .min(cell3.point.y()) + // .min(cell4.point.y()); + + // if cell1.point.x() == min_x && cell1.point.y() == min_y { + // *cell_1 + // } else if cell2.point.x() == min_x && cell2.point.y() == min_y { + // *cell_2 + // } else if cell3.point.x() == min_x && cell3.point.y() == min_y { + // *cell_3 + // } else { + // *cell_4 + // } + // } - /// Interpolate the TEC at 4 points described by provided [Rect]angle, returning a new [MapCell]. - /// We use 4 neighboring [MapCell]s to apply the interpolation equation with the highest precision - /// on each corner of the [Rect]angle. This is to be used when the ROI does not align with the grid space. - pub fn interpolate_at_grouping4( - &self, - roi: Rect, - neighbor_1: &Self, - neighbor_2: &Self, - neighbor_3: &Self, - ) -> Result { - let ( - (roi_ne_lat, roi_ne_long), - (roi_se_lat, roi_se_long), - (roi_sw_lat, roi_sw_long), - (roi_nw_lat, roi_nw_long), - ) = rectangle_to_cardinals(roi); - - // determines NE, SE, NW, NE cells conveniently - // so we tolerate a random order (but they need to be neighboring cells) - let ne_cell = Self::northeasternmost_cell4(self, neighbor_1, neighbor_2, neighbor_3); - let se_cell = Self::southasternmost_cell4(self, neighbor_1, neighbor_2, neighbor_3); - let nw_cell = Self::northwesternmost_cell4(self, neighbor_1, neighbor_2, neighbor_3); - let ne_cell = Self::northeasternmost_cell4(self, neighbor_1, neighbor_2, neighbor_3); - - // verifies they are all neighboring cells - if !nw_cell.is_western_neighbor(ne_cell) { - return Err(); - } - if !ne_cell.is_northern_neighbor(se_cell) { - return Err(); - } - if !se_cell.is_eastern_neighbor(sw_cell) { - return Err(); - } - if !nw_cell.is_northern_neighbor(sw_cell) { - return Err(); - } - } + // /// Interpolate the TEC at 4 points described by provided [Rect]angle, returning a new [MapCell]. + // /// We use 4 neighboring [MapCell]s to apply the interpolation equation with the highest precision + // /// on each corner of the [Rect]angle. This is to be used when the ROI does not align with the grid space. + // pub fn interpolate_at_grouping4( + // &self, + // roi: Rect, + // neighbor_1: &Self, + // neighbor_2: &Self, + // neighbor_3: &Self, + // ) -> Result { + // let ( + // (roi_ne_lat, roi_ne_long), + // (roi_se_lat, roi_se_long), + // (roi_sw_lat, roi_sw_long), + // (roi_nw_lat, roi_nw_long), + // ) = rectangle_to_cardinals(roi); + + // // determines NE, SE, NW, NE cells conveniently + // // so we tolerate a random order (but they need to be neighboring cells) + // let ne_cell = Self::northeasternmost_cell4(self, neighbor_1, neighbor_2, neighbor_3); + // let se_cell = Self::southasternmost_cell4(self, neighbor_1, neighbor_2, neighbor_3); + // let nw_cell = Self::northwesternmost_cell4(self, neighbor_1, neighbor_2, neighbor_3); + // let ne_cell = Self::northeasternmost_cell4(self, neighbor_1, neighbor_2, neighbor_3); + + // // verifies they are all neighboring cells + // if !nw_cell.is_western_neighbor(ne_cell) { + // return Err(); + // } + // if !ne_cell.is_northern_neighbor(se_cell) { + // return Err(); + // } + // if !se_cell.is_eastern_neighbor(sw_cell) { + // return Err(); + // } + // if !nw_cell.is_northern_neighbor(sw_cell) { + // return Err(); + // } + // } - /// Merges two neighboring [MapCell]s forming a new (upscaled) [MapCell]. - /// Both cells must be synchronous. - pub fn merge_neighbors(&self, rhs: &Self) -> Result { - if !self.temporal_match(rhs) { - return Err(Error::TemporalMismatch); - } + // /// Merges two neighboring [MapCell]s forming a new (upscaled) [MapCell]. + // /// Both cells must be synchronous. + // pub fn merge_neighbors(&self, rhs: &Self) -> Result { + // if !self.temporal_match(rhs) { + // return Err(Error::TemporalMismatch); + // } - // 1: determine the matching corners. The center - // point of the matching line is to become the new center. - // 2: form a new cell, of the boundary rect, interpolate - // both TEC at the new center - - if self.is_northern_neighbor(rhs) { - Ok(Self::default()) // TODO - } else if self.is_southern_neighbor(rhs) { - Ok(Self::default()) // TODO - } else if self.is_western_neighbor(rhs) { - Ok(Self::default()) // TODO - } else if self.is_eastern_neighbor(rhs) { - Ok(Self::default()) // TODO - } else { - Err(Error::SpatialMismatch) - } - } + // // 1: determine the matching corners. The center + // // point of the matching line is to become the new center. + // // 2: form a new cell, of the boundary rect, interpolate + // // both TEC at the new center + + // if self.is_northern_neighbor(rhs) { + // Ok(Self::default()) // TODO + // } else if self.is_southern_neighbor(rhs) { + // Ok(Self::default()) // TODO + // } else if self.is_western_neighbor(rhs) { + // Ok(Self::default()) // TODO + // } else if self.is_eastern_neighbor(rhs) { + // Ok(Self::default()) // TODO + // } else { + // Err(Error::SpatialMismatch) + // } + // } // /// Interpolates two [MapCell]s that must describe the same area, // /// but a different point in time. @@ -642,25 +748,25 @@ impl MapCell { // Ok(Self { // epoch, - // north_east: MapPoint { + // north_east: TecPoint { // point: self.north_east.point, // tec: TEC::from_tecu( // num_1 * ne_1 /dt + num_2 * ne_2 /dt // ), // }, - // north_west: MapPoint { + // north_west: TecPoint { // point: self.north_west.point, // tec: TEC::from_tecu( // num_1 * nw_1 /dt + num_2 * nw_2 /dt // ), // }, - // south_east: MapPoint { + // south_east: TecPoint { // point: self.south_east.point, // tec: TEC::from_tecu( // num_1 * se_1 /dt + num_2 * se_2 /dt // ), // }, - // south_west: MapPoint { + // south_west: TecPoint { // point: self.south_west.point, // tec: TEC::from_tecu( // num_1 * sw_1 /dt + num_2 * sw_2 /dt diff --git a/src/lib.rs b/src/lib.rs index 938b5a8..944a305 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -60,7 +60,7 @@ use flate2::{read::GzDecoder, write::GzEncoder, Compression as GzCompression}; use hifitime::prelude::Epoch; use crate::{ - cell::{Cardinal, MapCell, MapPoint}, + cell::{Cardinal, Cell9, MapCell, TecPoint}, coordinates::QuantizedCoordinates, error::{Error, FormattingError, ParsingError}, file_attributes::{FileAttributes, Region}, @@ -76,7 +76,7 @@ pub mod prelude { // export pub use crate::{ bias::BiasSource, - cell::{Cardinal, MapCell}, + cell::{Cardinal, Cell9, MapCell}, error::{Error, FormattingError, ParsingError}, file_attributes::*, grid::{Axis, Grid}, @@ -686,15 +686,16 @@ impl IONEX { } /// Modify the grid spacing (quantization) while preserving the dimensions, - /// and interpolates the TEC values. + /// and interpolates the new TEC values. /// /// This only modify the grid quantization, the map dimensions are preserved. /// /// ## Input /// - axis: [Axis] to be resampled /// - factor: - /// - > 1.0: upscaling - /// - < 1.0: downscaling + /// - > 1.0: spatial upscaling. The grid precision increases, + /// we interpolate the TEC at new intermiediate coordinates. + /// - < 1.0: downscaling. The grid precision decreases. /// - 0.0: invalid pub fn spatially_resampled(&self, axis: Axis, factor: f64) -> Result { let mut s = self.clone(); @@ -768,6 +769,10 @@ impl IONEX { Box::new(self.record.map.keys().map(|k| k.epoch).unique().sorted()) } + /// Iterates this [IONEX] by a group of 9 neighboring [MapCell]s, + /// which is particularly convenient for accurate interpolation and upscaling. + pub fn iter_cell9(&self) -> Box + '_> {} + /// Modify the grid spacing (quantization) while preserving the dimensions, /// and interpolates the TEC values. /// @@ -802,8 +807,21 @@ impl IONEX { }, } + let lat_pairs = self.header.grid.latitude.quantize().tuple_windows(); + let long_pairs = self.header.grid.longitude.quantize().tuple_windows(); + + // applies the stretching + // synchronous spatial interpolation - for epoch in self.epoch_iter() {} + for epoch in self.epoch_iter() { + if factor > 1.0 { + // upscaling + } else { + // downscaling + } + } + + // saturate the possible overflows to the world map proj Ok(()) } @@ -867,7 +885,6 @@ impl IONEX { let timeseries = self.header.timeseries(); let lat_pairs = self.header.grid.latitude.quantize().tuple_windows(); - let long_pairs = self.header.grid.longitude.quantize().tuple_windows(); let fixed_altitude_km = self.header.grid.altitude.start; @@ -923,19 +940,19 @@ impl IONEX { Some(MapCell { epoch, - north_east: MapPoint { + north_east: TecPoint { tec: *northeast, point: Point::new(long1.real_value(), lat1.real_value()), }, - north_west: MapPoint { + north_west: TecPoint { tec: *northwest, point: Point::new(long2.real_value(), lat1.real_value()), }, - south_east: MapPoint { + south_east: TecPoint { tec: *southeast, point: Point::new(long1.real_value(), lat2.real_value()), }, - south_west: MapPoint { + south_west: TecPoint { tec: *southwest, point: Point::new(long2.real_value(), lat2.real_value()), }, From 2504a8dfea7c13aab16118db6d68dcf02214b1b6 Mon Sep 17 00:00:00 2001 From: "Guillaume W. Bres" Date: Mon, 10 Nov 2025 18:13:39 +0100 Subject: [PATCH 21/25] up Signed-off-by: Guillaume W. Bres --- src/cell.rs | 176 ++++++++++++++++++++++++++++++--------------------- src/cell9.rs | 88 ++++++++++++++++++++++++++ src/lib.rs | 4 +- 3 files changed, 196 insertions(+), 72 deletions(-) create mode 100644 src/cell9.rs diff --git a/src/cell.rs b/src/cell.rs index d01bb86..70b05b8 100644 --- a/src/cell.rs +++ b/src/cell.rs @@ -65,77 +65,6 @@ pub struct MapCell { pub south_west: TecPoint, } -/// A structure holding 9 synchronous neighboring [MapCells] -#[derive(Debug, Default, Clone, PartialEq)] -pub struct Cell9 { - /// The central [MapCell] - pub central: MapCell, - - /// 8 Neighboring [MapCells] - pub neighbors: HashMap, -} - -impl Cell9 { - /// Returns true if this [Cell9] definition is "complete" - pub(crate) fn is_complete(&self) -> bool { - self.neighbors.len() == 8 - } - - /// Builds a new updated [Cell9] with updated central cell. - pub fn with_central_cell(mut self, center: MapCell) -> Self { - self.central = center; - self - } - - /// Builds a new [Cell9] conveniently from 9 [MapCells], correctly - /// electing the central ROI, and the 8 neighbors. - /// Returns None if no central ROI exists or could not determine 8 neighbors. - pub fn from_slice(cells: [MapCell; 9]) -> Option { - for i in 0..9 { - // determine whether the ith cell is the potential center - let mut all_neighbors = true; - // must be neighbor with all other 8 ROIs - for j in 0..9 { - if j != i { - if !cells[i].is_neighbor(&cells[j]) { - all_neighbors = false; - break; - } - } - } - - if all_neighbors { - let mut ret = Self::default().with_central_cell(cells[i]); - - // i is the central ROI, order other ROIs correctly & exit - for j in 0..9 { - if j != i { - if cells[j].is_northwestern_neighbor(&cells[i]) { - ret.neighbors.insert(Cardinal::NorthWest, cells[j]); - } else if cells[j].is_northeastern_neighbor(&cells[i]) { - ret.neighbors.insert(Cardinal::NorthEast, cells[j]); - } else if cells[j].is_northern_neighbor(&cells[i]) { - ret.neighbors.insert(Cardinal::North, cells[j]); - } else if cells[j].is_southwestern_neighbor(&cells[i]) { - ret.neighbors.insert(Cardinal::SouthWest, cells[j]); - } else if cells[j].is_southeastern_neighbor(&cells[i]) { - ret.neighbors.insert(Cardinal::SouthEast, cells[j]); - } else if cells[j].is_southern_neighbor(&cells[i]) { - ret.neighbors.insert(Cardinal::South, cells[j]); - } - } - } - - if ret.is_complete() { - return Some(ret); - } - } - } - - None - } -} - impl MapCell { /// Define a new [MapCell] from 4 (latitude_ddeg, longitude_ddeg) cardinal tuples and /// associated TEC values. @@ -519,6 +448,111 @@ impl MapCell { Ok(TEC::from_tecu(tecu)) } + /// Returns a downscaled (resized minized) new ROI applying the interpolation equation + /// on each corners. + pub fn downscaled(&self, factor: f64) -> Result { + if !factor.is_normal() { + return Err(Error::InvalidStretchFactor); + } + + if factor > 1.0 { + return Err(Error::InvalidStretchFactor); + } + + // apply interpolation eq. at 4 coordinates + let (north_east, north_west, south_east, south_west) = ( + Point::new( + self.north_east.point.x() * factor, + self.north_east.point.y() * factor, + ), + Point::new( + self.north_west.point.x() * factor, + self.north_west.point.y() * factor, + ), + Point::new( + self.south_east.point.x() * factor, + self.south_east.point.y() * factor, + ), + Point::new( + self.south_west.point.x() * factor, + self.south_west.point.y() * factor, + ), + ); + + Ok(MapCell { + north_east: TecPoint { + tec: self.spatial_tec_interp(north_east), + point: north_east, + }, + north_west: TecPoint { + tec: self.spatial_tec_interp(north_west), + point: north_west, + }, + south_west: TecPoint { + tec: self.spatial_tec_interp(south_west), + point: south_west, + }, + south_east: TecPoint { + tec: self.spatial_tec_interp(south_east), + point: south_east, + }, + }) + } + + /// Returns a upscaled (resized upscaled) new ROI applying the interpolation equation + /// on each corners. This method should be used for very small stretch ratios. + /// For higher accuracy, you should work with a [Cell9] grouping and [Cell9::upscale] + /// in particular, which takes advantage of neighboring information on the map, to minimize + /// the error on the resulting ROI. + pub fn upscaled(&self, factor: f64) -> Result { + if !factor.is_normal() { + return Err(Error::InvalidStretchFactor); + } + + if factor < 1.0 { + return Err(Error::InvalidStretchFactor); + } + + // apply interpolation eq. at 4 coordinates + let (north_east, north_west, south_east, south_west) = ( + Point::new( + self.north_east.point.x() * factor, + self.north_east.point.y() * factor, + ), + Point::new( + self.north_west.point.x() * factor, + self.north_west.point.y() * factor, + ), + Point::new( + self.south_east.point.x() * factor, + self.south_east.point.y() * factor, + ), + Point::new( + self.south_west.point.x() * factor, + self.south_west.point.y() * factor, + ), + ); + + Ok(MapCell { + north_east: TecPoint { + tec: self.spatial_tec_interp(north_east), + point: north_east, + }, + north_west: TecPoint { + tec: self.spatial_tec_interp(north_west), + point: north_west, + }, + south_west: TecPoint { + tec: self.spatial_tec_interp(south_west), + point: south_west, + }, + south_east: TecPoint { + tec: self.spatial_tec_interp(south_east), + point: south_east, + }, + }) + } + // /// Determines the northeastern cell amongst a grouping of 4 // pub fn northeasternmost_cell4(cell1: &Self, cell2: &Self, cell3: &Self, cell4: &Self) -> Self { // let min_x = cell1 diff --git a/src/cell9.rs b/src/cell9.rs new file mode 100644 index 0000000..43e43eb --- /dev/null +++ b/src/cell9.rs @@ -0,0 +1,88 @@ +use std::collections::HashMap; + +use geo::{Contains, GeodesicArea, Geometry, Point, Rect}; + +use crate::prelude::{Cardinal, Epoch, Error, MapCell, TEC}; + +/// A structure holding 9 synchronous neighboring [MapCells] +#[derive(Debug, Default, Clone, PartialEq)] +pub struct Cell9 { + /// The central [MapCell] + pub central: MapCell, + + /// 8 Neighboring [MapCells] + pub neighbors: HashMap, +} + +impl Cell9 { + /// Returns true if this [Cell9] definition is "complete" + pub(crate) fn is_complete(&self) -> bool { + self.neighbors.len() == 8 + } + + /// Builds a new updated [Cell9] with updated central cell. + pub fn with_central_cell(mut self, center: MapCell) -> Self { + self.central = center; + self + } + + /// Builds a new [Cell9] conveniently from 9 [MapCells], correctly + /// electing the central ROI, and the 8 neighbors. + /// Returns None if no central ROI exists or could not determine 8 neighbors. + pub fn from_slice(cells: [MapCell; 9]) -> Option { + for i in 0..9 { + // determine whether the ith cell is the potential center + let mut all_neighbors = true; + // must be neighbor with all other 8 ROIs + for j in 0..9 { + if j != i { + if !cells[i].is_neighbor(&cells[j]) { + all_neighbors = false; + break; + } + } + } + + if all_neighbors { + let mut ret = Self::default().with_central_cell(cells[i]); + + // i is the central ROI, order other ROIs correctly & exit + for j in 0..9 { + if j != i { + if cells[j].is_northwestern_neighbor(&cells[i]) { + ret.neighbors.insert(Cardinal::NorthWest, cells[j]); + } else if cells[j].is_northeastern_neighbor(&cells[i]) { + ret.neighbors.insert(Cardinal::NorthEast, cells[j]); + } else if cells[j].is_northern_neighbor(&cells[i]) { + ret.neighbors.insert(Cardinal::North, cells[j]); + } else if cells[j].is_southwestern_neighbor(&cells[i]) { + ret.neighbors.insert(Cardinal::SouthWest, cells[j]); + } else if cells[j].is_southeastern_neighbor(&cells[i]) { + ret.neighbors.insert(Cardinal::SouthEast, cells[j]); + } else if cells[j].is_southern_neighbor(&cells[i]) { + ret.neighbors.insert(Cardinal::South, cells[j]); + } + } + } + + if ret.is_complete() { + return Some(ret); + } + } + } + + None + } + + /// Upscale or downscale the central ROI of this [Cell9] by a fractional positive scalar. + pub fn upscale(&self, factor: f64) -> Result { + // interpolate the 4 new central coordinates + // interpolate using the 4 neighbors + // take the average of both for each 4 coordinntes + } + + pub fn downscale(&self, factor: f64) -> Result { + // interpolate the central ROI + self.central.downscale(factor) + } +} diff --git a/src/lib.rs b/src/lib.rs index 944a305..3842c19 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -33,6 +33,7 @@ pub mod tec; pub mod version; mod cell; +mod cell9; mod coordinates; mod epoch; mod ionosphere; @@ -60,7 +61,8 @@ use flate2::{read::GzDecoder, write::GzEncoder, Compression as GzCompression}; use hifitime::prelude::Epoch; use crate::{ - cell::{Cardinal, Cell9, MapCell, TecPoint}, + cell::Cell9, + cell::{Cardinal, MapCell, TecPoint}, coordinates::QuantizedCoordinates, error::{Error, FormattingError, ParsingError}, file_attributes::{FileAttributes, Region}, From a3fc8d12c49020b061627b547a9c99a45e5c0ff6 Mon Sep 17 00:00:00 2001 From: "Guillaume W. Bres" Date: Mon, 10 Nov 2025 18:22:13 +0100 Subject: [PATCH 22/25] up Signed-off-by: Guillaume W. Bres --- src/lib.rs | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 3842c19..b8d4ba4 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -717,6 +717,8 @@ impl IONEX { /// - < 1.0: downscaling /// - 0.0: invalid pub fn spatial_stretching_mut(&mut self, axis: Axis, factor: f64) -> Result<(), Error> { + // the maximal factor supported for on iter is 2 + // Update the header specs match axis { Axis::Latitude => { @@ -761,7 +763,13 @@ impl IONEX { } // synchronous spatial interpolation - for epoch in self.epoch_iter() {} + for epoch in self.epoch_iter() { + for cell9 in self.cell9_iter().filter(|cell| cell.epoch == epoch) { + if factor 1.0 { + } else { + } + } + } Ok(()) } @@ -771,9 +779,11 @@ impl IONEX { Box::new(self.record.map.keys().map(|k| k.epoch).unique().sorted()) } - /// Iterates this [IONEX] by a group of 9 neighboring [MapCell]s, + /// Iterates this [IONEX] by a group of 9 synchronous neighboring [MapCell]s, /// which is particularly convenient for accurate interpolation and upscaling. - pub fn iter_cell9(&self) -> Box + '_> {} + pub fn cell9_iter(&self) -> Box + '_> { + + } /// Modify the grid spacing (quantization) while preserving the dimensions, /// and interpolates the TEC values. From c759741c3fea6e329dfb1e97dbfdb197d084178e Mon Sep 17 00:00:00 2001 From: "Guillaume W. Bres" Date: Tue, 11 Nov 2025 11:23:49 +0100 Subject: [PATCH 23/25] up Signed-off-by: Guillaume W. Bres --- src/{cell.rs => cell/mod.rs} | 145 +++++---------- src/cell/three_by_three.rs | 284 ++++++++++++++++++++++++++++++ src/cell9.rs | 88 ---------- src/lib.rs | 332 ++++++++++++++++++----------------- 4 files changed, 507 insertions(+), 342 deletions(-) rename src/{cell.rs => cell/mod.rs} (91%) create mode 100644 src/cell/three_by_three.rs delete mode 100644 src/cell9.rs diff --git a/src/cell.rs b/src/cell/mod.rs similarity index 91% rename from src/cell.rs rename to src/cell/mod.rs index 70b05b8..4f3c232 100644 --- a/src/cell.rs +++ b/src/cell/mod.rs @@ -7,33 +7,8 @@ use crate::{ rectangle_to_cardinals, }; -#[derive(Debug, Copy, Clone, PartialEq, Eq, Hash, PartialOrd)] -#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))] -pub enum Cardinal { - /// NE [Cardinal] - NorthEast, - - /// N [Cardinal] - North, - - /// NW [Cardinal] - NorthWest, - - /// W [Cardinal] - West, - - /// SW [Cardinal] - SouthWest, - - /// S [Cardinal] - South, - - /// SE [Cardinal] - SouthEast, - - /// E [Cardinal] - East, -} +mod three_by_three; +pub use three_by_three::Cell3x3; #[derive(Debug, Default, Clone, Copy, PartialEq)] pub struct TecPoint { @@ -171,6 +146,11 @@ impl MapCell { self.epoch == rhs.epoch } + /// Returns true if both [MapCell]s describe the same spatial region at the same instant. + pub fn spatial_temporal_match(&self, rhs: &Self) -> bool { + self.spatial_match(rhs) && self.temporal_match(rhs) + } + /// Returns true if self is the northern neighbor of provided (rhs) [MapCell]. pub fn is_northern_neighbor(&self, rhs: &Self) -> bool { rhs.north_east.point == self.south_east.point @@ -355,6 +335,12 @@ impl MapCell { self } + /// Copies and updates the temporal instant + pub fn with_epoch(mut self, epoch: Epoch) -> Self { + self.epoch = epoch; + self + } + /// Returns the (latitude, longitude) span of this [MapCell] /// as tuplet in degrees pub fn latitude_longitude_span_degrees(&self) -> (f64, f64) { @@ -448,17 +434,17 @@ impl MapCell { Ok(TEC::from_tecu(tecu)) } - /// Returns a downscaled (resized minized) new ROI applying the interpolation equation - /// on each corners. - pub fn downscaled(&self, factor: f64) -> Result { + /// Returns a stretched (either upscaled or downscaled, resized in dimension) ROI, + /// by applying the interpolation equation on each corners. + /// Although this operation may apply to any [MapCell], for best precision + /// it is recommend to restrict the upscaling factor to small values (<2), otherwise, + /// the [Cell3x3] should be prefered and will give more accurate results, by + /// taking into account the neighboring cells. + pub fn stretching_mut(&mut self, factor: f64) -> Result<(), Error> { if !factor.is_normal() { return Err(Error::InvalidStretchFactor); } - if factor > 1.0 { - return Err(Error::InvalidStretchFactor); - } - // apply interpolation eq. at 4 coordinates let (north_east, north_west, south_east, south_west) = ( Point::new( @@ -479,78 +465,43 @@ impl MapCell { ), ); - Ok(MapCell { - north_east: TecPoint { - tec: self.spatial_tec_interp(north_east), + let (north_east, north_west, south_east, south_west) = ( + TecPoint { point: north_east, + tec: self.spatial_tec_interp(north_east)?, }, - north_west: TecPoint { - tec: self.spatial_tec_interp(north_west), + TecPoint { point: north_west, + tec: self.spatial_tec_interp(north_west)?, }, - south_west: TecPoint { - tec: self.spatial_tec_interp(south_west), - point: south_west, - }, - south_east: TecPoint { - tec: self.spatial_tec_interp(south_east), + TecPoint { point: south_east, + tec: self.spatial_tec_interp(south_east)?, }, - }) - } - - /// Returns a upscaled (resized upscaled) new ROI applying the interpolation equation - /// on each corners. This method should be used for very small stretch ratios. - /// For higher accuracy, you should work with a [Cell9] grouping and [Cell9::upscale] - /// in particular, which takes advantage of neighboring information on the map, to minimize - /// the error on the resulting ROI. - pub fn upscaled(&self, factor: f64) -> Result { - if !factor.is_normal() { - return Err(Error::InvalidStretchFactor); - } + TecPoint { + point: south_west, + tec: self.spatial_tec_interp(south_west)?, + }, + ); - if factor < 1.0 { - return Err(Error::InvalidStretchFactor); - } + self.north_east = north_east; + self.north_west = north_west; + self.south_east = south_east; + self.south_west = south_west; - // apply interpolation eq. at 4 coordinates - let (north_east, north_west, south_east, south_west) = ( - Point::new( - self.north_east.point.x() * factor, - self.north_east.point.y() * factor, - ), - Point::new( - self.north_west.point.x() * factor, - self.north_west.point.y() * factor, - ), - Point::new( - self.south_east.point.x() * factor, - self.south_east.point.y() * factor, - ), - Point::new( - self.south_west.point.x() * factor, - self.south_west.point.y() * factor, - ), - ); + Ok(()) + } - Ok(MapCell { - north_east: TecPoint { - tec: self.spatial_tec_interp(north_east), - point: north_east, - }, - north_west: TecPoint { - tec: self.spatial_tec_interp(north_west), - point: north_west, - }, - south_west: TecPoint { - tec: self.spatial_tec_interp(south_west), - point: south_west, - }, - south_east: TecPoint { - tec: self.spatial_tec_interp(south_east), - point: south_east, - }, - }) + /// Returns a stretched (either upscaled or downscaled, resized in dimension) ROI, + /// by applying the interpolation equation on each corners. + /// Although this operation may apply to any [MapCell], for best precision + /// it is recommend to restrict the upscaling factor to small values (<2), otherwise, + /// the [Cell3x3] should be prefered and will give more accurate results, by + /// taking into account the neighboring cells. + pub fn stretched(&self, factor: f64) -> Result { + let mut s = self.clone(); + s.stretching_mut(factor)?; + Ok(s) } // /// Determines the northeastern cell amongst a grouping of 4 diff --git a/src/cell/three_by_three.rs b/src/cell/three_by_three.rs new file mode 100644 index 0000000..7d4a400 --- /dev/null +++ b/src/cell/three_by_three.rs @@ -0,0 +1,284 @@ +use std::collections::HashMap; + +use geo::{Contains, GeodesicArea, Geometry, Point, Rect}; + +use crate::prelude::{Epoch, Error, MapCell, TEC}; + +// #[derive(Debug, Copy, Clone, PartialEq, Eq, Hash, PartialOrd)] +// #[cfg_attr(feature = "serde", derive(Serialize, Deserialize))] +// pub enum Cardinal { +// /// NE [Cardinal] +// NorthEast, +// +// /// N [Cardinal] +// North, +// +// /// NW [Cardinal] +// NorthWest, +// +// /// W [Cardinal] +// West, +// +// /// SW [Cardinal] +// SouthWest, +// +// /// S [Cardinal] +// South, +// +// /// SE [Cardinal] +// SouthEast, +// +// /// E [Cardinal] +// East, +// } + +/// A synchronous 3x3 ROI made of a central [MapCell] element and 8 neighboring [MapCell]s. +#[derive(Debug, Default, Clone, PartialEq)] +pub struct Cell3x3 { + /// The central [MapCell] + pub center: MapCell, + + /// The northeastern neighboring [MapCell] + pub northeast: MapCell, + + /// The northern neighboring [MapCell] + pub north: MapCell, + + /// The northwestern neighboring [MapCell] + pub northwest: MapCell, + + /// The western neighboring [MapCell] + pub west: MapCell, + + /// The southwestern neighboring [MapCell] + pub southwest: MapCell, + + /// The southern neighboring [MapCell] + pub south: MapCell, + + /// The southeastern neighboring [MapCell] + pub southeast: MapCell, + + /// The eastern neighboring [MapCell] + pub east: MapCell, +} + +impl Cell3x3 { + /// Returns true if both [Cell3x3] cells represent the same spatial region + pub fn spatial_match(&self, rhs: &Self) -> bool { + self.center.spatial_match(&rhs.center) + && self.northeast.spatial_match(&rhs.northeast) + && self.north.spatial_match(&rhs.north) + && self.northwest.spatial_match(&rhs.northwest) + && self.west.spatial_match(&rhs.west) + && self.southwest.spatial_match(&rhs.southwest) + && self.south.spatial_match(&rhs.south) + && self.southeast.spatial_match(&rhs.southeast) + && self.east.spatial_match(&rhs.east) + } + + /// Returns true if both [Cell3x3] cells are synchronous. + pub fn temporal_match(&self, rhs: &Self) -> bool { + self.center.temporal_match(&rhs.center) + && self.northeast.temporal_match(&rhs.northeast) + && self.north.temporal_match(&rhs.north) + && self.northwest.temporal_match(&rhs.northwest) + && self.west.temporal_match(&rhs.west) + && self.southwest.temporal_match(&rhs.southwest) + && self.south.temporal_match(&rhs.south) + && self.southeast.temporal_match(&rhs.southeast) + && self.east.temporal_match(&rhs.east) + } + + /// Returns true if both [Cell3x3] cells represent the same spatial region + /// at the same instant. + pub fn spatial_temporal_match(&self, rhs: &Self) -> bool { + self.spatial_match(rhs) && self.temporal_match(rhs) + } + + /// Builds a new [Cell3x3] updated in time + pub fn with_epoch(mut self, epoch: Epoch) -> Self { + self.center.epoch = epoch; + self.northeast.epoch = epoch; + self.north.epoch = epoch; + self.northwest.epoch = epoch; + self.east.epoch = epoch; + self.southeast.epoch = epoch; + self.south.epoch = epoch; + self.southwest.epoch = epoch; + self.west.epoch = epoch; + self + } + + /// Builds a new [Cell3x3] with updated central cell, which must be synchronous. + pub fn with_central_cell(mut self, cell: MapCell) -> Result { + if cell.epoch == self.center.epoch { + self.center = cell; + Ok(self) + } else { + Err(Error::TemporalMismatch) + } + } + + /// Builds a new [Cell3x3] with updated western cell, which must be synchronous. + pub fn with_western_cell(mut self, cell: MapCell) -> Result { + if cell.epoch == self.west.epoch { + self.west = cell; + Ok(self) + } else { + Err(Error::TemporalMismatch) + } + } + + /// Builds a new [Cell3x3] with updated eastern cell, which must be synchronous. + pub fn with_eastern_cell(mut self, cell: MapCell) -> Result { + if cell.epoch == self.east.epoch { + self.east = cell; + Ok(self) + } else { + Err(Error::TemporalMismatch) + } + } + + /// Builds a new [Cell3x3] with updated southern cell, which must be synchronous. + pub fn with_southern_cell(mut self, cell: MapCell) -> Result { + if cell.epoch == self.south.epoch { + self.south = cell; + Ok(self) + } else { + Err(Error::TemporalMismatch) + } + } + + /// Builds a new [Cell3x3] with updated northern cell, which must be synchronous. + pub fn with_northern_cell(mut self, cell: MapCell) -> Result { + if cell.epoch == self.north.epoch { + self.north = cell; + Ok(self) + } else { + Err(Error::TemporalMismatch) + } + } + + /// Builds a new [Cell3x3] with updated southwestern cell, which must be synchronous. + pub fn with_southwestern_cell(mut self, cell: MapCell) -> Result { + if cell.epoch == self.southwest.epoch { + self.southwest = cell; + Ok(self) + } else { + Err(Error::TemporalMismatch) + } + } + + /// Builds a new [Cell3x3] with updated southeastern cell, which must be synchronous. + pub fn with_southeastern_cell(mut self, cell: MapCell) -> Result { + if cell.epoch == self.southeast.epoch { + self.southeast = cell; + Ok(self) + } else { + Err(Error::TemporalMismatch) + } + } + + /// Builds a new [Cell3x3] with updated northwestern cell, which must be synchronous. + pub fn with_northwestern_cell(mut self, cell: MapCell) -> Result { + if cell.epoch == self.northwest.epoch { + self.northwest = cell; + Ok(self) + } else { + Err(Error::TemporalMismatch) + } + } + + /// Builds a new [Cell3x3] with updated northeastern cell, which must be synchronous. + pub fn with_northeastern_cell(mut self, cell: MapCell) -> Result { + if cell.epoch == self.northeast.epoch { + self.northeast = cell; + Ok(self) + } else { + Err(Error::TemporalMismatch) + } + } + + /// Builds a new [Cell3x3] from a slice of 9 unordered [MapCell]s, by electing a central element (if feasible) + /// and the neighboring cells, which must all be synchronous. + pub fn from_slice(cells: [MapCell; 9]) -> Option { + for i in 0..9 { + // determine whether the ith cell is the potential center + + // must be neighbor with all other 8 ROIs and must be synchronous + let mut all_synchronous_neighbors = true; + + for j in 0..9 { + if j != i { + if !cells[i].is_neighbor(&cells[j]) { + all_synchronous_neighbors = false; + break; + } + if !cells[i].temporal_match(&cells[j]) { + all_synchronous_neighbors = false; + break; + } + } + } + + if all_synchronous_neighbors { + let mut count = 0; + + let mut ret = Self::default() + .with_epoch(cells[i].epoch) + .with_central_cell(cells[i]) + .unwrap(); // infaillible + + // i is the central ROI, order other ROIs correctly & exit + for j in 0..9 { + if j != i { + if cells[j].is_northwestern_neighbor(&cells[i]) { + count += 1; + ret = ret.with_northwestern_cell(cells[j]).unwrap(); + // infaillible + } else if cells[j].is_northeastern_neighbor(&cells[i]) { + count += 1; + ret = ret.with_northeastern_cell(cells[j]).unwrap(); + // infaillible + } else if cells[j].is_northern_neighbor(&cells[i]) { + count += 1; + ret = ret.with_northern_cell(cells[j]).unwrap(); // infaillible + } else if cells[j].is_southwestern_neighbor(&cells[i]) { + count += 1; + ret = ret.with_southwestern_cell(cells[j]).unwrap(); + // infaillible + } else if cells[j].is_southeastern_neighbor(&cells[i]) { + count += 1; + ret = ret.with_southeastern_cell(cells[j]).unwrap(); + // infaillible + } else if cells[j].is_southern_neighbor(&cells[i]) { + count += 1; + ret = ret.with_southern_cell(cells[j]).unwrap(); // infaillible + } else if cells[j].is_eastern_neighbor(&cells[i]) { + count += 1; + ret = ret.with_eastern_cell(cells[j]).unwrap(); // infaillible + } else if cells[j].is_western_neighbor(&cells[i]) { + count += 1; + ret = ret.with_western_cell(cells[j]).unwrap(); // infaillible + } + } + } + + if count == 8 { + // 3x3 roi completed + return Some(ret); + } + } + } + + None + } + + /// Returns a stretched (spatially upscaled or downscaled) [MapCell] by + /// stretching the central element and taking the relative neighboring values into + /// account. + pub fn stretched(&self, factor: f64) -> Result { + Err(Error::OutsideSpatialBoundaries) // TODO + } +} diff --git a/src/cell9.rs b/src/cell9.rs deleted file mode 100644 index 43e43eb..0000000 --- a/src/cell9.rs +++ /dev/null @@ -1,88 +0,0 @@ -use std::collections::HashMap; - -use geo::{Contains, GeodesicArea, Geometry, Point, Rect}; - -use crate::prelude::{Cardinal, Epoch, Error, MapCell, TEC}; - -/// A structure holding 9 synchronous neighboring [MapCells] -#[derive(Debug, Default, Clone, PartialEq)] -pub struct Cell9 { - /// The central [MapCell] - pub central: MapCell, - - /// 8 Neighboring [MapCells] - pub neighbors: HashMap, -} - -impl Cell9 { - /// Returns true if this [Cell9] definition is "complete" - pub(crate) fn is_complete(&self) -> bool { - self.neighbors.len() == 8 - } - - /// Builds a new updated [Cell9] with updated central cell. - pub fn with_central_cell(mut self, center: MapCell) -> Self { - self.central = center; - self - } - - /// Builds a new [Cell9] conveniently from 9 [MapCells], correctly - /// electing the central ROI, and the 8 neighbors. - /// Returns None if no central ROI exists or could not determine 8 neighbors. - pub fn from_slice(cells: [MapCell; 9]) -> Option { - for i in 0..9 { - // determine whether the ith cell is the potential center - let mut all_neighbors = true; - // must be neighbor with all other 8 ROIs - for j in 0..9 { - if j != i { - if !cells[i].is_neighbor(&cells[j]) { - all_neighbors = false; - break; - } - } - } - - if all_neighbors { - let mut ret = Self::default().with_central_cell(cells[i]); - - // i is the central ROI, order other ROIs correctly & exit - for j in 0..9 { - if j != i { - if cells[j].is_northwestern_neighbor(&cells[i]) { - ret.neighbors.insert(Cardinal::NorthWest, cells[j]); - } else if cells[j].is_northeastern_neighbor(&cells[i]) { - ret.neighbors.insert(Cardinal::NorthEast, cells[j]); - } else if cells[j].is_northern_neighbor(&cells[i]) { - ret.neighbors.insert(Cardinal::North, cells[j]); - } else if cells[j].is_southwestern_neighbor(&cells[i]) { - ret.neighbors.insert(Cardinal::SouthWest, cells[j]); - } else if cells[j].is_southeastern_neighbor(&cells[i]) { - ret.neighbors.insert(Cardinal::SouthEast, cells[j]); - } else if cells[j].is_southern_neighbor(&cells[i]) { - ret.neighbors.insert(Cardinal::South, cells[j]); - } - } - } - - if ret.is_complete() { - return Some(ret); - } - } - } - - None - } - - /// Upscale or downscale the central ROI of this [Cell9] by a fractional positive scalar. - pub fn upscale(&self, factor: f64) -> Result { - // interpolate the 4 new central coordinates - // interpolate using the 4 neighbors - // take the average of both for each 4 coordinntes - } - - pub fn downscale(&self, factor: f64) -> Result { - // interpolate the central ROI - self.central.downscale(factor) - } -} diff --git a/src/lib.rs b/src/lib.rs index b8d4ba4..62abf2e 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -33,7 +33,6 @@ pub mod tec; pub mod version; mod cell; -mod cell9; mod coordinates; mod epoch; mod ionosphere; @@ -61,8 +60,7 @@ use flate2::{read::GzDecoder, write::GzEncoder, Compression as GzCompression}; use hifitime::prelude::Epoch; use crate::{ - cell::Cell9, - cell::{Cardinal, MapCell, TecPoint}, + cell::{Cell3x3, MapCell, TecPoint}, coordinates::QuantizedCoordinates, error::{Error, FormattingError, ParsingError}, file_attributes::{FileAttributes, Region}, @@ -78,7 +76,7 @@ pub mod prelude { // export pub use crate::{ bias::BiasSource, - cell::{Cardinal, Cell9, MapCell}, + cell::{Cell3x3, MapCell}, error::{Error, FormattingError, ParsingError}, file_attributes::*, grid::{Axis, Grid}, @@ -687,23 +685,24 @@ impl IONEX { Ok(s) } - /// Modify the grid spacing (quantization) while preserving the dimensions, - /// and interpolates the new TEC values. - /// - /// This only modify the grid quantization, the map dimensions are preserved. - /// - /// ## Input - /// - axis: [Axis] to be resampled - /// - factor: - /// - > 1.0: spatial upscaling. The grid precision increases, - /// we interpolate the TEC at new intermiediate coordinates. - /// - < 1.0: downscaling. The grid precision decreases. - /// - 0.0: invalid - pub fn spatially_resampled(&self, axis: Axis, factor: f64) -> Result { - let mut s = self.clone(); - s.spatial_resampling_mut(axis, factor)?; - Ok(s) - } + // TODO + // /// Modify the grid spacing (quantization) while preserving the dimensions, + // /// and interpolates the new TEC values. + // /// + // /// This only modify the grid quantization, the map dimensions are preserved. + // /// + // /// ## Input + // /// - axis: [Axis] to be resampled + // /// - factor: + // /// - > 1.0: spatial upscaling. The grid precision increases, + // /// we interpolate the TEC at new intermiediate coordinates. + // /// - < 1.0: downscaling. The grid precision decreases. + // /// - 0.0: invalid + // pub fn spatially_resampled(&self, axis: Axis, factor: f64) -> Result { + // let mut s = self.clone(); + // s.spatial_resampling_mut(axis, factor)?; + // Ok(s) + // } /// Modify the grid dimensions by a positive, possibly fractional number, /// and interpolates the TEC values. @@ -764,10 +763,13 @@ impl IONEX { // synchronous spatial interpolation for epoch in self.epoch_iter() { - for cell9 in self.cell9_iter().filter(|cell| cell.epoch == epoch) { - if factor 1.0 { - } else { - } + for cell3x3 in self + .cell3x3_iter() + .filter(|cell| cell.center.epoch == epoch) + { + // if factor 1.0 { + // } else { + // } } } @@ -779,64 +781,63 @@ impl IONEX { Box::new(self.record.map.keys().map(|k| k.epoch).unique().sorted()) } - /// Iterates this [IONEX] by a group of 9 synchronous neighboring [MapCell]s, + /// Iterates this [IONEX] by a 3x3 group of neighboring [MapCell]s expressed as [Cell3x3], /// which is particularly convenient for accurate interpolation and upscaling. - pub fn cell9_iter(&self) -> Box + '_> { - - } + pub fn cell3x3_iter(&self) -> Box + '_> {} - /// Modify the grid spacing (quantization) while preserving the dimensions, - /// and interpolates the TEC values. - /// - /// This only modify the grid quantization, the map dimensions are preserved. - /// - /// ## Input - /// - axis: [Axis] to be resampled - /// - factor: - /// - > 1.0: upscaling - /// - < 1.0: downscaling - /// - 0.0: invalid - pub fn spatial_resampling_mut(&mut self, axis: Axis, factor: f64) -> Result<(), Error> { - // Stretch header specs - match axis { - Axis::Latitude => { - self.header.grid.latitude.resample_mut(factor)?; - }, - Axis::Longitude => { - self.header.grid.longitude.resample_mut(factor)?; - }, - Axis::Altitude => { - self.header.grid.altitude.resample_mut(factor)?; - }, - Axis::Planar => { - self.header.grid.latitude.resample_mut(factor)?; - self.header.grid.longitude.resample_mut(factor)?; - }, - Axis::All => { - self.header.grid.latitude.resample_mut(factor)?; - self.header.grid.longitude.resample_mut(factor)?; - self.header.grid.altitude.resample_mut(factor)?; - }, - } + // TODO + // /// Modify the grid spacing (quantization) while preserving the dimensions, + // /// and interpolates the TEC values. + // /// + // /// This only modify the grid quantization, the map dimensions are preserved. + // /// + // /// ## Input + // /// - axis: [Axis] to be resampled + // /// - factor: + // /// - > 1.0: upscaling + // /// - < 1.0: downscaling + // /// - 0.0: invalid + // pub fn spatial_resampling_mut(&mut self, axis: Axis, factor: f64) -> Result<(), Error> { + // // Stretch header specs + // match axis { + // Axis::Latitude => { + // self.header.grid.latitude.resample_mut(factor)?; + // }, + // Axis::Longitude => { + // self.header.grid.longitude.resample_mut(factor)?; + // }, + // Axis::Altitude => { + // self.header.grid.altitude.resample_mut(factor)?; + // }, + // Axis::Planar => { + // self.header.grid.latitude.resample_mut(factor)?; + // self.header.grid.longitude.resample_mut(factor)?; + // }, + // Axis::All => { + // self.header.grid.latitude.resample_mut(factor)?; + // self.header.grid.longitude.resample_mut(factor)?; + // self.header.grid.altitude.resample_mut(factor)?; + // }, + // } - let lat_pairs = self.header.grid.latitude.quantize().tuple_windows(); - let long_pairs = self.header.grid.longitude.quantize().tuple_windows(); + // let lat_pairs = self.header.grid.latitude.quantize().tuple_windows(); + // let long_pairs = self.header.grid.longitude.quantize().tuple_windows(); - // applies the stretching + // // applies the stretching - // synchronous spatial interpolation - for epoch in self.epoch_iter() { - if factor > 1.0 { - // upscaling - } else { - // downscaling - } - } + // // synchronous spatial interpolation + // for epoch in self.epoch_iter() { + // if factor > 1.0 { + // // upscaling + // } else { + // // downscaling + // } + // } - // saturate the possible overflows to the world map proj + // // saturate the possible overflows to the world map proj - Ok(()) - } + // Ok(()) + // } /// Upscale (upsample) or Downscale (downsample) this mutable [IONEX], /// modifying the stretch on the temporal axis. @@ -973,25 +974,53 @@ impl IONEX { ) } - /// Obtain a [MapCell] at specified point in time, wrapping this ROI. - /// The temporal instant must be within the current temporal axis. - /// The returned [MapCell] is a bounding rectangle, of 4 corners - /// correctly interpolated in space and time, than you can further manipulate. - /// Depending on the ROI dimensions, this is either: - /// - a unitary [MapCell]: the smallest region we can represent in this [IONEX]. - /// Prefer [Self::unitary_roi_at] in this case, to not bother designing the unit [Rect]. - /// - or the closest fitting [MapCell] (wrapping boundaries) that we could fit - pub fn roi_at(&self, epoch: Epoch, roi: Rect) -> Result { + /// Form a possibly interpolated [MapCell] at specified point in time, wrapping this ROI. + /// The ROI is described by a [Geometry], which can be complex. + /// Unlike [Self::unitary_roi_at], the returned [MapCell] may not be a map quantization cell, + /// but might result of the combination of several, in case the ROI does not fit in a single quantization step. + /// + /// ## Input + /// - epoch: the returned cell will be synchronous to this [Epoch]. + /// If this instant lines up with the temporal axis, the process will not involve temporal interpolation. + /// When the instant does not line up with the temporal axis, we will use temporal interpolation, + /// to deduce the most precise and correct values. + /// In any case, this instant must fit within the temporal axis, after the first [Epoch] + /// and before the last [Epoch] described in [Header]. + /// + /// - roi: [Geometry] defining the local region we want to wrap (fully contained by returned cell). + pub fn roi_at(&self, epoch: Epoch, roi: Geometry) -> Result { // determine whether this is within the temporal axis if epoch < self.header.epoch_of_first_map || epoch > self.header.epoch_of_last_map { return Err(Error::OutsideTemporalBoundaries); } - // determine this is within the map - if !self.bounding_rect_degrees().contains(&Geometry::Rect(roi)) { + // convert the ROI to its bounding rectangle + let roi = roi.bounding_rect().ok_or(Error::UndefinedBoundaries)?; + + // apply the modulo operation in case provided coordinates are "incorrect". + // in case of worldwide ionex, this makes this operation work at all times, conveniently. + let roi = Rect::new( + coord!(x: roi.min().x, y: roi.min().y), + coord!(x: roi.max().x, y: roi.max().y), + ); + + // check the ROI is within the map. + // in case of regional IONEX, we might not be able to return + if !self.bounding_rect_degrees().contains(&roi) { return Err(Error::OutsideSpatialBoundaries); } + // determine whether this lies within a single element or not + let mut unitary_element = true; + + if roi.width() > self.header.grid.longitude.spacing { + unitary_element = false; + } + + if roi.height() > self.header.grid.latitude.spacing { + unitary_element = false; + } + // determine whether we need temporal interpolation or not let mut needs_temporal_interp = true; let mut t = self.header.epoch_of_first_map; @@ -1005,98 +1034,87 @@ impl IONEX { t += self.header.sampling_period; } - // we don't need spatial interpolation if the ROI height and width - // are a perfect multiple of the grid dimensions + // determine whether we need spatial interpolation or not. + // we don't need spatial interpolation when the ROI height and width + // perfectly line up with the IONEX grid. let (width, height) = (roi.width(), roi.height()); - let mut needs_lat_interpolation = - height.rem_euclid(self.header.grid.latitude.spacing) == 0.0; - let mut needs_long_interpolation = - width.rem_euclid(self.header.grid.longitude.spacing) == 0.0; + let (needs_lat_interp, needs_long_interp) = ( + height.rem_euclid(self.header.grid.latitude.spacing) == 0.0, + width.rem_euclid(self.header.grid.longitude.spacing) == 0.0, + ); - let needs_spatial_interpolation = needs_lat_interpolation || needs_long_interpolation; - if needs_spatial_interpolation { - panic!("not supported yet"); + if unitary_element { + // this can be reduced to a unitary MapCell + } else { + if needs_temporal_interp { + if needs_lat_interp || needs_long_interp { + panic!("not supported yet"); + } else { + } + } else { + if needs_lat_interp || needs_long_interp { + panic!("not supported yet"); + } else { + } + } } // TODO Err(Error::TemporalMismatch) } - /// Obtain the smallest [MapCell] (map quantization) at provided point in time - /// and coordinates. + /// Obtain the [MapCell] (smallest map ROI) at provided point in time and containing provided coordinates. + /// We will select the synchronous [MapCell] that contains the given coordinates. /// /// ## Input /// - epoch: [Epoch] that must fit within the temporal axis. /// If this instant aligns with the sampling rate, this process will not require temporal interpolation. - /// - coordinates: 2D [Point] that must lie within the map. + /// If the instant does not line up with the sampling rate, we will return a temporally interpolated ROI + /// at the specified point in time. + /// - contains: [Geometry] to be fully contained by the returned /// If the coordinates align with the grid, this process will not require spatial interpolation. - /// - card: [Cardinal] defining which corner these coordinates will represent in the returned [MapCell]. - /// For example: - /// - [Cardinal::NorthEast] means the northeastern component of fitted ROI [MapCell] will lie at these coordinates. - pub fn unitary_roi_at( - &self, - epoch: Epoch, - coordinates: Point, - card: Cardinal, - ) -> Result { - // builds the unitary ROI (as Rect) matching those coordinates and cardinal - let (lat_pairs, long_pairs) = ( - self.header.grid.latitude.quantize().tuple_windows(), - self.header.grid.longitude.quantize().tuple_windows(), - ); - - let (lat_spacing, long_spacing) = ( - self.header.grid.latitude.spacing, - self.header.grid.longitude.spacing, - ); - - for (lat1, lat2) in lat_pairs { - for (long1, long2) in long_pairs.clone() { - let (lat1, lat2, long1, long2) = ( - lat1.real_value(), - lat2.real_value(), - long1.real_value(), - long2.real_value(), - ); - match card { - Cardinal::NorthEast => { - if lat1 == coordinates.y() && long1 == coordinates.x() { - let roi = - Rect::new(coord!(x: long1, y: lat1), coord!(x: long2, y: lat2)); - - return self.roi_at(epoch, roi); - } - }, - Cardinal::SouthEast => { - if lat1 == coordinates.y() && long1 == coordinates.x() { - let roi = - Rect::new(coord!(x: long1, y: lat1), coord!(x: long2, y: lat2)); + pub fn unitary_roi_at(&self, epoch: Epoch, coordinates: Point) -> Option { + // determine whether we need temporal interpolation or not + let mut needs_temporal_interp = true; + let mut t = self.header.epoch_of_first_map; - return self.roi_at(epoch, roi); - } - }, - Cardinal::SouthWest => { - if lat1 == coordinates.y() && long1 == coordinates.x() { - let roi = - Rect::new(coord!(x: long1, y: lat1), coord!(x: long2, y: lat2)); + while t < self.header.epoch_of_last_map { + if t == epoch { + needs_temporal_interp = false; + break; + } - return self.roi_at(epoch, roi); - } - }, - Cardinal::NorthWest => { - if lat1 == coordinates.y() && long1 == coordinates.x() { - let roi = - Rect::new(coord!(x: long1, y: lat1), coord!(x: long2, y: lat2)); + t += self.header.sampling_period; + } - return self.roi_at(epoch, roi); + let coordinates = Geometry::Point(coordinates); + + if needs_temporal_interp { + for (t0, t1) in self.epoch_iter().tuple_windows() { + if t0 < epoch && t1 > epoch { + for (cell0, cell1) in self + .synchronous_map_cell_iter(t0) + .zip(self.synchronous_map_cell_iter(t1)) + { + if cell0.contains(&coordinates) && cell1.contains(&coordinates) { + // TODO + // if let Ok(interpolated) = cell0.temporally_interpolated(epoch, &cell1)? { + // return Some(interpolated); + // } } - }, + } + } + } + } else { + for cell in self.synchronous_map_cell_iter(epoch) { + if cell.contains(&coordinates) { + return Some(cell); } } } - Err(Error::OutsideSpatialBoundaries) + None } /// Obtain the best suited [MapCell] spatially wrapping this Geometry that contains following [Geometry]. From 8e23b24074ebdef1df31c453fe424bc1dd05fae4 Mon Sep 17 00:00:00 2001 From: "Guillaume W. Bres" Date: Tue, 11 Nov 2025 13:17:01 +0100 Subject: [PATCH 24/25] up Signed-off-by: Guillaume W. Bres --- src/lib.rs | 353 +++++++++++++++++++++++----------------------- src/record/mod.rs | 101 ++++++++----- 2 files changed, 244 insertions(+), 210 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 62abf2e..aa1b55f 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -57,7 +57,7 @@ use geo::{coord, BoundingRect, Contains, Geometry, LineString, Point, Polygon, R #[cfg(feature = "flate2")] use flate2::{read::GzDecoder, write::GzEncoder, Compression as GzCompression}; -use hifitime::prelude::Epoch; +use hifitime::prelude::{Epoch, TimeSeries}; use crate::{ cell::{Cell3x3, MapCell, TecPoint}, @@ -599,14 +599,19 @@ impl IONEX { /// spatial dimensions are modified by this operation. /// [Polygon::bounding_rect] must be defined for this operation to work correctly. pub fn to_regional_ionex(&self, roi: Polygon) -> Result { - let mut ionex = IONEX::default(); + let mut ionex = IONEX::default().with_header(self.header.clone()); - let bounding_rect = roi.bounding_rect().ok_or(Error::UndefinedBoundaries)?; + // simply restrict to wrapping area + let roi = Geometry::Polygon(roi); - let (min_long, min_lat) = (bounding_rect.min().x, bounding_rect.min().y); - let (max_long, max_lat) = (bounding_rect.max().x, bounding_rect.max().y); + let cells = self + .map_cell_iter() + .filter(|cell| cell.contains(&roi)) + .collect::>(); + + ionex.record = Record::from_map_cells(&cells, self.header.grid.altitude.start); - // attributes + // update attributes match &self.attributes { Some(attributes) => { ionex.attributes = Some(attributes.clone().regionalized()); @@ -617,73 +622,45 @@ impl IONEX { }, } - // copy header - ionex.header = self.header.clone(); - - // upscale so the grid boundaries become a perfect multiple - let long_upscaling = 0.0; - let lat_upscaling = 0.0; - - // spatial upscaling - ionex.spatial_stretching_mut(Axis::Latitude, lat_upscaling)?; - ionex.spatial_stretching_mut(Axis::Longitude, long_upscaling)?; - - // simply restrict to ROI area - let geometry = Geometry::Polygon(roi); - - let cells = self - .map_cell_iter() - .filter(|cell| cell.contains(&geometry)) - .collect::>(); - - ionex.record = Record::from_map_cells( - self.header.grid.altitude.start, - min_lat, - max_lat, - min_long, - max_long, - &cells, - ); - Ok(ionex) } - /// Modify the grid dimensions by a positive, possibly fractional number, - /// and interpolates the TEC values. - /// - /// This does not modify the grid quantization, only the dimensions. - /// - /// ## Input - /// - axis: [Axis] to be stretched - /// - factor: - /// - > 1.0: upscaling - /// - < 1.0: downscaling - /// - 0.0: invalid - /// - /// Example 1: spatial upscaling - /// ``` - /// use ionex::prelude::{ - /// IONEX, - /// Axis, - /// }; - /// - /// // grab a local dataset, upscaling hardly applies to worldwide maps.. - /// // this example is a data/IONEX/V1/CKMG0020.22I.gz that went through a x0.5 downscale. - /// let ionex = IONEX::from_gzip_file("data/IONEX/V1/CKMR0020.22I.gz").unwrap_or_else(|e| { - /// panic!("Failed to parse CKMR0020: {}", e); - /// }); - /// - /// // spatial dimensions upscale including interpolation: upscales to worldwide dimensions - /// let upscaled = ionex.spatially_stretched(Axis::Planar, 2.0) - /// .unwrap(); - /// - /// // testbench: compares these results to the original data - /// ``` - pub fn spatially_stretched(&self, axis: Axis, factor: f64) -> Result { - let mut s = self.clone(); - s.spatial_stretching_mut(axis, factor)?; - Ok(s) - } + // /// Modify the grid dimensions by a positive, possibly fractional number, + // /// and interpolates the TEC values. + // /// + // /// This does not modify the grid quantization, only the dimensions. + // /// + // /// ## Input + // /// - axis: [Axis] to be stretched + // /// - factor: + // /// - > 1.0: upscaling + // /// - < 1.0: downscaling + // /// - 0.0: invalid + // /// + // /// Example 1: spatial upscaling + // /// ``` + // /// use ionex::prelude::{ + // /// IONEX, + // /// Axis, + // /// }; + // /// + // /// // grab a local dataset, upscaling hardly applies to worldwide maps.. + // /// // this example is a data/IONEX/V1/CKMG0020.22I.gz that went through a x0.5 downscale. + // /// let ionex = IONEX::from_gzip_file("data/IONEX/V1/CKMR0020.22I.gz").unwrap_or_else(|e| { + // /// panic!("Failed to parse CKMR0020: {}", e); + // /// }); + // /// + // /// // spatial dimensions upscale including interpolation: upscales to worldwide dimensions + // /// let upscaled = ionex.spatially_stretched(Axis::Planar, 2.0) + // /// .unwrap(); + // /// + // /// // testbench: compares these results to the original data + // /// ``` + // pub fn spatially_stretched(&self, axis: Axis, factor: f64) -> Result { + // let mut s = self.clone(); + // s.spatial_stretching_mut(axis, factor)?; + // Ok(s) + // } // TODO // /// Modify the grid spacing (quantization) while preserving the dimensions, @@ -704,87 +681,92 @@ impl IONEX { // Ok(s) // } - /// Modify the grid dimensions by a positive, possibly fractional number, - /// and interpolates the TEC values. - /// - /// This does not modify the grid quantization, only the dimensions. - /// - /// ## Input - /// - axis: [Axis] to be stretched - /// - factor: - /// - > 1.0: upscaling - /// - < 1.0: downscaling - /// - 0.0: invalid - pub fn spatial_stretching_mut(&mut self, axis: Axis, factor: f64) -> Result<(), Error> { - // the maximal factor supported for on iter is 2 - - // Update the header specs - match axis { - Axis::Latitude => { - self.header.grid.latitude.stretch_mut(factor)?; - }, - Axis::Longitude => { - self.header.grid.longitude.stretch_mut(factor)?; - }, - Axis::Altitude => { - self.header.grid.altitude.stretch_mut(factor)?; - }, - Axis::Planar => { - self.header.grid.latitude.stretch_mut(factor)?; - self.header.grid.longitude.stretch_mut(factor)?; - }, - Axis::All => { - self.header.grid.latitude.stretch_mut(factor)?; - self.header.grid.longitude.stretch_mut(factor)?; - self.header.grid.altitude.stretch_mut(factor)?; - }, - } + // /// Modify the grid dimensions by a positive, possibly fractional number, + // /// and interpolates the TEC values. + // /// + // /// This does not modify the grid quantization, only the dimensions. + // /// + // /// ## Input + // /// - axis: [Axis] to be stretched + // /// - factor: + // /// - > 1.0: upscaling + // /// - < 1.0: downscaling + // /// - 0.0: invalid + // pub fn spatial_stretching_mut(&mut self, axis: Axis, factor: f64) -> Result<(), Error> { + // if factor > 2.0 { + // // the maximal factor supported per iter is 2 + // // let mut fact = factor.max(2.0); + // // while fact > 1.0 { + // // self.spatial_stretching_mut(fact); + // // fact = factor /2.0; + // // } + // Err(Error::TemporalMismatch) // TODO + // } else { + + // // Update header specs + // match axis { + // Axis::Latitude => { + // self.header.grid.latitude.stretch_mut(factor)?; + // }, + // Axis::Longitude => { + // self.header.grid.longitude.stretch_mut(factor)?; + // }, + // Axis::Altitude => { + // self.header.grid.altitude.stretch_mut(factor)?; + // }, + // Axis::Planar => { + // self.header.grid.latitude.stretch_mut(factor)?; + // self.header.grid.longitude.stretch_mut(factor)?; + // }, + // Axis::All => { + // self.header.grid.latitude.stretch_mut(factor)?; + // self.header.grid.longitude.stretch_mut(factor)?; + // self.header.grid.altitude.stretch_mut(factor)?; + // }, + // } - // update the attributes - match self.attributes { - Some(ref mut attributes) => { - if self.header.grid.is_worldwide() { - attributes.region = Region::Worldwide; - } else { - attributes.region = Region::Regional; - } - }, - None => { - let mut attributes = FileAttributes::default(); - if self.header.grid.is_worldwide() { - attributes.region = Region::Worldwide; - } else { - attributes.region = Region::Regional; - } + // // update the attributes + // match self.attributes { + // Some(ref mut attributes) => { + // if self.header.grid.is_worldwide() { + // attributes.region = Region::Worldwide; + // } else { + // attributes.region = Region::Regional; + // } + // }, + // None => { + // let mut attributes = FileAttributes::default(); + // if self.header.grid.is_worldwide() { + // attributes.region = Region::Worldwide; + // } else { + // attributes.region = Region::Regional; + // } - self.attributes = Some(attributes); - }, - } + // self.attributes = Some(attributes); + // }, + // } - // synchronous spatial interpolation - for epoch in self.epoch_iter() { - for cell3x3 in self - .cell3x3_iter() - .filter(|cell| cell.center.epoch == epoch) - { - // if factor 1.0 { - // } else { - // } - } - } + // // synchronous spatial interpolation + // for epoch in self.epoch_iter() { + // for cell3x3 in self + // .map_cell3x3_iter() + // .filter(|cell| cell.center.epoch == epoch) + // { + // // if factor 1.0 { + // // } else { + // // } + // } + // } - Ok(()) - } + // Ok(()) + // } + // } /// Returns an iterator over [Epoch]s in chronological order. pub fn epoch_iter(&self) -> Box + '_> { Box::new(self.record.map.keys().map(|k| k.epoch).unique().sorted()) } - /// Iterates this [IONEX] by a 3x3 group of neighboring [MapCell]s expressed as [Cell3x3], - /// which is particularly convenient for accurate interpolation and upscaling. - pub fn cell3x3_iter(&self) -> Box + '_> {} - // TODO // /// Modify the grid spacing (quantization) while preserving the dimensions, // /// and interpolates the TEC values. @@ -853,50 +835,60 @@ impl IONEX { return Err(Error::InvalidStretchFactor); } + let new_dt = self.header.sampling_period * factor; + + if factor > 1.0 { + // temporal upscaling + } else { + // temporal decimation + } + // update header - self.header.sampling_period = self.header.sampling_period / factor; + self.header.sampling_period = new_dt; - // TODO: (t1, t2) iter Ok(()) } - /// Stretchs this mutable [IONEX] both in temporal and spatial dimensions, - /// modifying both dimensions. - /// When factor > 1.0 this is an upscaling operation, when factor < 1.0, - /// this is a downscaling operation. - /// This uses both spatial and temporal interpolation to form valid data points. - pub fn temporal_spatial_stretching_mut( - &mut self, - temporal_factor: f64, - axis: Axis, - spatial_factor: f64, - ) -> Result<(), Error> { - self.spatial_stretching_mut(axis, spatial_factor)?; - self.temporal_stretching_mut(temporal_factor)?; - Ok(()) - } + // /// Stretchs this mutable [IONEX] both in temporal and spatial dimensions, + // /// modifying both dimensions. + // /// When factor > 1.0 this is an upscaling operation, when factor < 1.0, + // /// this is a downscaling operation. + // /// This uses both spatial and temporal interpolation to form valid data points. + // pub fn temporal_spatial_stretching_mut( + // &mut self, + // temporal_factor: f64, + // axis: Axis, + // spatial_factor: f64, + // ) -> Result<(), Error> { + // self.spatial_stretching_mut(axis, spatial_factor)?; + // self.temporal_stretching_mut(temporal_factor)?; + // Ok(()) + // } - /// Stretchs this [IONEX] into a new [IONEX], both temporally and spatially stretched, - /// modifying both dimensions. - /// When factor > 1.0 this is an upscaling operation, when factor < 1.0, - /// this is a downscaling operation. - /// This uses both spatial and temporal interpolation to form valid data points. - pub fn temporally_spatially_stretched( - &self, - temporal_factor: f64, - axis: Axis, - spatial_factor: f64, - ) -> Result { - let mut s = self.clone(); - s.temporal_spatial_stretching_mut(temporal_factor, axis, spatial_factor)?; - Ok(s) + // /// Stretchs this [IONEX] into a new [IONEX], both temporally and spatially stretched, + // /// modifying both dimensions. + // /// When factor > 1.0 this is an upscaling operation, when factor < 1.0, + // /// this is a downscaling operation. + // /// This uses both spatial and temporal interpolation to form valid data points. + // pub fn temporally_spatially_stretched( + // &self, + // temporal_factor: f64, + // axis: Axis, + // spatial_factor: f64, + // ) -> Result { + // let mut s = self.clone(); + // s.temporal_spatial_stretching_mut(temporal_factor, axis, spatial_factor)?; + // Ok(s) + // } + + /// Produces a [TimeSeries] from this [IONEX], describing the temporal axis. + pub fn timeseries(&self) -> TimeSeries { + self.header.timeseries() } - /// Designs a [MapCell] iterator (micro rectangle region made of 4 local points) - /// that can then be interpolated. + /// Designs a [MapCell] iterator (micro ROI following the grid quantization) + /// that allows micro interpolation. pub fn map_cell_iter(&self) -> Box + '_> { - let timeseries = self.header.timeseries(); - let lat_pairs = self.header.grid.latitude.quantize().tuple_windows(); let long_pairs = self.header.grid.longitude.quantize().tuple_windows(); @@ -904,7 +896,7 @@ impl IONEX { let fixed_altitude_q = Quantized::auto_scaled(fixed_altitude_km); Box::new( - timeseries + self.timeseries() .cartesian_product(lat_pairs.cartesian_product(long_pairs)) .filter_map(move |(epoch, ((lat1, lat2), (long1, long2)))| { let northeast = Key { @@ -974,6 +966,13 @@ impl IONEX { ) } + // /// Returns a [Cell3x3] iterator for each 3x3 region at a specific point in time + // pub fn synchronous_map_cell3x3_iter(&self, epoch: Epoch) -> Box + '_> { + // Box::new( + // self.map_cell3x3_iter().filter(move |cell| cell.center.epoch == epoch) + // ) + // } + /// Form a possibly interpolated [MapCell] at specified point in time, wrapping this ROI. /// The ROI is described by a [Geometry], which can be complex. /// Unlike [Self::unitary_roi_at], the returned [MapCell] may not be a map quantization cell, diff --git a/src/record/mod.rs b/src/record/mod.rs index 224cd2f..c0fb7fc 100644 --- a/src/record/mod.rs +++ b/src/record/mod.rs @@ -59,13 +59,15 @@ impl Record { )) } - /// Obtain [TEC] value from IONEX [Record], at specified spatial and temporal coordinates that must exist. + /// Obtain [TEC] (single point) from IONEX [Record], at specified spatial and temporal coordinates that must exist. /// This is an indexing method, not an interpolation method. - /// For interpolation, use the [MapCell] API. pub fn get(&self, key: &Key) -> Option<&TEC> { self.map.get(key) } + /// Obtain a [MapCell] (4 single points) from IONEX [Record], at specified point in time and coordinates. + /// The coordinates + /// Obtain mutable [TEC] reference from IONEX [Record], at specified spatial and temporal coordinates that must exist. /// This is an indexing method, not an interpolation method. /// For interpolation, use the [MapCell] API. @@ -73,37 +75,70 @@ impl Record { self.map.get_mut(key) } - /// Collect a IONEX [Record] from a list of [MapCell]. - /// NB: work in progress, not validated yet. - pub fn from_map_cells( - fixed_altitude_km: f64, - min_latitude_ddeg: f64, - max_latitude_ddeg: f64, - min_longitude_ddeg: f64, - max_longitude_ddeg: f64, - cells: &[MapCell], - ) -> Self { - let mut map = Default::default(); - //let (mut new_lat, mut new_long) = (true, true); - //let (mut prev_lat, mut prev_long) = (0.0_f64, 0.0_f64); - - // for cell in cells.iter() { - // SW bound is always introduced - // let sw_key = Key::from_decimal_degrees_km( - // cell.epoch, - // cell.south_west.point.y(), - // cell.south_west.point.x(), - // fixed_altitude_km, - // ); - - // SE bound is introduced for any but first cell - // if !new_lat { - - // } - - //if cell.north_east.point.y() == max_latitude { - //} - // } + /// Collect IONEX [Record] from a list of [MapCell]s. + /// This is particularly useful to reconstruct a [IONEX] file from a possibly processed + /// and modified slice of [MapCell]s. + /// + /// ## Input + /// - slice: slice of [MapCell]s that must have identical dimensions, + /// otherwise this operation will result in corrupt/illegal content, + /// and we do not verify it! + /// - fixed_altitude_km: the fixed altitude in kilometers, + /// use to represent the IONEX plane from the slice of planar [MapCell]s + pub fn from_map_cells(slice: &[MapCell], fixed_altitude_km: f64) -> Self { + let mut map = BTreeMap::::default(); + + for cell in slice.iter() { + // for each cell, we can produce 4 points + // we take advantage of the map to avoid replicated points + let epoch = cell.epoch; + + let (ne_point, nw_point, sw_point, se_point) = ( + cell.north_east.point, + cell.north_west.point, + cell.south_west.point, + cell.south_east.point, + ); + + let (ne_tec, nw_tec, sw_tec, se_tec) = ( + cell.north_east.tec, + cell.north_west.tec, + cell.south_west.tec, + cell.south_east.tec, + ); + + let (ne_key, nw_key, sw_key, se_key) = ( + Key::from_decimal_degrees_km( + cell.epoch, + ne_point.y(), + ne_point.x(), + fixed_altitude_km, + ), + Key::from_decimal_degrees_km( + cell.epoch, + nw_point.y(), + nw_point.x(), + fixed_altitude_km, + ), + Key::from_decimal_degrees_km( + cell.epoch, + sw_point.y(), + sw_point.x(), + fixed_altitude_km, + ), + Key::from_decimal_degrees_km( + cell.epoch, + se_point.y(), + se_point.x(), + fixed_altitude_km, + ), + ); + + map.insert(ne_key, ne_tec); + map.insert(nw_key, nw_tec); + map.insert(sw_key, sw_tec); + map.insert(se_key, se_tec); + } Self { map } } From ee074791c263a56f99204c21bacae8d9a619bb70 Mon Sep 17 00:00:00 2001 From: "Guillaume W. Bres" Date: Tue, 11 Nov 2025 14:29:56 +0100 Subject: [PATCH 25/25] up Signed-off-by: Guillaume W. Bres --- src/cell/mod.rs | 5 +---- src/lib.rs | 44 +++++++++++++++++++++++++++++++++++++++----- src/record/mod.rs | 18 ++---------------- src/tests/mod.rs | 2 +- src/tests/v1.rs | 12 ++++++++++++ 5 files changed, 55 insertions(+), 26 deletions(-) diff --git a/src/cell/mod.rs b/src/cell/mod.rs index 4f3c232..f1eeab9 100644 --- a/src/cell/mod.rs +++ b/src/cell/mod.rs @@ -2,10 +2,7 @@ use std::collections::HashMap; use geo::{Contains, GeodesicArea, Geometry, Point, Rect}; -use crate::{ - prelude::{Epoch, Error, TEC}, - rectangle_to_cardinals, -}; +use crate::prelude::{Epoch, Error, TEC}; mod three_by_three; pub use three_by_three::Cell3x3; diff --git a/src/lib.rs b/src/lib.rs index aa1b55f..5c83adc 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -116,7 +116,7 @@ fn div_ceil(value: usize, divider: usize) -> usize { } /// Converts a geo [Rect]angle to NE, SE, SW, NW (latitude, longitude) quadruplets -pub(crate) fn rectangle_to_cardinals( +pub(crate) fn rectangle_quadrant_decomposition( rect: Rect, ) -> ((f64, f64), (f64, f64), (f64, f64), (f64, f64)) { let (min, width, height) = (rect.min(), rect.width(), rect.height()); @@ -129,6 +129,15 @@ pub(crate) fn rectangle_to_cardinals( ) } +/// Converts a quadruplet (NE, SE, SW, NW) (latitude, longitude) coordinates to a [Rect]angle in degrees +pub(crate) fn quadrant_to_rectangle( + quadrant: ((f64, f64), (f64, f64), (f64, f64), (f64, f64)), +) -> Rect { + let (ne_lat, ne_long) = quadrant.0; + let (se_lat, se_long) = quadrant.1; + Rect::new(coord!(x: se_long, y: se_lat), coord!(x: ne_long, y: ne_lat)) +} + /// macro to format one header line or a comment pub(crate) fn fmt_ionex(content: &str, marker: &str) -> String { if content.len() < 60 { @@ -598,11 +607,28 @@ impl IONEX { /// The quantization (both spatial and temporal) is preserved, only the /// spatial dimensions are modified by this operation. /// [Polygon::bounding_rect] must be defined for this operation to work correctly. + /// + /// Assuming the original map and '*' marking the complex ROI you want to draw, + /// the resuling IONEX is the smallest rectangle the ROI fits in, because this library + /// (and correct/valid IONEX) is limited to rectangles. + /// + /// origin + /// ---------------------------------- + /// | returned roi___ | + /// | | * * * | | + /// | | * * | | + /// | | * specs * | | + /// | | * * | | + /// | | ****** * | | + /// | ______________| | + /// |---------------------------------| pub fn to_regional_ionex(&self, roi: Polygon) -> Result { let mut ionex = IONEX::default().with_header(self.header.clone()); + let boundaries = roi.bounding_rect().ok_or(Error::UndefinedBoundaries)?; + // simply restrict to wrapping area - let roi = Geometry::Polygon(roi); + let roi = Geometry::Rect(boundaries); let cells = self .map_cell_iter() @@ -622,6 +648,14 @@ impl IONEX { }, } + // update header + ionex.header.grid.latitude.start = self.header.grid.latitude.start.min(boundaries.max().y); + ionex.header.grid.longitude.start = + self.header.grid.longitude.start.max(boundaries.min().x); + + ionex.header.grid.latitude.end = self.header.grid.latitude.end.max(boundaries.min().y); + ionex.header.grid.longitude.end = self.header.grid.longitude.end.min(boundaries.max().x); + Ok(ionex) } @@ -1294,7 +1328,7 @@ impl gnss_qc_traits::Merge for IONEX { #[cfg(test)] mod test { - use crate::{div_ceil, fmt_comment, prelude::*, rectangle_to_cardinals}; + use crate::{div_ceil, fmt_comment, prelude::*, rectangle_quadrant_decomposition}; #[test] fn fmt_comments_singleline() { @@ -1335,7 +1369,7 @@ mod test { } #[test] - fn rect_to_cardinals() { + fn rectangle_decomposition() { for (rect, ((lat11, long11), (lat12, long12), (lat21, long21), (lat22, long22))) in [ ( Rect::new(coord!(x: -30.0, y: -30.0), coord!(x: 30.0, y: 30.0)), @@ -1357,7 +1391,7 @@ mod test { ), ] { assert_eq!( - rectangle_to_cardinals(rect), + rectangle_quadrant_decomposition(rect), ( (lat11, long11), (lat12, long12), diff --git a/src/record/mod.rs b/src/record/mod.rs index c0fb7fc..30dfcc8 100644 --- a/src/record/mod.rs +++ b/src/record/mod.rs @@ -169,14 +169,7 @@ mod test { let map_cells = ionex.map_cell_iter().collect::>(); // build from scratch - let record = Record::from_map_cells( - ionex.header.grid.altitude.start, - ionex.header.grid.latitude.end, - ionex.header.grid.latitude.start, - ionex.header.grid.longitude.start, - ionex.header.grid.longitude.end, - &map_cells, - ); + let record = Record::from_map_cells(&map_cells, ionex.header.grid.altitude.start); // reciprocal assert_eq!(record, ionex.record); @@ -193,14 +186,7 @@ mod test { let map_cells = ionex.map_cell_iter().collect::>(); // build from scratch - let record = Record::from_map_cells( - ionex.header.grid.altitude.start, - ionex.header.grid.latitude.end, - ionex.header.grid.latitude.start, - ionex.header.grid.longitude.start, - ionex.header.grid.longitude.end, - &map_cells, - ); + let record = Record::from_map_cells(&map_cells, ionex.header.grid.altitude.start); // reciprocal assert_eq!(record, ionex.record); diff --git a/src/tests/mod.rs b/src/tests/mod.rs index b77be60..2ec56a5 100644 --- a/src/tests/mod.rs +++ b/src/tests/mod.rs @@ -5,7 +5,7 @@ mod filename; mod parsing; mod qc; mod roi; -mod stretching; +// mod stretching; mod v1; diff --git a/src/tests/v1.rs b/src/tests/v1.rs index c3e505b..03d1274 100644 --- a/src/tests/v1.rs +++ b/src/tests/v1.rs @@ -20,6 +20,12 @@ fn parse_ckmg0020() { assert!(ionex.is_2d()); assert!(!ionex.is_3d()); + // worldwide boundaries + assert_eq!( + ionex.bounding_rect_degrees(), + Rect::new(coord!(x: -180.0, y: -87.5), coord!(x: 180.0, y: 87.5),) + ); + let testpoints = vec![ TestPoint { epoch_str: "2022-01-02T00:00:00 UTC", @@ -268,6 +274,12 @@ fn parse_jplg() { assert!(ionex.is_2d()); assert!(!ionex.is_3d()); + // worldwide boundaries + assert_eq!( + ionex.bounding_rect_degrees(), + Rect::new(coord!(x: -180.0, y: -87.5), coord!(x: 180.0, y: 87.5),) + ); + let testpoints = vec![ TestPoint { epoch_str: "2017-01-01T00:00:00 UTC",