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/Cargo.toml b/Cargo.toml index 58b637f..6feb7cd 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -22,28 +22,49 @@ 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", +] + +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"] [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"] } +maud = { version = "0.26", optional = true } 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 +77,4 @@ harness = false [[bench]] name = "formatting" harness = false + diff --git a/README.md b/README.md index dcdb53b..68f719e 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 @@ -75,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/cell.rs b/src/cell.rs deleted file mode 100644 index 3c96a4a..0000000 --- a/src/cell.rs +++ /dev/null @@ -1,520 +0,0 @@ -use geo::{Contains, GeodesicArea, Geometry, Point, Rect}; - -use crate::prelude::{Epoch, TEC}; - -#[derive(Debug, Default, Clone, Copy, PartialEq)] -pub struct MapPoint { - /// [Point] - pub point: Point, - - /// TEC - pub tec: TEC, -} - -/// [MapCell] describing a region that we can then interpolate. -#[derive(Debug, Default, Clone, Copy, PartialEq)] -pub struct MapCell { - /// Epoch of observation - pub epoch: Epoch, - - /// North East [MapPoint] - pub north_east: MapPoint, - - /// North West [MapPoint] - pub north_west: MapPoint, - - /// South East [MapPoint] - pub south_east: MapPoint, - - /// South West [MapPoint] - pub south_west: MapPoint, -} - -impl MapCell { - /// Define a new [MapCell] from 4 (latitude_ddeg, longitude_ddeg) cardinal tuples and - /// associated TEC values. - pub fn from_lat_long_degrees( - epoch: Epoch, - northeast_ddeg: (f64, f64), - northeast_tec: TEC, - northwest_ddeg: (f64, f64), - northwest_tec: TEC, - southeast_ddeg: (f64, f64), - southeast_tec: TEC, - southwest_ddeg: (f64, f64), - southwest_tec: TEC, - ) -> Self { - Self { - epoch, - north_east: MapPoint { - point: Point::new(northeast_ddeg.0, northeast_ddeg.1), - tec: northeast_tec, - }, - north_west: MapPoint { - point: Point::new(northwest_ddeg.0, northwest_ddeg.1), - tec: northwest_tec, - }, - south_east: MapPoint { - point: Point::new(southeast_ddeg.0, southeast_ddeg.1), - tec: southeast_tec, - }, - south_west: MapPoint { - point: Point::new(southwest_ddeg.0, southwest_ddeg.1), - tec: southwest_tec, - }, - } - } - - /// Define a new [MapCell] from 4 (latitude_rad, longitude_rad) cardinal tuples and - /// associated TEC values. - pub fn from_lat_long_radians( - epoch: Epoch, - northeast_rad: (f64, f64), - northeast_tec: TEC, - northwest_rad: (f64, f64), - northwest_tec: TEC, - southeast_rad: (f64, f64), - southeast_tec: TEC, - southwest_rad: (f64, f64), - southwest_tec: TEC, - ) -> Self { - Self { - epoch, - north_east: MapPoint { - point: Point::new(northeast_rad.0.to_degrees(), northeast_rad.1.to_degrees()), - tec: northeast_tec, - }, - north_west: MapPoint { - point: Point::new(northwest_rad.0.to_degrees(), northwest_rad.1.to_degrees()), - tec: northwest_tec, - }, - south_east: MapPoint { - point: Point::new(southeast_rad.0.to_degrees(), southeast_rad.1.to_degrees()), - tec: southeast_tec, - }, - south_west: MapPoint { - 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]. - pub fn from_cardinal_points( - epoch: Epoch, - north_east: MapPoint, - north_west: MapPoint, - south_east: MapPoint, - south_west: MapPoint, - ) -> Self { - Self { - epoch, - north_east, - north_west, - south_east, - south_west, - } - } - - /// 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 - /// - (x/long=1, y/lat=0) is the SE corner - /// - (x/long=0, y/lat=1) is the NW corner - /// - (x/long=1, y/lat=1) is the NE corner - /// - /// ``` - /// use ionex::prelude::{MapCell, Epoch, TEC}; - /// - /// let epoch = Epoch::default(); - /// let tec = TEC::from_tecu(1.0); - /// let cell = MapCell::from_unitary_tec(epoch, tec, tec, tec, tec); - /// - /// // 1.0° unitary ECEF span - /// assert_eq!(cell.latitude_longitude_span_degrees(), (1.0, 1.0)); - /// - /// // ECEF! - /// assert!((cell.geodesic_perimeter() - 443770.0).abs() < 1.0); - /// assert!((cell.geodesic_area() - 12308778361.0).abs() < 1.0); - /// ``` - pub fn from_unitary_tec( - epoch: Epoch, - northeast_tec: TEC, - northwest_tec: TEC, - southeast_tec: TEC, - southwest_tec: TEC, - ) -> Self { - Self::from_lat_long_degrees( - epoch, - (1.0, 1.0), - northeast_tec, - (0.0, 1.0), - northwest_tec, - (1.0, 0.0), - southeast_tec, - (0.0, 0.0), - southwest_tec, - ) - } - - /// Defines a unitary [MapCell] ((0,0), (0,1), (1,0), (1,1)) with Null TEC values, - /// where - /// - (x/long=0, y/lat=0) is the SW corner - /// - (x/long=1, y/lat=0) is the SE corner - /// - (x/long=0, y/lat=1) is the NW corner - /// - (x/long=1, y/lat=1) is the NE corner - /// - /// ``` - /// use ionex::prelude::{MapCell, Epoch, TEC}; - /// - /// let epoch = Epoch::default(); - /// let tec = TEC::from_tecu(1.0); - /// let cell = MapCell::unitary_null_tec(epoch); - /// - /// // Null values! - /// assert_eq!(cell.north_east.tec.tecu(), 0.0); - /// assert_eq!(cell.north_west.tec.tecu(), 0.0); - /// assert_eq!(cell.south_east.tec.tecu(), 0.0); - /// assert_eq!(cell.south_west.tec.tecu(), 0.0); - /// - /// // 1.0° unitary ECEF span - /// assert_eq!(cell.latitude_longitude_span_degrees(), (1.0, 1.0)); - /// - /// // ECEF! - /// assert!((cell.geodesic_perimeter() - 443770.0).abs() < 1.0); - /// assert!((cell.geodesic_area() - 12308778361.0).abs() < 1.0); - /// ``` - pub fn unitary_null_tec(epoch: Epoch) -> Self { - let null_tec = TEC::default(); - Self::from_unitary_tec(epoch, null_tec, null_tec, null_tec, null_tec) - } - - /// Returns central [Point] of this [MapCell]. - pub fn center(&self) -> Point { - geo::Point(self.borders().center()) - } - - /// Returns borders of this [MapCell] expressed as a [Rect]angle. - pub fn borders(&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() - } - - /// Returns geodesic area (in squared meters) of this [MapCell]. - pub fn geodesic_area(&self) -> f64 { - self.borders().geodesic_area_unsigned() - } - - /// Returns true if following [Geometry] is contained within this [MapCell]. - pub fn contains(&self, geometry: &Geometry) -> bool { - self.borders().contains(geometry) - } - - /// Copies and updates the Northeastern TEC component - pub fn with_northeastern_tec(mut self, tec: TEC) -> Self { - self.north_east.tec = tec; - self - } - - /// Copies and updates the Northwestern TEC component - pub fn with_northwestern_tec(mut self, tec: TEC) -> Self { - self.north_west.tec = tec; - self - } - - /// Copies and updates the Southeastern TEC component - pub fn with_southeastern_tec(mut self, tec: TEC) -> Self { - self.south_east.tec = tec; - self - } - - /// Copies and updates the Southwestern TEC component - pub fn with_southwestern_tec(mut self, tec: TEC) -> Self { - self.south_west.tec = tec; - self - } - - /// Returns the (latitude, longitude) span of this [MapCell] - /// as tuplet in degrees - pub fn latitude_longitude_span_degrees(&self) -> (f64, f64) { - (self.latitude_span_degrees(), self.longitude_span_degrees()) - } - - /// Returns latitude span of this [MapCell] in degrees - pub fn latitude_span_degrees(&self) -> f64 { - let borders = self.borders(); - 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(); - borders.max().x - borders.min().x - } - - /// Spatial interpolation of the [TEC] value using planery equation - /// and 4 boundaries of this [MapCell]. [MapCell::contains] should be true - /// for the proposed geometry for this to be correct. - /// This method does not verify this assertion, it is up to you to use valid coordinates here. - /// - /// Example: unitary cell - /// ``` - /// use ionex::prelude::{MapCell, Epoch, Point, TEC, Unit}; - /// - /// // create unitary cell with simple values - /// let t0 = Epoch::default(); - /// let one_tec = TEC::from_tecu(1.0); - /// - /// let cell = MapCell::from_unitary_tec(t0, one_tec, one_tec, one_tec, one_tec); - /// - /// // central point - /// let center = Point::new(0.5, 0.5); - /// let tec = cell.spatial_interpolation(center); - /// assert_eq!(tec.tecu(), 1.0); - /// ``` - /// - /// Example: North East gradient - /// ``` - /// use ionex::prelude::{MapCell, Epoch, Point, TEC, Unit}; - /// - /// // create unitary cell with simple values - /// let t0 = Epoch::default(); - /// - /// let gradient = ( - /// TEC::from_tecu(0.0), - /// TEC::from_tecu(0.0), - /// TEC::from_tecu(0.0), - /// TEC::from_tecu(1.0), - /// ); - /// - /// 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)); - /// assert_eq!(tec.tecu(), 0.25); - /// - /// // SW boundary - /// let tec = cell.spatial_interpolation(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)); - /// assert_eq!(tec.tecu(), 0.81); - /// - /// // SWwern point - /// 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 { - let (latitude_span, longitude_span) = self.latitude_longitude_span_degrees(); - - let (p, q) = (point.y() / latitude_span, point.x() / longitude_span); - - let (e00, e10, e01, e11) = ( - self.south_west.tec.tecu(), - self.south_east.tec.tecu(), - self.north_west.tec.tecu(), - self.north_east.tec.tecu(), - ); - - 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) - } - - /// 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 - /// for the results to be correct, but this is not verified here: it is up to you - /// to use valid coordinates here. - /// Proposed [Epoch] should lie within both observation instants, otherwise this method - /// returns None. - /// - /// ``` - /// use ionex::prelude::{MapCell, Epoch, Point, TEC, Unit}; - /// - /// // create two unitary cells with simple values - /// let t0 = Epoch::default(); - /// let t1 = t0 + 30.0 * Unit::Second; - /// let t_ok = t0 + 15.0 * Unit::Second; // within interval - /// let t_nok = t0 + 45.0 * Unit::Second; // outside interval - /// - /// let one_tec = TEC::from_tecu(1.0); - /// - /// let center = Point::new(0.5, 0.5); // unitary cell - /// let cell0 = MapCell::from_unitary_tec(t0, one_tec, one_tec, one_tec, one_tec); - /// 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); - /// assert_eq!(central_tec0.tecu(), 1.0); - /// - /// // verify central point value - /// let central_tec1 = cell1.spatial_interpolation(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()); - /// - /// // spatial + temporal interpolation - /// let tec = cell0.temporal_spatial_interpolation(t_ok, center, &cell1).unwrap(); - /// assert_eq!(tec.tecu(), 1.0); - /// ``` - pub fn temporal_spatial_interpolation( - &self, - epoch: Epoch, - point: Point, - rhs: &Self, - ) -> Option { - // interpolate at exact coordinates - let (tecu_0, tecu_1) = ( - self.spatial_interpolation(point).tecu(), - rhs.spatial_interpolation(point).tecu(), - ); - - if epoch >= self.epoch && epoch < rhs.epoch { - // forward - let dt = (rhs.epoch - self.epoch).to_seconds(); - - let tecu = (rhs.epoch - epoch).to_seconds() / dt * tecu_0 - + (epoch - self.epoch).to_seconds() / dt * tecu_1; - - Some(TEC::from_tecu(tecu)) - } else if epoch >= rhs.epoch && epoch < self.epoch { - // backwards - let dt = (self.epoch - rhs.epoch).to_seconds(); - - let tecu = (self.epoch - epoch).to_seconds() / dt * tecu_1 - + (epoch - rhs.epoch).to_seconds() / dt * tecu_0; - - Some(TEC::from_tecu(tecu)) - } else { - None - } - } -} - -#[cfg(test)] -mod test { - use super::*; - - use crate::prelude::{Epoch, Geometry, Point, Unit, TEC}; - - #[test] - fn spatial_unitary_interpolation() { - let epoch = Epoch::default(); - - let northeast_tec = TEC::from_tecu(1.0); - let northwest_tec = TEC::from_tecu(1.0); - let southeast_tec = TEC::from_tecu(1.0); - let southwest_tec = TEC::from_tecu(1.0); - - let cell = MapCell::from_unitary_tec( - epoch, - northeast_tec, - northwest_tec, - southeast_tec, - southwest_tec, - ); - - assert_eq!(cell.latitude_longitude_span_degrees(), (1.0, 1.0)); - - assert!((cell.geodesic_perimeter() - 443770.0).abs() < 1.0); - assert!((cell.geodesic_area() - 12308778361.0).abs() < 1.0); - - let center = Point::new(0.5, 0.5); - let outside = Point::new(1.5, 0.5); - - assert!(cell.contains(&Geometry::Point(center))); - assert!(!cell.contains(&Geometry::Point(outside))); - assert_eq!(center, cell.center()); - - let interpolated = cell.spatial_interpolation(center); - assert_eq!(interpolated.tecu(), 1.0); - } - - #[test] - fn spatial_south_west_gradient_interpolation() { - let epoch = Epoch::default(); - - let northeast_tec = TEC::from_tecu(0.0); - let northwest_tec = TEC::from_tecu(0.0); - let southeast_tec = TEC::from_tecu(0.0); - let southwest_tec = TEC::from_tecu(1.0); - - let cell = MapCell::from_unitary_tec( - epoch, - northeast_tec, - northwest_tec, - southeast_tec, - southwest_tec, - ); - - for (x_deg, y_deg, tecu) in [ - (0.5, 0.5, 0.25), - (0.1, 0.1, 0.81), - (0.01, 0.01, 0.9801), - (0.0, 0.0, 1.0), - ] { - let point = Point::new(x_deg, y_deg); - let interpolated = cell.spatial_interpolation(point).tecu(); - assert_eq!(interpolated, tecu, "failed at (x={}, y={})", x_deg, y_deg); - } - } - - #[test] - fn temporal_interpolation() { - let t0 = Epoch::default(); - let t1 = t0 + 1.0 * Unit::Second; - - let t_ok = t0 + 0.5 * Unit::Second; - let t_nok = t1 + 2.0 * Unit::Second; - - let center = Point::new(0.5, 0.5); - - let northeast_tec_0 = TEC::from_tecu(1.0); - let northwest_tec_0 = TEC::from_tecu(1.0); - let southeast_tec_0 = TEC::from_tecu(1.0); - let southwest_tec_0 = TEC::from_tecu(1.0); - - let cell0 = MapCell::from_unitary_tec( - t0, - northeast_tec_0, - northwest_tec_0, - southeast_tec_0, - southwest_tec_0, - ); - - let northeast_tec_1 = TEC::from_tecu(1.0); - let northwest_tec_1 = TEC::from_tecu(1.0); - let southeast_tec_1 = TEC::from_tecu(1.0); - let southwest_tec_1 = TEC::from_tecu(1.0); - - let cell1 = MapCell::from_unitary_tec( - t1, - northeast_tec_1, - northwest_tec_1, - southeast_tec_1, - southwest_tec_1, - ); - - assert!( - cell0 - .temporal_spatial_interpolation(t_nok, center, &cell1) - .is_none(), - "interpolation is temporally incorrect!" - ); - - let tec = cell0 - .temporal_spatial_interpolation(t_ok, center, &cell1) - .unwrap(); - - assert_eq!(tec.tecu(), 1.0); - } -} diff --git a/src/cell/mod.rs b/src/cell/mod.rs new file mode 100644 index 0000000..f1eeab9 --- /dev/null +++ b/src/cell/mod.rs @@ -0,0 +1,962 @@ +use std::collections::HashMap; + +use geo::{Contains, GeodesicArea, Geometry, Point, Rect}; + +use crate::prelude::{Epoch, Error, TEC}; + +mod three_by_three; +pub use three_by_three::Cell3x3; + +#[derive(Debug, Default, Clone, Copy, PartialEq)] +pub struct TecPoint { + /// TEC + pub tec: TEC, + + /// [Point] + pub point: Point, +} + +/// [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 + pub epoch: Epoch, + + /// 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, +} + +impl MapCell { + /// Define a new [MapCell] from 4 (latitude_ddeg, longitude_ddeg) cardinal tuples and + /// associated TEC values. + pub fn from_lat_long_degrees( + epoch: Epoch, + northeast_ddeg: (f64, f64), + northeast_tec: TEC, + northwest_ddeg: (f64, f64), + northwest_tec: TEC, + southeast_ddeg: (f64, f64), + southeast_tec: TEC, + southwest_ddeg: (f64, f64), + southwest_tec: TEC, + ) -> Self { + Self { + epoch, + north_east: TecPoint { + point: Point::new(northeast_ddeg.0, northeast_ddeg.1), + tec: northeast_tec, + }, + north_west: TecPoint { + point: Point::new(northwest_ddeg.0, northwest_ddeg.1), + tec: northwest_tec, + }, + south_east: TecPoint { + point: Point::new(southeast_ddeg.0, southeast_ddeg.1), + tec: southeast_tec, + }, + south_west: TecPoint { + point: Point::new(southwest_ddeg.0, southwest_ddeg.1), + tec: southwest_tec, + }, + } + } + + /// Define a new [MapCell] from 4 (latitude_rad, longitude_rad) cardinal tuples and + /// associated TEC values. + pub fn from_lat_long_radians( + epoch: Epoch, + northeast_rad: (f64, f64), + northeast_tec: TEC, + northwest_rad: (f64, f64), + northwest_tec: TEC, + southeast_rad: (f64, f64), + southeast_tec: TEC, + southwest_rad: (f64, f64), + southwest_tec: TEC, + ) -> Self { + Self { + epoch, + north_east: TecPoint { + point: Point::new(northeast_rad.0.to_degrees(), northeast_rad.1.to_degrees()), + tec: northeast_tec, + }, + north_west: TecPoint { + point: Point::new(northwest_rad.0.to_degrees(), northwest_rad.1.to_degrees()), + tec: northwest_tec, + }, + south_east: TecPoint { + point: Point::new(southeast_rad.0.to_degrees(), southeast_rad.1.to_degrees()), + tec: southeast_tec, + }, + 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 [TecPoint]s describing each corner at this [Epoch]. + pub fn from_cardinal_points( + epoch: Epoch, + north_east: TecPoint, + north_west: TecPoint, + south_east: TecPoint, + south_west: TecPoint, + ) -> Self { + Self { + epoch, + north_east, + north_west, + south_east, + south_west, + } + } + + /// 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 + } + + /// 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 + && rhs.north_west.point == self.south_west.point + } + + 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 + } + + 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]. + 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_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] 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())) + } + + /// 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 + /// - (x/long=1, y/lat=0) is the SE corner + /// - (x/long=0, y/lat=1) is the NW corner + /// - (x/long=1, y/lat=1) is the NE corner + /// + /// ``` + /// use ionex::prelude::{MapCell, Epoch, TEC}; + /// + /// let epoch = Epoch::default(); + /// let tec = TEC::from_tecu(1.0); + /// let cell = MapCell::from_unitary_tec(epoch, tec, tec, tec, tec); + /// + /// // 1.0° unitary ECEF span + /// assert_eq!(cell.latitude_longitude_span_degrees(), (1.0, 1.0)); + /// + /// // ECEF! + /// assert!((cell.geodesic_perimeter() - 443770.0).abs() < 1.0); + /// assert!((cell.geodesic_area() - 12308778361.0).abs() < 1.0); + /// ``` + pub fn from_unitary_tec( + epoch: Epoch, + northeast_tec: TEC, + northwest_tec: TEC, + southeast_tec: TEC, + southwest_tec: TEC, + ) -> Self { + Self::from_lat_long_degrees( + epoch, + (1.0, 1.0), + northeast_tec, + (0.0, 1.0), + northwest_tec, + (1.0, 0.0), + southeast_tec, + (0.0, 0.0), + southwest_tec, + ) + } + + /// Defines a unitary [MapCell] ((0,0), (0,1), (1,0), (1,1)) with Null TEC values, + /// where + /// - (x/long=0, y/lat=0) is the SW corner + /// - (x/long=1, y/lat=0) is the SE corner + /// - (x/long=0, y/lat=1) is the NW corner + /// - (x/long=1, y/lat=1) is the NE corner + /// + /// ``` + /// use ionex::prelude::{MapCell, Epoch, TEC}; + /// + /// let epoch = Epoch::default(); + /// let tec = TEC::from_tecu(1.0); + /// let cell = MapCell::unitary_null_tec(epoch); + /// + /// // Null values! + /// assert_eq!(cell.north_east.tec.tecu(), 0.0); + /// assert_eq!(cell.north_west.tec.tecu(), 0.0); + /// assert_eq!(cell.south_east.tec.tecu(), 0.0); + /// assert_eq!(cell.south_west.tec.tecu(), 0.0); + /// + /// // 1.0° unitary ECEF span + /// assert_eq!(cell.latitude_longitude_span_degrees(), (1.0, 1.0)); + /// + /// // ECEF! + /// assert!((cell.geodesic_perimeter() - 443770.0).abs() < 1.0); + /// assert!((cell.geodesic_area() - 12308778361.0).abs() < 1.0); + /// ``` + pub fn unitary_null_tec(epoch: Epoch) -> Self { + let null_tec = TEC::default(); + Self::from_unitary_tec(epoch, null_tec, null_tec, null_tec, null_tec) + } + + /// Returns central [Point] of this [MapCell]. + pub fn center(&self) -> Point { + geo::Point(self.bounding_rect_degrees().center()) + } + + /// 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) + } + + /// Returns geodesic perimeter (in meters) of this [MapCell]. + pub fn geodesic_perimeter(&self) -> f64 { + self.bounding_rect_degrees().geodesic_perimeter() + } + + /// Returns geodesic area (in squared meters) of this [MapCell]. + pub fn geodesic_area(&self) -> f64 { + self.bounding_rect_degrees().geodesic_area_unsigned() + } + + /// Returns true if following [Geometry], expressed in decimal degrees, + /// is contained within this [MapCell]. + pub fn contains(&self, geometry: &Geometry) -> bool { + self.bounding_rect_degrees().contains(geometry) + } + + /// Copies and updates the Northeastern TEC component + pub fn with_northeastern_tec(mut self, tec: TEC) -> Self { + self.north_east.tec = tec; + self + } + + /// Copies and updates the Northwestern TEC component + pub fn with_northwestern_tec(mut self, tec: TEC) -> Self { + self.north_west.tec = tec; + self + } + + /// Copies and updates the Southeastern TEC component + pub fn with_southeastern_tec(mut self, tec: TEC) -> Self { + self.south_east.tec = tec; + self + } + + /// Copies and updates the Southwestern TEC component + pub fn with_southwestern_tec(mut self, tec: TEC) -> Self { + self.south_west.tec = tec; + 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) { + (self.latitude_span_degrees(), self.longitude_span_degrees()) + } + + /// Returns latitude span of this [MapCell] in degrees + pub fn latitude_span_degrees(&self) -> f64 { + 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.bounding_rect_degrees(); + borders.max().x - borders.min().x + } + + /// Spatial interpolation of the [TEC] value using planery equation + /// and 4 boundaries of this [MapCell]. [MapCell::contains] should be true + /// for the proposed geometry for this to be correct. + /// This method does not verify this assertion, it is up to you to use valid coordinates here. + /// + /// Example: unitary cell + /// ``` + /// use ionex::prelude::{MapCell, Epoch, Point, TEC, Unit}; + /// + /// // create unitary cell with simple values + /// let t0 = Epoch::default(); + /// let one_tec = TEC::from_tecu(1.0); + /// + /// let cell = MapCell::from_unitary_tec(t0, one_tec, one_tec, one_tec, one_tec); + /// + /// // central point + /// let center = Point::new(0.5, 0.5); + /// let tec = cell.spatial_tec_interp(center); + /// assert_eq!(tec.tecu(), 1.0); + /// ``` + /// + /// Example: North East gradient + /// ``` + /// use ionex::prelude::{MapCell, Epoch, Point, TEC, Unit}; + /// + /// // create unitary cell with simple values + /// let t0 = Epoch::default(); + /// + /// let gradient = ( + /// TEC::from_tecu(0.0), + /// TEC::from_tecu(0.0), + /// TEC::from_tecu(0.0), + /// TEC::from_tecu(1.0), + /// ); + /// + /// let cell = MapCell::from_unitary_tec(t0, gradient.0, gradient.1, gradient.2, gradient.3); + /// + /// // central point + /// let tec = cell.spatial_tec_interp(Point::new(0.5, 0.5)); + /// assert_eq!(tec.tecu(), 0.25); + /// + /// // SW boundary + /// let tec = cell.spatial_tec_interp(Point::new(0.0, 0.0)); + /// assert_eq!(tec.tecu(), 1.0); + /// + /// // SWern point + /// let tec = cell.spatial_tec_interp(Point::new(0.1, 0.1)); + /// assert_eq!(tec.tecu(), 0.81); + /// + /// // SWwern point + /// 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 { + if !self.contains(&Geometry::Point(point)) { + return Err(Error::OutsideSpatialBoundaries); + } + + let (latitude_span, longitude_span) = self.latitude_longitude_span_degrees(); + + let (p, q) = (point.y() / latitude_span, point.x() / longitude_span); + + let (e00, e10, e01, e11) = ( + self.south_west.tec.tecu(), + self.south_east.tec.tecu(), + self.north_west.tec.tecu(), + self.north_east.tec.tecu(), + ); + + let tecu = + (1.0 - p) * (1.0 - q) * e00 + p * (1.0 - q) * e10 + q * (1.0 - p) * e01 + p * q * e11; + + Ok(TEC::from_tecu(tecu)) + } + + /// 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); + } + + // 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, + ), + ); + + let (north_east, north_west, south_east, south_west) = ( + TecPoint { + point: north_east, + tec: self.spatial_tec_interp(north_east)?, + }, + TecPoint { + point: north_west, + tec: self.spatial_tec_interp(north_west)?, + }, + TecPoint { + point: south_east, + tec: self.spatial_tec_interp(south_east)?, + }, + TecPoint { + point: south_west, + tec: self.spatial_tec_interp(south_west)?, + }, + ); + + self.north_east = north_east; + self.north_west = north_west; + self.south_east = south_east; + self.south_west = south_west; + + Ok(()) + } + + /// 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 + // 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) { + // 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, + // /// 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: TecPoint { + // point: self.north_east.point, + // tec: TEC::from_tecu( + // num_1 * ne_1 /dt + num_2 * ne_2 /dt + // ), + // }, + // north_west: TecPoint { + // point: self.north_west.point, + // tec: TEC::from_tecu( + // num_1 * nw_1 /dt + num_2 * nw_2 /dt + // ), + // }, + // south_east: TecPoint { + // point: self.south_east.point, + // tec: TEC::from_tecu( + // num_1 * se_1 /dt + num_2 * se_2 /dt + // ), + // }, + // south_west: TecPoint { + // 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 + /// for the results to be correct, but this is not verified here: it is up to you + /// to use valid coordinates here. + /// Proposed [Epoch] should lie within both observation instants, otherwise this method + /// returns None. + /// + /// ``` + /// use ionex::prelude::{MapCell, Epoch, Point, TEC, Unit}; + /// + /// // create two unitary cells with simple values + /// let t0 = Epoch::default(); + /// let t1 = t0 + 30.0 * Unit::Second; + /// let t_ok = t0 + 15.0 * Unit::Second; // within interval + /// let t_nok = t0 + 45.0 * Unit::Second; // outside interval + /// + /// let one_tec = TEC::from_tecu(1.0); + /// + /// let center = Point::new(0.5, 0.5); // unitary cell + /// let cell0 = MapCell::from_unitary_tec(t0, one_tec, one_tec, one_tec, one_tec); + /// let cell1 = MapCell::from_unitary_tec(t1, one_tec, one_tec, one_tec, one_tec); + /// + /// // verify central point value + /// let central_tec0 = cell0.spatial_tec_interp(center); + /// assert_eq!(central_tec0.tecu(), 1.0); + /// + /// // verify central point value + /// 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_tec_interp(t_nok, center, &cell1).is_none()); + /// + /// // spatial + temporal interpolation + /// let tec = cell0.temporal_spatial_tec_interp(t_ok, center, &cell1).unwrap(); + /// assert_eq!(tec.tecu(), 1.0); + /// ``` + pub fn temporal_spatial_tec_interp( + &self, + epoch: Epoch, + coordinates: Point, + rhs: &Self, + ) -> Result { + // interpolate at exact coordinates + let (tecu_0, tecu_1) = ( + self.spatial_tec_interp(coordinates)?.tecu(), + rhs.spatial_tec_interp(coordinates)?.tecu(), + ); + + if epoch >= self.epoch && epoch < rhs.epoch { + // forward + let dt = (rhs.epoch - self.epoch).to_seconds(); + + let tecu = (rhs.epoch - epoch).to_seconds() / dt * tecu_0 + + (epoch - self.epoch).to_seconds() / dt * tecu_1; + + Ok(TEC::from_tecu(tecu)) + } else if epoch >= rhs.epoch && epoch < self.epoch { + // backwards + let dt = (self.epoch - rhs.epoch).to_seconds(); + + let tecu = (self.epoch - epoch).to_seconds() / dt * tecu_1 + + (epoch - rhs.epoch).to_seconds() / dt * tecu_0; + + Ok(TEC::from_tecu(tecu)) + } else { + Err(Error::TemporalMismatch) + } + } +} + +#[cfg(test)] +mod test { + use super::*; + + use crate::prelude::{Epoch, Geometry, Point, Unit, TEC}; + + #[test] + fn spatial_unitary_interpolation() { + let epoch = Epoch::default(); + + let northeast_tec = TEC::from_tecu(1.0); + let northwest_tec = TEC::from_tecu(1.0); + let southeast_tec = TEC::from_tecu(1.0); + let southwest_tec = TEC::from_tecu(1.0); + + let cell = MapCell::from_unitary_tec( + epoch, + northeast_tec, + northwest_tec, + southeast_tec, + southwest_tec, + ); + + assert_eq!(cell.latitude_longitude_span_degrees(), (1.0, 1.0)); + + assert!((cell.geodesic_perimeter() - 443770.0).abs() < 1.0); + assert!((cell.geodesic_area() - 12308778361.0).abs() < 1.0); + + let center = Point::new(0.5, 0.5); + let outside = Point::new(1.5, 0.5); + + assert!(cell.contains(&Geometry::Point(center))); + assert!(!cell.contains(&Geometry::Point(outside))); + assert_eq!(center, cell.center()); + + let interpolated = cell.spatial_tec_interp(center).unwrap_or_else(|e| { + panic!("should have been feasible: {}", e); + }); + + assert_eq!(interpolated.tecu(), 1.0); + } + + #[test] + fn spatial_south_west_gradient_interpolation() { + let epoch = Epoch::default(); + + let northeast_tec = TEC::from_tecu(0.0); + let northwest_tec = TEC::from_tecu(0.0); + let southeast_tec = TEC::from_tecu(0.0); + let southwest_tec = TEC::from_tecu(1.0); + + let cell = MapCell::from_unitary_tec( + epoch, + northeast_tec, + northwest_tec, + southeast_tec, + southwest_tec, + ); + + for (x_deg, y_deg, tecu) in [ + (0.5, 0.5, 0.25), + (0.1, 0.1, 0.81), + (0.01, 0.01, 0.9801), + (0.0, 0.0, 1.0), + ] { + let point = Point::new(x_deg, y_deg); + + 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); + } + } + + #[test] + fn temporal_interpolation() { + let t0 = Epoch::default(); + let t1 = t0 + 1.0 * Unit::Second; + + let t_ok = t0 + 0.5 * Unit::Second; + let t_nok = t1 + 2.0 * Unit::Second; + + let center = Point::new(0.5, 0.5); + + let northeast_tec_0 = TEC::from_tecu(1.0); + let northwest_tec_0 = TEC::from_tecu(1.0); + let southeast_tec_0 = TEC::from_tecu(1.0); + let southwest_tec_0 = TEC::from_tecu(1.0); + + let cell0 = MapCell::from_unitary_tec( + t0, + northeast_tec_0, + northwest_tec_0, + southeast_tec_0, + southwest_tec_0, + ); + + let northeast_tec_1 = TEC::from_tecu(1.0); + let northwest_tec_1 = TEC::from_tecu(1.0); + let southeast_tec_1 = TEC::from_tecu(1.0); + let southwest_tec_1 = TEC::from_tecu(1.0); + + let cell1 = MapCell::from_unitary_tec( + t1, + northeast_tec_1, + northwest_tec_1, + southeast_tec_1, + southwest_tec_1, + ); + + assert!( + cell0 + .temporal_spatial_tec_interp(t_nok, center, &cell1) + .is_err(), + "interpolation is temporally incorrect!" + ); + + let tec = cell0 + .temporal_spatial_tec_interp(t_ok, center, &cell1) + .unwrap_or_else(|e| { + panic!("should have been feasible! {}", e); + }); + + assert_eq!(tec.tecu(), 1.0); + } +} 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/error.rs b/src/error.rs index 1f92308..ffe51bd 100644 --- a/src/error.rs +++ b/src/error.rs @@ -91,6 +91,30 @@ pub enum ParsingError { ExponentScaling, } +#[derive(Error, Debug)] +pub enum Error { + #[error("strech factor must be positive finite number")] + InvalidStretchFactor, + + #[error("undefined ROI exterior boundaries")] + UndefinedBoundaries, + + #[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, + + #[error("both regions are not synchronous in time")] + TemporalMismatch, + + #[error("invalid temporal interpolation instant")] + InvalidTemporalPoint, +} + /// Errors that may rise during Formatting process #[derive(Error, Debug)] pub enum FormattingError { diff --git a/src/production.rs b/src/file_attributes.rs similarity index 63% rename from src/production.rs rename to src/file_attributes.rs index 3848337..cca551a 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'), } } } @@ -30,15 +29,15 @@ 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)] -pub struct ProductionAttributes { - /// Production agency +#[derive(Debug, Clone, PartialEq)] +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,34 @@ pub struct ProductionAttributes { pub gzip_compressed: bool, } -impl std::fmt::Display for ProductionAttributes { +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 { + 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 { let len = std::cmp::min(self.agency.len(), 3); @@ -86,7 +112,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 +133,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 }; @@ -127,6 +153,33 @@ impl std::str::FromStr for ProductionAttributes { } } +#[cfg(feature = "qc")] +impl gnss_qc_traits::Merge for FileAttributes { + 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::Worldwide { + self.region = Region::Worldwide; + } + + #[cfg(feature = "flate2")] + if self.gzip_compressed && !rhs.gzip_compressed { + self.gzip_compressed = false; + } + + Ok(()) + } +} + #[cfg(test)] mod test { use super::*; @@ -135,15 +188,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); @@ -159,12 +212,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/grid.rs b/src/grid.rs index 1b72604..f24061c 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)] @@ -19,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/header/mod.rs b/src/header/mod.rs index 3d5215f..eaff34a 100644 --- a/src/header/mod.rs +++ b/src/header/mod.rs @@ -1,6 +1,9 @@ mod formatting; mod parsing; +#[cfg(feature = "qc")] +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..c755780 --- /dev/null +++ b/src/header/qc.rs @@ -0,0 +1,97 @@ +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() { + if !self.comments.contains(&comment) { + self.comments.push(comment.clone()); + } + } + + // 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 2336d34..5c83adc 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. @@ -23,12 +22,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; @@ -53,20 +52,21 @@ 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}; -use hifitime::prelude::Epoch; +use hifitime::prelude::{Epoch, TimeSeries}; use crate::{ - cell::{MapCell, MapPoint}, + cell::{Cell3x3, MapCell, TecPoint}, coordinates::QuantizedCoordinates, - error::{FormattingError, ParsingError}, + error::{Error, FormattingError, ParsingError}, + file_attributes::{FileAttributes, Region}, + grid::{Axis, Grid}, header::Header, key::Key, - production::{ProductionAttributes, Region}, quantized::Quantized, record::Record, tec::TEC, @@ -76,15 +76,15 @@ pub mod prelude { // export pub use crate::{ bias::BiasSource, - cell::MapCell, - error::{FormattingError, ParsingError}, - grid::Grid, + cell::{Cell3x3, MapCell}, + error::{Error, FormattingError, ParsingError}, + file_attributes::*, + grid::{Axis, Grid}, header::Header, ionosphere::IonosphereParameters, key::Key, linspace::Linspace, mapf::MappingFunction, - production::*, record::Record, system::ReferenceSystem, tec::TEC, @@ -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}; @@ -104,13 +104,49 @@ 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 + } +} + +/// Converts a geo [Rect]angle to NE, SE, SW, NW (latitude, longitude) quadruplets +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()); + let (x0, y0) = (min.x, min.y); + ( + (x0, y0), + (x0 + width, y0), + (x0 + width, y0 + height), + (x0, y0 + height), + ) +} + +/// 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 { format!("{: 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(); /// @@ -161,14 +197,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(); @@ -196,9 +232,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 { @@ -207,7 +242,7 @@ impl IONEX { Self { header, record, - production: None, + attributes: None, comments: Default::default(), } } @@ -218,7 +253,7 @@ impl IONEX { header, record: self.record.clone(), comments: self.comments.clone(), - production: self.production.clone(), + attributes: self.attributes.clone(), } } @@ -233,7 +268,7 @@ impl IONEX { record, header: self.header.clone(), comments: self.comments.clone(), - production: self.production.clone(), + attributes: self.attributes.clone(), } } @@ -259,21 +294,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 { "" @@ -285,12 +320,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; } @@ -299,16 +335,20 @@ 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 = if self.header.grid.is_worldwide() { + Region::Worldwide + } else { + Region::Regional + }; - 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 @@ -331,7 +371,7 @@ impl IONEX { header, record, comments, - production: Default::default(), + attributes: Default::default(), }) } @@ -358,12 +398,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(); @@ -371,7 +411,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 @@ -384,7 +424,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) } @@ -403,9 +443,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)?; @@ -453,7 +493,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 @@ -468,14 +508,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::*; @@ -486,10 +527,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> { @@ -512,22 +557,38 @@ impl IONEX { false } - /// Returns map borders as [Retc]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(); - 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 @@ -541,95 +602,335 @@ impl IONEX { ionex } - /// 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 { - let mut ionex = IONEX::default(); + /// 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. + /// + /// 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::Rect(boundaries); - 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); + let cells = self + .map_cell_iter() + .filter(|cell| cell.contains(&roi)) + .collect::>(); - // copy attributes - ionex.production = self.production.clone(); + ionex.record = Record::from_map_cells(&cells, self.header.grid.altitude.start); - if let Some(production) = &mut ionex.production { - production.region = Region::Regional; + // update attributes + match &self.attributes { + Some(attributes) => { + ionex.attributes = Some(attributes.clone().regionalized()); + }, + None => { + let mut attributes = FileAttributes::default().regionalized(); + ionex.attributes = Some(attributes); + }, } - // copy & rework header - ionex.header = self.header.clone(); + // 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.start = max_lat; - ionex.header.grid.latitude.end = min_lat; - ionex.header.grid.longitude.start = min_long; - ionex.header.grid.longitude.end = max_long; + 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); - let region = Geometry::Polygon(region); + Ok(ionex) + } - // restrict area - let cells = self - .map_cell_iter() - .filter(|cell| cell.contains(®ion)) - .collect::>(); + // /// 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, + // /// 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. + // /// + // /// 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; + // } + + // self.attributes = Some(attributes); + // }, + // } + + // // 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(()) + // } + // } + + /// 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()) + } - ionex.record = Record::from_map_cells( - self.header.grid.altitude.start, - min_lat, - max_lat, - min_long, - max_long, - &cells, - ); + // 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(); + + // // applies the stretching + + // // 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 + + // 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); + } - Some(ionex) - } + let new_dt = self.header.sampling_period * factor; - /// 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 { - let mut s = self.clone(); - s.grid_precision_stretch_mut(stretch_factor); - s - } + if factor > 1.0 { + // temporal upscaling + } else { + // temporal decimation + } - /// 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; + // update header + self.header.sampling_period = new_dt; - // // update map - // for epoch in self.record.epochs_iter() { - // } + Ok(()) } - /// Designs a [MapCell] iterator (micro rectangle region made of 4 local points) - /// that can then be interpolated. - pub fn map_cell_iter(&self) -> Box + '_> { - let timeseries = self.header.timeseries(); + // /// 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) + // } + + /// Produces a [TimeSeries] from this [IONEX], describing the temporal axis. + pub fn timeseries(&self) -> TimeSeries { + self.header.timeseries() + } + /// Designs a [MapCell] iterator (micro ROI following the grid quantization) + /// that allows micro interpolation. + pub fn map_cell_iter(&self) -> Box + '_> { 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; 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 { @@ -678,19 +979,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()), }, @@ -699,7 +1000,157 @@ impl IONEX { ) } - /// Returns the [MapCell] that contains following [Geometry]. + // /// 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, + /// 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); + } + + // 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; + + while t < self.header.epoch_of_last_map { + if t == epoch { + needs_temporal_interp = false; + break; + } + + t += self.header.sampling_period; + } + + // 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 (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, + ); + + 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 [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. + /// 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. + 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; + + while t < self.header.epoch_of_last_map { + if t == epoch { + needs_temporal_interp = false; + break; + } + + t += self.header.sampling_period; + } + + 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); + } + } + } + + None + } + + /// Obtain the best suited [MapCell] spatially wrapping this Geometry that contains following [Geometry]. /// /// ## Input /// - epoch: [Epoch] that must exist in this [IONEX] @@ -725,115 +1176,159 @@ 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.map_borders_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; + // /// 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 + // } +} - // for all requested coordinates - for point in area.points() { - let geo = Geometry::Point(point); +/// 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) + } - let key = Key::from_decimal_degrees_km(epoch, point.y(), point.x(), fixed_altitude_km); + fn merge_mut(&mut self, rhs: &Self) -> Result<(), gnss_qc_traits::MergeError> { + self.header.merge_mut(&rhs.header)?; + self.record.merge_mut(&rhs.record)?; - for cell in self.synchronous_map_cell_iter(epoch) { - if cell.contains(&geo) { - let tec = cell.spatial_interpolation(point); - values.insert(key, tec); + match self.attributes { + Some(ref mut prods) => { + if let Some(rhs) = &rhs.attributes { + prods.merge_mut(rhs)?; } - } + }, + None => { + if let Some(rhs) = &rhs.attributes { + self.attributes = Some(rhs.clone()); + } + }, } - 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); - } - } - } - } + // add new comments + for comment in rhs.comments.iter() { + if !self.comments.contains(&comment) { + self.comments.push(comment.clone()); } } - values + Ok(()) } } #[cfg(test)] mod test { - use crate::fmt_comment; - use crate::prelude::*; + use crate::{div_ceil, fmt_comment, prelude::*, rectangle_quadrant_decomposition}; #[test] fn fmt_comments_singleline() { @@ -861,7 +1356,7 @@ mod test { fn fmt_wrapped_comments() { for desc in ["just trying to form a very lengthy comment that will overflow since it does not fit in a single line", "just trying to form a very very lengthy comment that will overflow since it does fit on three very meaningful lines. Imazdmazdpoakzdpoakzpdokpokddddddddddddddddddaaaaaaaaaaaaaaaaaaaaaaa"] { - let nb_lines = num_integer::div_ceil(desc.len(), 60); + let nb_lines = div_ceil(desc.len(), 60); let comments = fmt_comment(desc); assert_eq!(comments.lines().count(), nb_lines); @@ -874,22 +1369,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 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)), + ( + (-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_quadrant_decomposition(rect), + ( + (lat11, long11), + (lat12, long12), + (lat21, long21), + (lat22, long22) + ), + "failed for {:?}", + rect + ); + } } } diff --git a/src/linspace.rs b/src/linspace.rs index 2a54218..ac81314 100644 --- a/src/linspace.rs +++ b/src/linspace.rs @@ -1,6 +1,9 @@ use std::ops::Rem; -use crate::{error::ParsingError, quantized::Quantized}; +use crate::{ + error::{Error, ParsingError}, + quantized::Quantized, +}; /// Quantized Linspace for iteration #[derive(Debug, Copy, Clone, Default, PartialEq, PartialOrd)] @@ -67,6 +70,51 @@ impl Linspace { (self.min(), self.max()) } + /// 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_normal() { + return Err(Error::InvalidStretchFactor); + } + self.start *= factor; + self.end *= factor; + Ok(()) + } + + /// 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 + /// 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)?; + Ok(s) + } + + /// Resample this mutable [Linspace] so the point spacing is modified. + /// 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_normal() { + return Err(Error::InvalidStretchFactor); + } + 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 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(); + 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 +237,74 @@ 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/record/mod.rs b/src/record/mod.rs index dc47da9..30dfcc8 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; @@ -56,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. @@ -70,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 } } @@ -131,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); @@ -155,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/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/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()) - } - } -} 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 { 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); } } diff --git a/src/tests/mod.rs b/src/tests/mod.rs index 7baf321..2ec56a5 100644 --- a/src/tests/mod.rs +++ b/src/tests/mod.rs @@ -3,6 +3,9 @@ pub mod toolkit; mod filename; mod parsing; +mod qc; +mod roi; +// mod stretching; mod v1; diff --git a/src/tests/qc.rs b/src/tests/qc.rs new file mode 100644 index 0000000..09522c7 --- /dev/null +++ b/src/tests/qc.rs @@ -0,0 +1,103 @@ +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); + }); + + assert!(merged.is_merged(), "merged file not declared as merged!"); + + // 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 mut parsed = IONEX::from_file("merged.txt").unwrap_or_else(|e| { + panic!("failed to parse back CKMG V1: {}", e); + }); + + 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); +} + +#[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); + }); + + // 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); + + 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); + }); + + // 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()); +} diff --git a/src/tests/roi.rs b/src/tests/roi.rs new file mode 100644 index 0000000..a7a2b24 --- /dev/null +++ b/src/tests/roi.rs @@ -0,0 +1,37 @@ +use crate::{ + prelude::{coord, Rect, IONEX}, + tests::init_logger, +}; + +#[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/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" + ); + } +} diff --git a/src/tests/v1.rs b/src/tests/v1.rs index da7c50f..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", @@ -220,7 +226,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)) ); @@ -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",