From 8fdbda2a23a9d176e2c332ac6fbaee1a946967ff Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Tue, 19 Mar 2024 17:42:31 +0100 Subject: [PATCH 01/10] A CoordTrait experiment --- src/coordinate/mod.rs | 237 +++++++++++++++++++++++++++++++++++++ src/ellipsoid/geodesics.rs | 39 +++--- 2 files changed, 261 insertions(+), 15 deletions(-) diff --git a/src/coordinate/mod.rs b/src/coordinate/mod.rs index 7ce29bdc..25b327c6 100644 --- a/src/coordinate/mod.rs +++ b/src/coordinate/mod.rs @@ -104,3 +104,240 @@ pub trait CoordinateSet: CoordinateMetadata { } } } + +// An experiment with an extended version of Kyle Barron's CoordTrait PR +// over at https://github.com/georust/geo/pull/1157 + +pub trait CoordNum {} +impl CoordNum for f32 {} +impl CoordNum for f64 {} + +/// A trait for accessing data from a generic Coord. +pub trait CoordTrait { + type T: CoordNum; + + /// Accessors for the coordinate tuple components + fn first(&self) -> Self::T; + fn second(&self) -> Self::T; + fn third(&self) -> Self::T; + fn fourth(&self) -> Self::T; + fn measure(&self) -> Self::T; + + /// Accessors for the coordinate tuple components + fn first_as_f64(&self) -> f64; + fn second_as_f64(&self) -> f64; + fn third_as_f64(&self) -> f64; + fn fourth_as_f64(&self) -> f64; + fn measure_as_f64(&self) -> f64; + + /// x component of this coord + fn x(&self) -> Self::T { + self.first() + } + + /// y component of this coord + fn y(&self) -> Self::T { + self.second() + } + + /// z component of this coord + fn z(&self) -> Self::T { + self.third() + } + + /// t component of this coord + fn t(&self) -> Self::T { + self.fourth() + } + + /// Returns a tuple that contains the two first components of the coord. + fn xy(&self) -> (Self::T, Self::T) { + (self.first(), self.second()) + } + + /// Returns a tuple that contains the three first components of the coord. + fn xyz(&self) -> (Self::T, Self::T, Self::T) { + (self.first(), self.second(), self.third()) + } + + /// Returns a tuple that contains the three first components of the coord. + fn xyzt(&self) -> (Self::T, Self::T, Self::T, Self::T) { + (self.first(), self.second(), self.third(), self.fourth()) + } + + /// Returns a tuple that contains the two first components of the coord converted to f64. + fn xy_as_f64(&self) -> (f64, f64) { + (self.first_as_f64(), self.second_as_f64()) + } + + /// Returns a tuple that contains the three first components of the coord converted to f64. + fn xyz_as_f64(&self) -> (f64, f64, f64) { + ( + self.first_as_f64(), + self.second_as_f64(), + self.third_as_f64(), + ) + } + + /// Returns a tuple that contains the three first components of the coord converted to f64. + fn xyzt_as_f64(&self) -> (f64, f64, f64, f64) { + ( + self.first_as_f64(), + self.second_as_f64(), + self.third_as_f64(), + self.fourth_as_f64(), + ) + } +} + +impl CoordTrait for Coor2D { + type T = f64; + + /// Accessors for the coordinate tuple components + fn first(&self) -> Self::T { + self.0[0] + } + fn second(&self) -> Self::T { + self.0[1] + } + fn third(&self) -> Self::T { + f64::NAN + } + fn fourth(&self) -> Self::T { + f64::NAN + } + fn measure(&self) -> Self::T { + f64::NAN + } + + /// Accessors for the coordinate tuple components + fn first_as_f64(&self) -> f64 { + self.0[0] + } + fn second_as_f64(&self) -> f64 { + self.0[1] + } + fn third_as_f64(&self) -> f64 { + f64::NAN + } + fn fourth_as_f64(&self) -> f64 { + f64::NAN + } + fn measure_as_f64(&self) -> f64 { + f64::NAN + } +} + +impl CoordTrait for Coor32 { + type T = f32; + + /// Accessors for the coordinate tuple components + fn first(&self) -> Self::T { + self.0[0] + } + fn second(&self) -> Self::T { + self.0[1] + } + fn third(&self) -> Self::T { + f32::NAN + } + fn fourth(&self) -> Self::T { + f32::NAN + } + fn measure(&self) -> Self::T { + f32::NAN + } + + /// Accessors for the coordinate tuple components + fn first_as_f64(&self) -> f64 { + self.0[0] as f64 + } + fn second_as_f64(&self) -> f64 { + self.0[1] as f64 + } + fn third_as_f64(&self) -> f64 { + f64::NAN + } + fn fourth_as_f64(&self) -> f64 { + f64::NAN + } + fn measure_as_f64(&self) -> f64 { + f64::NAN + } +} + +impl CoordTrait for Coor3D { + type T = f64; + + /// Accessors for the coordinate tuple components + fn first(&self) -> Self::T { + self.0[0] + } + fn second(&self) -> Self::T { + self.0[1] + } + fn third(&self) -> Self::T { + self.0[2] + } + fn fourth(&self) -> Self::T { + f64::NAN + } + fn measure(&self) -> Self::T { + f64::NAN + } + + /// Accessors for the coordinate tuple components + fn first_as_f64(&self) -> f64 { + self.0[0] + } + fn second_as_f64(&self) -> f64 { + self.0[1] + } + fn third_as_f64(&self) -> f64 { + self.0[2] + } + fn fourth_as_f64(&self) -> f64 { + f64::NAN + } + fn measure_as_f64(&self) -> f64 { + f64::NAN + } +} + +impl CoordTrait for Coor4D { + type T = f64; + + /// Accessors for the coordinate tuple components + fn first(&self) -> Self::T { + self.0[0] + } + fn second(&self) -> Self::T { + self.0[1] + } + fn third(&self) -> Self::T { + self.0[2] + } + fn fourth(&self) -> Self::T { + self.0[3] + } + fn measure(&self) -> Self::T { + f64::NAN + } + + /// Accessors for the coordinate tuple components + fn first_as_f64(&self) -> f64 { + self.0[0] + } + fn second_as_f64(&self) -> f64 { + self.0[1] + } + fn third_as_f64(&self) -> f64 { + self.0[2] + } + fn fourth_as_f64(&self) -> f64 { + self.0[3] + } + fn measure_as_f64(&self) -> f64 { + f64::NAN + } +} diff --git a/src/ellipsoid/geodesics.rs b/src/ellipsoid/geodesics.rs index 102bdf33..c298bb16 100644 --- a/src/ellipsoid/geodesics.rs +++ b/src/ellipsoid/geodesics.rs @@ -1,3 +1,7 @@ +use crate::coordinate::CoordTrait; + +// Now using an extended version of Kyle Barron's CoordTrait, cf. src/coordinate/mod.rs + use super::*; // ----- Geodesics ------------------------------------------------------------- @@ -14,10 +18,12 @@ impl Ellipsoid { /// Federico Dolce and Michael Kirk, provides a Rust implementation of Karney's algorithm. #[must_use] #[allow(non_snake_case)] - pub fn geodesic_fwd(&self, from: &Coor4D, azimuth: f64, distance: f64) -> Coor4D { + pub fn geodesic_fwd(&self, from: &T, azimuth: f64, distance: f64) -> Coor4D + where + T: CoordTrait, + { // Coordinates of the point of origin, P1 - let B1 = from[1]; - let L1 = from[0]; + let (L1, B1) = from.xy_as_f64(); // The latitude of P1 projected onto the auxiliary sphere let U1 = self.latitude_geographic_to_reduced(B1); @@ -93,13 +99,13 @@ impl Ellipsoid { /// See [`geodesic_fwd`](crate::Ellipsoid::geodesic_fwd) #[must_use] #[allow(non_snake_case)] // So we can use the mathematical notation from the original text - pub fn geodesic_inv(&self, from: &Coor4D, to: &Coor4D) -> Coor4D { - let B1 = from[1]; - let B2 = to[1]; + pub fn geodesic_inv(&self, from: &T, to: &T) -> Coor4D + where + T: CoordTrait, + { + let (L1, B1) = from.xy_as_f64(); + let (L2, B2) = to.xy_as_f64(); let B = B2 - B1; - - let L1 = from[0]; - let L2 = to[0]; let L = L2 - L1; // Below the micrometer level, we don't care about directions @@ -191,14 +197,17 @@ impl Ellipsoid { /// // Compute the distance between Copenhagen and Paris /// use geodesy::prelude::*; /// if let Ok(ellps) = Ellipsoid::named("GRS80") { - /// let p0 = Coor4D::geo(55., 12., 0., 0.); - /// let p1 = Coor4D::geo(49., 2., 0., 0.); + /// let p0 = Coor2D::geo(55., 12.); + /// let p1 = Coor2D::geo(49., 2.); /// let d = ellps.distance(&p0, &p1); /// assert!((d - 956_066.231_959).abs() < 1e-5); /// } /// ``` #[must_use] - pub fn distance(&self, from: &Coor4D, to: &Coor4D) -> f64 { + pub fn distance(&self, from: &T, to: &T) -> f64 + where + T: CoordTrait, + { self.geodesic_inv(from, to)[2] } } @@ -216,8 +225,8 @@ mod tests { // Copenhagen (Denmark)--Paris (France) // Expect distance good to 0.01 mm, azimuths to a nanodegree - let p1 = Coor4D::gis(12., 55., 0., 0.); - let p2 = Coor4D::gis(2., 49., 0., 0.); + let p1 = Coor2D::gis(12., 55.); + let p2 = Coor2D::gis(2., 49.); let d = ellps.geodesic_inv(&p1, &p2); assert!((d[0].to_degrees() - (-130.15406042072)).abs() < 1e-9); @@ -231,7 +240,7 @@ mod tests { // Copenhagen (Denmark)--Rabat (Morocco) // Expect distance good to 0.1 mm, azimuths to a nanodegree - let p2 = Coor4D::gis(7., 34., 0., 0.); + let p2 = Coor2D::gis(7., 34.); let d = ellps.geodesic_inv(&p1, &p2); assert!((d[0].to_degrees() - (-168.48914418666)).abs() < 1e-9); From 8800929f2e2a96553b1244ec14c5c4aa2477b77b Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Wed, 20 Mar 2024 14:00:01 +0100 Subject: [PATCH 02/10] Simplified, and added dimension+measure --- examples/01-geometric_geodesy.rs | 4 +- .../08-user_defined_operators_using_proj.rs | 5 +- src/coordinate/coor2d.rs | 12 +- src/coordinate/coor32.rs | 12 +- src/coordinate/coor3d.rs | 21 +-- src/coordinate/coor4d.rs | 17 +- src/coordinate/mod.rs | 165 ++++++++---------- src/ellipsoid/geodesics.rs | 5 +- src/inner_op/cart.rs | 3 +- src/inner_op/geodesic.rs | 6 +- src/inner_op/iso6709.rs | 7 +- src/inner_op/molodensky.rs | 11 +- 12 files changed, 111 insertions(+), 157 deletions(-) diff --git a/examples/01-geometric_geodesy.rs b/examples/01-geometric_geodesy.rs index 7b248ae1..30c57f53 100644 --- a/examples/01-geometric_geodesy.rs +++ b/examples/01-geometric_geodesy.rs @@ -46,8 +46,8 @@ fn main() -> Result<(), Error> { // surface of the ellipsoid. Let's compute the distance and // azimuth between the approximate locations of the airports // of Copenhagen (CPH) and Paris (CDG). - let CPH = Coor4D::geo(55., 12., 0., 0.); - let CDG = Coor4D::geo(49., 2., 0., 0.); + let CPH = Coor2D::geo(55., 12.); + let CDG = Coor2D::geo(49., 2.); // By historical convention the "from A to B" situation is considered // the inverse sense of the geodesic problem - hence `geodesic_inv`: diff --git a/examples/08-user_defined_operators_using_proj.rs b/examples/08-user_defined_operators_using_proj.rs index 63f9dfdb..1f4fe03b 100644 --- a/examples/08-user_defined_operators_using_proj.rs +++ b/examples/08-user_defined_operators_using_proj.rs @@ -173,6 +173,7 @@ pub fn proj_constructor(parameters: &RawParameters, _ctx: &dyn Context) -> Resul fn main() -> anyhow::Result<()> { let mut prv = geodesy::Minimal::new(); prv.register_op("proj", OpConstructor(proj_constructor)); + let e = Ellipsoid::default(); // Check that we can access the `proj` binary - if not, just ignore if Command::new("proj").stderr(Stdio::piped()).spawn().is_err() { @@ -195,7 +196,7 @@ fn main() -> anyhow::Result<()> { println!("projected: {:?}", geo[0]); ctx.apply(op, Inv, &mut geo)?; - assert!(rtp[0].default_ellps_dist(&geo[0]) < 1e-5); + assert!(e.distance(&rtp[0], &geo[0]) < 1e-5); println!("roundtrip: {:?}", geo[0].to_degrees()); // Inverted invocation - note "proj inv ..." @@ -208,7 +209,7 @@ fn main() -> anyhow::Result<()> { // Now, we get the inverse utm projection when calling the operator in the Fwd direction ctx.apply(op, Fwd, &mut utm)?; - assert!(geo[0].default_ellps_dist(&utm[0]) < 1e-5); + assert!(e.distance(&utm[0], &geo[0]) < 1e-5); // ...and roundtrip back to utm ctx.apply(op, Inv, &mut utm)?; assert!(rtp[0].hypot2(&utm[0]) < 1e-5); diff --git a/src/coordinate/coor2d.rs b/src/coordinate/coor2d.rs index 3f0ed8e3..0923c6dd 100644 --- a/src/coordinate/coor2d.rs +++ b/src/coordinate/coor2d.rs @@ -164,16 +164,9 @@ impl Coor2D { (self[0] - other[0]).hypot(self[1] - other[1]) } - /// The Geodesic distance on the default ellipsoid. Mostly a shortcut - /// for test authoring - pub fn default_ellps_dist(&self, other: &Self) -> f64 { - Ellipsoid::default().distance( - &Coor4D([self[0], self[1], 0., 0.]), - &Coor4D([other[0], other[1], 0., 0.]), - ) - } } + impl From for Coor4D { fn from(c: Coor2D) -> Self { Coor4D([c[0], c[1], 0.0, 0.0]) @@ -193,11 +186,12 @@ mod tests { use super::*; #[test] fn distances() { + let e = Ellipsoid::default(); let lat = angular::dms_to_dd(55, 30, 36.); let lon = angular::dms_to_dd(12, 45, 36.); let dms = Coor2D::geo(lat, lon); let geo = Coor2D::geo(55.51, 12.76); - assert!(geo.default_ellps_dist(&dms) < 1e-10); + assert!(e.distance(&geo, &dms) < 1e-10); } #[test] diff --git a/src/coordinate/coor32.rs b/src/coordinate/coor32.rs index cf2398c2..bc735667 100644 --- a/src/coordinate/coor32.rs +++ b/src/coordinate/coor32.rs @@ -177,15 +177,6 @@ impl Coor32 { pub fn hypot2(&self, other: &Self) -> f64 { (self[0] as f64 - other[0] as f64).hypot(self[1] as f64 - other[1] as f64) } - - /// The Geodesic distance on the default ellipsoid. Mostly a shortcut - /// for test authoring - pub fn default_ellps_dist(&self, other: &Self) -> f64 { - Ellipsoid::default().distance( - &Coor4D([self[0] as f64, self[1] as f64, 0., 0.]), - &Coor4D([other[0] as f64, other[1] as f64, 0., 0.]), - ) - } } // ----- T E S T S --------------------------------------------------- @@ -199,7 +190,8 @@ mod tests { let lon = angular::dms_to_dd(12, 45, 36.); let dms = Coor32::geo(lat, lon); let geo = Coor32::geo(55.51, 12.76); - assert!(geo.default_ellps_dist(&dms) < 1e-10); + let e = &Ellipsoid::default(); + assert!(e.distance(&dms, &geo) < 1e-9); } #[test] diff --git a/src/coordinate/coor3d.rs b/src/coordinate/coor3d.rs index 3cbc28ff..df7486b5 100644 --- a/src/coordinate/coor3d.rs +++ b/src/coordinate/coor3d.rs @@ -263,23 +263,6 @@ impl Coor3D { .hypot(self[2] - other[2]) } - /// The 3D distance between two points given as internal angular - /// coordinates. Mostly a shortcut for test authoring - pub fn default_ellps_3d_dist(&self, other: &Self) -> f64 { - let e = Ellipsoid::default(); - let from = Coor4D([self[0], self[1], self[2], 0.]); - let to = Coor4D([other[0], other[1], other[2], 0.]); - - e.cartesian(&from).hypot3(&e.cartesian(&to)) - } - - /// The Geodesic distance on the default ellipsoid. Mostly a shortcut - /// for test authoring - pub fn default_ellps_dist(&self, other: &Self) -> f64 { - let from = Coor4D([self[0], self[1], self[2], 0.]); - let to = Coor4D([other[0], other[1], other[2], 0.]); - Ellipsoid::default().distance(&from, &to) - } } // ----- T E S T S --------------------------------------------------- @@ -287,13 +270,15 @@ impl Coor3D { #[cfg(test)] mod tests { use super::*; + #[test] fn distances() { + let e = Ellipsoid::default(); let lat = angular::dms_to_dd(55, 30, 36.); let lon = angular::dms_to_dd(12, 45, 36.); let dms = Coor3D::geo(lat, lon, 0.); let geo = Coor3D::geo(55.51, 12.76, 0.); - assert!(geo.default_ellps_dist(&dms) < 1e-10); + assert!(e.distance(&geo, &dms) < 1e-10); } #[test] diff --git a/src/coordinate/coor4d.rs b/src/coordinate/coor4d.rs index f6969260..2661a6eb 100644 --- a/src/coordinate/coor4d.rs +++ b/src/coordinate/coor4d.rs @@ -269,19 +269,6 @@ impl Coor4D { .hypot(self[1] - other[1]) .hypot(self[2] - other[2]) } - - /// The 3D distance between two points given as internal angular - /// coordinates. Mostly a shortcut for test authoring - pub fn default_ellps_3d_dist(&self, other: &Self) -> f64 { - let e = Ellipsoid::default(); - e.cartesian(self).hypot3(&e.cartesian(other)) - } - - /// The Geodesic distance on the default ellipsoid. Mostly a shortcut - /// for test authoring - pub fn default_ellps_dist(&self, other: &Self) -> f64 { - Ellipsoid::default().distance(self, other) - } } // ----- T E S T S --------------------------------------------------- @@ -289,13 +276,15 @@ impl Coor4D { #[cfg(test)] mod tests { use super::*; + #[test] fn distances() { let lat = angular::dms_to_dd(55, 30, 36.); let lon = angular::dms_to_dd(12, 45, 36.); let dms = Coor4D::geo(lat, lon, 0., 2020.); let geo = Coor4D::geo(55.51, 12.76, 0., 2020.); - assert!(geo.default_ellps_dist(&dms) < 1e-10); + let e = Ellipsoid::default(); + assert!(e.distance(&geo, &dms) < 1e-10); } #[test] diff --git a/src/coordinate/mod.rs b/src/coordinate/mod.rs index 25b327c6..908c6d22 100644 --- a/src/coordinate/mod.rs +++ b/src/coordinate/mod.rs @@ -105,7 +105,7 @@ pub trait CoordinateSet: CoordinateMetadata { } } -// An experiment with an extended version of Kyle Barron's CoordTrait PR +// An experiment with an extended version of Kyle Barron's CoordTrait PR // over at https://github.com/georust/geo/pull/1157 pub trait CoordNum {} @@ -115,229 +115,216 @@ impl CoordNum for f64 {} /// A trait for accessing data from a generic Coord. pub trait CoordTrait { type T: CoordNum; + const DIMENSION: usize; + const MEASURE: bool; - /// Accessors for the coordinate tuple components - fn first(&self) -> Self::T; - fn second(&self) -> Self::T; - fn third(&self) -> Self::T; - fn fourth(&self) -> Self::T; - fn measure(&self) -> Self::T; + // Required implementations /// Accessors for the coordinate tuple components - fn first_as_f64(&self) -> f64; - fn second_as_f64(&self) -> f64; - fn third_as_f64(&self) -> f64; - fn fourth_as_f64(&self) -> f64; - fn measure_as_f64(&self) -> f64; + fn x(&self) -> Self::T; + fn y(&self) -> Self::T; + fn z(&self) -> Self::T; + fn t(&self) -> Self::T; + fn m(&self) -> Self::T; - /// x component of this coord - fn x(&self) -> Self::T { - self.first() - } + /// Accessors for the coordinate tuple components converted to f64 + fn x_as_f64(&self) -> f64; + fn y_as_f64(&self) -> f64; + fn z_as_f64(&self) -> f64; + fn t_as_f64(&self) -> f64; + fn m_as_f64(&self) -> f64; - /// y component of this coord - fn y(&self) -> Self::T { - self.second() - } - - /// z component of this coord - fn z(&self) -> Self::T { - self.third() - } - - /// t component of this coord - fn t(&self) -> Self::T { - self.fourth() - } + // Provided implementations /// Returns a tuple that contains the two first components of the coord. fn xy(&self) -> (Self::T, Self::T) { - (self.first(), self.second()) + (self.x(), self.y()) } /// Returns a tuple that contains the three first components of the coord. fn xyz(&self) -> (Self::T, Self::T, Self::T) { - (self.first(), self.second(), self.third()) + (self.x(), self.y(), self.z()) } /// Returns a tuple that contains the three first components of the coord. fn xyzt(&self) -> (Self::T, Self::T, Self::T, Self::T) { - (self.first(), self.second(), self.third(), self.fourth()) + (self.x(), self.y(), self.z(), self.t()) } /// Returns a tuple that contains the two first components of the coord converted to f64. fn xy_as_f64(&self) -> (f64, f64) { - (self.first_as_f64(), self.second_as_f64()) + (self.x_as_f64(), self.y_as_f64()) } /// Returns a tuple that contains the three first components of the coord converted to f64. fn xyz_as_f64(&self) -> (f64, f64, f64) { - ( - self.first_as_f64(), - self.second_as_f64(), - self.third_as_f64(), - ) + (self.x_as_f64(), self.y_as_f64(), self.z_as_f64()) } /// Returns a tuple that contains the three first components of the coord converted to f64. fn xyzt_as_f64(&self) -> (f64, f64, f64, f64) { - ( - self.first_as_f64(), - self.second_as_f64(), - self.third_as_f64(), - self.fourth_as_f64(), - ) + (self.x_as_f64(), self.y_as_f64(), self.z_as_f64(), self.t_as_f64()) } } impl CoordTrait for Coor2D { type T = f64; + const DIMENSION: usize = 2; + const MEASURE: bool = false; /// Accessors for the coordinate tuple components - fn first(&self) -> Self::T { + fn x(&self) -> Self::T { self.0[0] } - fn second(&self) -> Self::T { + fn y(&self) -> Self::T { self.0[1] } - fn third(&self) -> Self::T { + fn z(&self) -> Self::T { f64::NAN } - fn fourth(&self) -> Self::T { + fn t(&self) -> Self::T { f64::NAN } - fn measure(&self) -> Self::T { + fn m(&self) -> Self::T { f64::NAN } - /// Accessors for the coordinate tuple components - fn first_as_f64(&self) -> f64 { + /// Accessors for the coordinate tuple components converted to f64 + fn x_as_f64(&self) -> f64 { self.0[0] } - fn second_as_f64(&self) -> f64 { + fn y_as_f64(&self) -> f64 { self.0[1] } - fn third_as_f64(&self) -> f64 { + fn z_as_f64(&self) -> f64 { f64::NAN } - fn fourth_as_f64(&self) -> f64 { + fn t_as_f64(&self) -> f64 { f64::NAN } - fn measure_as_f64(&self) -> f64 { + fn m_as_f64(&self) -> f64 { f64::NAN } } + + impl CoordTrait for Coor32 { type T = f32; + const DIMENSION: usize = 2; + const MEASURE: bool = false; /// Accessors for the coordinate tuple components - fn first(&self) -> Self::T { + fn x(&self) -> Self::T { self.0[0] } - fn second(&self) -> Self::T { + fn y(&self) -> Self::T { self.0[1] } - fn third(&self) -> Self::T { + fn z(&self) -> Self::T { f32::NAN } - fn fourth(&self) -> Self::T { + fn t(&self) -> Self::T { f32::NAN } - fn measure(&self) -> Self::T { + fn m(&self) -> Self::T { f32::NAN } - /// Accessors for the coordinate tuple components - fn first_as_f64(&self) -> f64 { + /// Accessors for the coordinate tuple components converted to f64 + fn x_as_f64(&self) -> f64 { self.0[0] as f64 } - fn second_as_f64(&self) -> f64 { + fn y_as_f64(&self) -> f64 { self.0[1] as f64 } - fn third_as_f64(&self) -> f64 { + fn z_as_f64(&self) -> f64 { f64::NAN } - fn fourth_as_f64(&self) -> f64 { + fn t_as_f64(&self) -> f64 { f64::NAN } - fn measure_as_f64(&self) -> f64 { + fn m_as_f64(&self) -> f64 { f64::NAN } } impl CoordTrait for Coor3D { type T = f64; + const DIMENSION: usize = 3; + const MEASURE: bool = false; /// Accessors for the coordinate tuple components - fn first(&self) -> Self::T { + fn x(&self) -> Self::T { self.0[0] } - fn second(&self) -> Self::T { + fn y(&self) -> Self::T { self.0[1] } - fn third(&self) -> Self::T { + fn z(&self) -> Self::T { self.0[2] } - fn fourth(&self) -> Self::T { + fn t(&self) -> Self::T { f64::NAN } - fn measure(&self) -> Self::T { + fn m(&self) -> Self::T { f64::NAN } - /// Accessors for the coordinate tuple components - fn first_as_f64(&self) -> f64 { + /// Accessors for the coordinate tuple components converted to f64 + fn x_as_f64(&self) -> f64 { self.0[0] } - fn second_as_f64(&self) -> f64 { + fn y_as_f64(&self) -> f64 { self.0[1] } - fn third_as_f64(&self) -> f64 { + fn z_as_f64(&self) -> f64 { self.0[2] } - fn fourth_as_f64(&self) -> f64 { + fn t_as_f64(&self) -> f64 { f64::NAN } - fn measure_as_f64(&self) -> f64 { + fn m_as_f64(&self) -> f64 { f64::NAN } } impl CoordTrait for Coor4D { type T = f64; + const DIMENSION: usize = 4; + const MEASURE: bool = false; /// Accessors for the coordinate tuple components - fn first(&self) -> Self::T { + fn x(&self) -> Self::T { self.0[0] } - fn second(&self) -> Self::T { + fn y(&self) -> Self::T { self.0[1] } - fn third(&self) -> Self::T { + fn z(&self) -> Self::T { self.0[2] } - fn fourth(&self) -> Self::T { + fn t(&self) -> Self::T { self.0[3] } - fn measure(&self) -> Self::T { + fn m(&self) -> Self::T { f64::NAN } - /// Accessors for the coordinate tuple components - fn first_as_f64(&self) -> f64 { + /// Accessors for the coordinate tuple components converted to f64 + fn x_as_f64(&self) -> f64 { self.0[0] } - fn second_as_f64(&self) -> f64 { + fn y_as_f64(&self) -> f64 { self.0[1] } - fn third_as_f64(&self) -> f64 { + fn z_as_f64(&self) -> f64 { self.0[2] } - fn fourth_as_f64(&self) -> f64 { + fn t_as_f64(&self) -> f64 { self.0[3] } - fn measure_as_f64(&self) -> f64 { + fn m_as_f64(&self) -> f64 { f64::NAN } } diff --git a/src/ellipsoid/geodesics.rs b/src/ellipsoid/geodesics.rs index c298bb16..33048ed1 100644 --- a/src/ellipsoid/geodesics.rs +++ b/src/ellipsoid/geodesics.rs @@ -206,9 +206,9 @@ impl Ellipsoid { #[must_use] pub fn distance(&self, from: &T, to: &T) -> f64 where - T: CoordTrait, + T: CoordTrait { - self.geodesic_inv(from, to)[2] + self.geodesic_inv::(&from, &to)[2] } } @@ -217,6 +217,7 @@ impl Ellipsoid { #[cfg(test)] mod tests { use super::*; + use crate::Coor2D; #[test] fn geodesics() -> Result<(), Error> { let ellps = Ellipsoid::named("GRS80")?; diff --git a/src/inner_op/cart.rs b/src/inner_op/cart.rs index 9f33213b..9206a4b5 100644 --- a/src/inner_op/cart.rs +++ b/src/inner_op/cart.rs @@ -169,6 +169,7 @@ mod tests { ), ]; + let e = Ellipsoid::default(); // Forward let mut operands = geo; ctx.apply(op, Fwd, &mut operands)?; @@ -179,7 +180,7 @@ mod tests { // Inverse ctx.apply(op, Inv, &mut operands)?; for i in 0..5 { - assert!(operands[i].default_ellps_3d_dist(&geo[i]) < 10e-9); + assert!(e.distance(&operands[i], &geo[i]) < 1e-8); } Ok(()) diff --git a/src/inner_op/geodesic.rs b/src/inner_op/geodesic.rs index 3e9f50c5..2b250ffe 100644 --- a/src/inner_op/geodesic.rs +++ b/src/inner_op/geodesic.rs @@ -12,7 +12,7 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { let mut successes = 0_usize; for i in sliced { let args = operands.get_coord(i); - let origin = Coor4D::geo(args[0], args[1], 0.0, 0.0); + let origin = Coor2D::geo(args[0], args[1]); let azimuth = args[2].to_radians(); let distance = args[3]; @@ -43,8 +43,8 @@ fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { let mut successes = 0_usize; for i in sliced { - let mut from = Coor4D::origin(); - let mut to = Coor4D::origin(); + let mut from = Coor2D::origin(); + let mut to = Coor2D::origin(); let coord = operands.get_coord(i); from[0] = coord[1].to_radians(); diff --git a/src/inner_op/iso6709.rs b/src/inner_op/iso6709.rs index ffb12c54..d542a549 100644 --- a/src/inner_op/iso6709.rs +++ b/src/inner_op/iso6709.rs @@ -144,11 +144,12 @@ mod tests { let mut ctx = Minimal::default(); let op = ctx.op("dms")?; - let mut operands = [Coor4D::raw(553036., -124509., 0., 0.)]; - let geo = Coor4D::geo(55.51, -12.7525, 0., 0.); + let mut operands = [Coor2D::raw(553036., -124509.)]; + let geo = Coor2D::geo(55.51, -12.7525); ctx.apply(op, Fwd, &mut operands)?; - assert!(operands[0].default_ellps_dist(&geo) < 1e-10); + let e = Ellipsoid::default(); + assert!(e.distance(&operands[0], &geo) < 1e-10); ctx.apply(op, Inv, &mut operands)?; assert!((operands[0][0] - 553036.).abs() < 1e-10); diff --git a/src/inner_op/molodensky.rs b/src/inner_op/molodensky.rs index f926193e..978b8230 100644 --- a/src/inner_op/molodensky.rs +++ b/src/inner_op/molodensky.rs @@ -209,9 +209,12 @@ mod tests { use super::*; use crate::math::angular; + #[test] fn molodensky() -> Result<(), Error> { let mut ctx = Minimal::default(); + let e = Ellipsoid::default(); + // --------------------------------------------------------------------------- // Test case from OGP Publication 373-7-2: Geomatics Guidance Note number 7, // part 2: Transformation from WGS84 to ED50. @@ -244,14 +247,14 @@ mod tests { // within 5 mm in the plane and the elevation. let mut operands = [WGS84]; ctx.apply(op, Fwd, &mut operands)?; - assert!(ED50.default_ellps_dist(&operands[0]) < 0.005); + assert!(e.distance(&ED50, &operands[0]) < 0.005); assert!((ED50[2] - operands[0][2]).abs() < 0.005); // The same holds in the reverse unabridged case, where // additionally the elevation is even better let mut operands = [ED50]; ctx.apply(op, Inv, &mut operands)?; - assert!(WGS84.default_ellps_3d_dist(&operands[0]) < 0.005); + assert!(e.distance(&WGS84, &operands[0]) < 0.005); assert!((WGS84[2] - operands[0][2]).abs() < 0.001); // The abridged case. Same test point. Both plane coordinates and @@ -264,13 +267,13 @@ mod tests { let mut operands = [WGS84]; ctx.apply(op, Fwd, &mut operands)?; - assert!(ED50.default_ellps_dist(&operands[0]) < 0.1); + assert!(e.distance(&ED50, &operands[0]) < 0.1); // Heights are worse in the abridged case assert!((ED50[2] - operands[0][2]).abs() < 0.075); let mut operands = [ED50]; ctx.apply(op, Inv, &mut operands)?; - assert!(WGS84.default_ellps_dist(&operands[0]) < 0.1); + assert!(e.distance(&WGS84, &operands[0]) < 0.1); // Heights are worse in the abridged case assert!((WGS84[2] - operands[0][2]).abs() < 0.075); Ok(()) From 1e27657719ab2b92b00a80f7e9dcc4829d4bb445 Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Wed, 20 Mar 2024 14:47:00 +0100 Subject: [PATCH 03/10] Cleanup --- src/coordinate/coor2d.rs | 2 -- src/coordinate/coor3d.rs | 1 - src/coordinate/mod.rs | 9 ++++++--- src/ellipsoid/geodesics.rs | 4 ++-- src/inner_op/molodensky.rs | 1 - 5 files changed, 8 insertions(+), 9 deletions(-) diff --git a/src/coordinate/coor2d.rs b/src/coordinate/coor2d.rs index 0923c6dd..255b5039 100644 --- a/src/coordinate/coor2d.rs +++ b/src/coordinate/coor2d.rs @@ -163,10 +163,8 @@ impl Coor2D { pub fn hypot2(&self, other: &Self) -> f64 { (self[0] - other[0]).hypot(self[1] - other[1]) } - } - impl From for Coor4D { fn from(c: Coor2D) -> Self { Coor4D([c[0], c[1], 0.0, 0.0]) diff --git a/src/coordinate/coor3d.rs b/src/coordinate/coor3d.rs index df7486b5..f18feca2 100644 --- a/src/coordinate/coor3d.rs +++ b/src/coordinate/coor3d.rs @@ -262,7 +262,6 @@ impl Coor3D { .hypot(self[1] - other[1]) .hypot(self[2] - other[2]) } - } // ----- T E S T S --------------------------------------------------- diff --git a/src/coordinate/mod.rs b/src/coordinate/mod.rs index 908c6d22..92d7217c 100644 --- a/src/coordinate/mod.rs +++ b/src/coordinate/mod.rs @@ -163,7 +163,12 @@ pub trait CoordTrait { /// Returns a tuple that contains the three first components of the coord converted to f64. fn xyzt_as_f64(&self) -> (f64, f64, f64, f64) { - (self.x_as_f64(), self.y_as_f64(), self.z_as_f64(), self.t_as_f64()) + ( + self.x_as_f64(), + self.y_as_f64(), + self.z_as_f64(), + self.t_as_f64(), + ) } } @@ -207,8 +212,6 @@ impl CoordTrait for Coor2D { } } - - impl CoordTrait for Coor32 { type T = f32; const DIMENSION: usize = 2; diff --git a/src/ellipsoid/geodesics.rs b/src/ellipsoid/geodesics.rs index 33048ed1..895f5188 100644 --- a/src/ellipsoid/geodesics.rs +++ b/src/ellipsoid/geodesics.rs @@ -206,9 +206,9 @@ impl Ellipsoid { #[must_use] pub fn distance(&self, from: &T, to: &T) -> f64 where - T: CoordTrait + T: CoordTrait, { - self.geodesic_inv::(&from, &to)[2] + self.geodesic_inv::(from, to)[2] } } diff --git a/src/inner_op/molodensky.rs b/src/inner_op/molodensky.rs index 978b8230..92c85b22 100644 --- a/src/inner_op/molodensky.rs +++ b/src/inner_op/molodensky.rs @@ -209,7 +209,6 @@ mod tests { use super::*; use crate::math::angular; - #[test] fn molodensky() -> Result<(), Error> { let mut ctx = Minimal::default(); From 80120f55c6e2d23f7625469017cdf548532338d5 Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Thu, 21 Mar 2024 00:57:12 +0100 Subject: [PATCH 04/10] Much cleanup, num_traits new dependency --- Cargo.toml | 1 + src/coordinate/mod.rs | 243 +++++++++++-------------------------- src/ellipsoid/geodesics.rs | 20 ++- 3 files changed, 81 insertions(+), 183 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 3bb8e578..65d1cbe8 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -18,6 +18,7 @@ default-run = "kp" [dependencies] # Library functionality uuid = { version = "1", features = ["v4"] } +num-traits = { version = "0.2" } # Command line program helpers clap = { version = "4.3.11", features = ["derive"], optional = true } diff --git a/src/coordinate/mod.rs b/src/coordinate/mod.rs index 92d7217c..5bbe37f3 100644 --- a/src/coordinate/mod.rs +++ b/src/coordinate/mod.rs @@ -105,21 +105,29 @@ pub trait CoordinateSet: CoordinateMetadata { } } -// An experiment with an extended version of Kyle Barron's CoordTrait PR +//----------------------------------------------------------------------- +// An experiment with an extended version of Kyle Barron's CoordTrait PR // over at https://github.com/georust/geo/pull/1157 - -pub trait CoordNum {} -impl CoordNum for f32 {} -impl CoordNum for f64 {} - -/// A trait for accessing data from a generic Coord. +//----------------------------------------------------------------------- + +// The next 8 lines are mostly copied from georust/geo/geo-types/src/lib.rs +// Although I added ToPrimitive in the first line +use core::fmt::Debug; +use num_traits::{Float, Num, NumCast, ToPrimitive}; +pub trait CoordinateType: Num + Copy + NumCast + PartialOrd + Debug {} +impl CoordinateType for T {} +pub trait CoordNum: CoordinateType + Debug {} +impl CoordNum for T {} +pub trait CoordFloat: CoordNum + Float {} +impl CoordFloat for T {} + +// And here is Kyle Barron's CoordTrait from https://github.com/georust/geo/pull/1157 +// extended with z(), t(), m(), and associated consts DIMENSION and MEASURE pub trait CoordTrait { type T: CoordNum; const DIMENSION: usize; const MEASURE: bool; - // Required implementations - /// Accessors for the coordinate tuple components fn x(&self) -> Self::T; fn y(&self) -> Self::T; @@ -127,207 +135,100 @@ pub trait CoordTrait { fn t(&self) -> Self::T; fn m(&self) -> Self::T; - /// Accessors for the coordinate tuple components converted to f64 - fn x_as_f64(&self) -> f64; - fn y_as_f64(&self) -> f64; - fn z_as_f64(&self) -> f64; - fn t_as_f64(&self) -> f64; - fn m_as_f64(&self) -> f64; - - // Provided implementations - /// Returns a tuple that contains the two first components of the coord. - fn xy(&self) -> (Self::T, Self::T) { + fn x_y(&self) -> (Self::T, Self::T) { (self.x(), self.y()) } +} - /// Returns a tuple that contains the three first components of the coord. - fn xyz(&self) -> (Self::T, Self::T, Self::T) { - (self.x(), self.y(), self.z()) - } +// The CoordTuples trait is blanket-implemented for anything that +// CoordTrait is implemented for +impl CoordTuples for C {} - /// Returns a tuple that contains the three first components of the coord. - fn xyzt(&self) -> (Self::T, Self::T, Self::T, Self::T) { - (self.x(), self.y(), self.z(), self.t()) - } +// And here the actual implementation, which takes any CoordTrait implementing +// data type, and lets us access the contents as geodesy-compatible f64 tuples +#[rustfmt::skip] +pub trait CoordTuples: CoordTrait { + /// Accessors for the coordinate tuple components converted to f64 + fn x_as_f64(&self) -> f64 { self.x().to_f64().unwrap_or(f64::NAN) } + fn y_as_f64(&self) -> f64 { self.y().to_f64().unwrap_or(f64::NAN) } + fn z_as_f64(&self) -> f64 { self.z().to_f64().unwrap_or(f64::NAN) } + fn t_as_f64(&self) -> f64 { self.t().to_f64().unwrap_or(f64::NAN) } + fn m_as_f64(&self) -> f64 { self.m().to_f64().unwrap_or(f64::NAN) } - /// Returns a tuple that contains the two first components of the coord converted to f64. fn xy_as_f64(&self) -> (f64, f64) { (self.x_as_f64(), self.y_as_f64()) } - /// Returns a tuple that contains the three first components of the coord converted to f64. fn xyz_as_f64(&self) -> (f64, f64, f64) { (self.x_as_f64(), self.y_as_f64(), self.z_as_f64()) } - /// Returns a tuple that contains the three first components of the coord converted to f64. fn xyzt_as_f64(&self) -> (f64, f64, f64, f64) { - ( - self.x_as_f64(), - self.y_as_f64(), - self.z_as_f64(), - self.t_as_f64(), - ) + (self.x_as_f64(), self.y_as_f64(), self.z_as_f64(), self.t_as_f64()) } } +// We must still implement the foundational CoordTrait trait for +// the Geodesy data types Coor2D, Coor32, Coor3D, Coor4D + +#[rustfmt::skip] impl CoordTrait for Coor2D { type T = f64; const DIMENSION: usize = 2; const MEASURE: bool = false; - - /// Accessors for the coordinate tuple components - fn x(&self) -> Self::T { - self.0[0] - } - fn y(&self) -> Self::T { - self.0[1] - } - fn z(&self) -> Self::T { - f64::NAN - } - fn t(&self) -> Self::T { - f64::NAN - } - fn m(&self) -> Self::T { - f64::NAN - } - - /// Accessors for the coordinate tuple components converted to f64 - fn x_as_f64(&self) -> f64 { - self.0[0] - } - fn y_as_f64(&self) -> f64 { - self.0[1] - } - fn z_as_f64(&self) -> f64 { - f64::NAN - } - fn t_as_f64(&self) -> f64 { - f64::NAN - } - fn m_as_f64(&self) -> f64 { - f64::NAN - } + fn x(&self) -> Self::T { self.0[0] } + fn y(&self) -> Self::T { self.0[1] } + fn z(&self) -> Self::T { f64::NAN } + fn t(&self) -> Self::T { f64::NAN } + fn m(&self) -> Self::T { f64::NAN } } +#[rustfmt::skip] impl CoordTrait for Coor32 { type T = f32; const DIMENSION: usize = 2; const MEASURE: bool = false; - - /// Accessors for the coordinate tuple components - fn x(&self) -> Self::T { - self.0[0] - } - fn y(&self) -> Self::T { - self.0[1] - } - fn z(&self) -> Self::T { - f32::NAN - } - fn t(&self) -> Self::T { - f32::NAN - } - fn m(&self) -> Self::T { - f32::NAN - } - - /// Accessors for the coordinate tuple components converted to f64 - fn x_as_f64(&self) -> f64 { - self.0[0] as f64 - } - fn y_as_f64(&self) -> f64 { - self.0[1] as f64 - } - fn z_as_f64(&self) -> f64 { - f64::NAN - } - fn t_as_f64(&self) -> f64 { - f64::NAN - } - fn m_as_f64(&self) -> f64 { - f64::NAN - } + fn x(&self) -> Self::T { self.0[0] } + fn y(&self) -> Self::T { self.0[1] } + fn z(&self) -> Self::T { f32::NAN } + fn t(&self) -> Self::T { f32::NAN } + fn m(&self) -> Self::T { f32::NAN } } +#[rustfmt::skip] impl CoordTrait for Coor3D { type T = f64; const DIMENSION: usize = 3; const MEASURE: bool = false; - - /// Accessors for the coordinate tuple components - fn x(&self) -> Self::T { - self.0[0] - } - fn y(&self) -> Self::T { - self.0[1] - } - fn z(&self) -> Self::T { - self.0[2] - } - fn t(&self) -> Self::T { - f64::NAN - } - fn m(&self) -> Self::T { - f64::NAN - } - - /// Accessors for the coordinate tuple components converted to f64 - fn x_as_f64(&self) -> f64 { - self.0[0] - } - fn y_as_f64(&self) -> f64 { - self.0[1] - } - fn z_as_f64(&self) -> f64 { - self.0[2] - } - fn t_as_f64(&self) -> f64 { - f64::NAN - } - fn m_as_f64(&self) -> f64 { - f64::NAN - } + fn x(&self) -> Self::T { self.0[0] } + fn y(&self) -> Self::T { self.0[1] } + fn z(&self) -> Self::T { self.0[2] } + fn t(&self) -> Self::T { f64::NAN } + fn m(&self) -> Self::T { f64::NAN } } +#[rustfmt::skip] impl CoordTrait for Coor4D { type T = f64; const DIMENSION: usize = 4; const MEASURE: bool = false; + fn x(&self) -> Self::T { self.0[0] } + fn y(&self) -> Self::T { self.0[1] } + fn z(&self) -> Self::T { self.0[2] } + fn t(&self) -> Self::T { self.0[3] } + fn m(&self) -> Self::T { f64::NAN } +} - /// Accessors for the coordinate tuple components - fn x(&self) -> Self::T { - self.0[0] - } - fn y(&self) -> Self::T { - self.0[1] - } - fn z(&self) -> Self::T { - self.0[2] - } - fn t(&self) -> Self::T { - self.0[3] - } - fn m(&self) -> Self::T { - f64::NAN - } - - /// Accessors for the coordinate tuple components converted to f64 - fn x_as_f64(&self) -> f64 { - self.0[0] - } - fn y_as_f64(&self) -> f64 { - self.0[1] - } - fn z_as_f64(&self) -> f64 { - self.0[2] - } - fn t_as_f64(&self) -> f64 { - self.0[3] - } - fn m_as_f64(&self) -> f64 { - f64::NAN - } +// And let's also implement it for a plain 2D f64 tuple +#[rustfmt::skip] +impl CoordTrait for (f64, f64) { + type T = f64; + const DIMENSION: usize = 2; + const MEASURE: bool = false; + fn x(&self) -> Self::T { self.0 } + fn y(&self) -> Self::T { self.1 } + fn z(&self) -> Self::T { f64::NAN } + fn t(&self) -> Self::T { f64::NAN } + fn m(&self) -> Self::T { f64::NAN } } diff --git a/src/ellipsoid/geodesics.rs b/src/ellipsoid/geodesics.rs index 895f5188..91d1378d 100644 --- a/src/ellipsoid/geodesics.rs +++ b/src/ellipsoid/geodesics.rs @@ -1,4 +1,6 @@ -use crate::coordinate::CoordTrait; +use num_traits::AsPrimitive; + +use crate::coordinate::{CoordTrait, CoordTuples}; // Now using an extended version of Kyle Barron's CoordTrait, cf. src/coordinate/mod.rs @@ -18,10 +20,7 @@ impl Ellipsoid { /// Federico Dolce and Michael Kirk, provides a Rust implementation of Karney's algorithm. #[must_use] #[allow(non_snake_case)] - pub fn geodesic_fwd(&self, from: &T, azimuth: f64, distance: f64) -> Coor4D - where - T: CoordTrait, - { + pub fn geodesic_fwd(&self, from: &G, azimuth: f64, distance: f64) -> Coor4D { // Coordinates of the point of origin, P1 let (L1, B1) = from.xy_as_f64(); @@ -99,10 +98,7 @@ impl Ellipsoid { /// See [`geodesic_fwd`](crate::Ellipsoid::geodesic_fwd) #[must_use] #[allow(non_snake_case)] // So we can use the mathematical notation from the original text - pub fn geodesic_inv(&self, from: &T, to: &T) -> Coor4D - where - T: CoordTrait, - { + pub fn geodesic_inv(&self, from: &G, to: &G) -> Coor4D { let (L1, B1) = from.xy_as_f64(); let (L2, B2) = to.xy_as_f64(); let B = B2 - B1; @@ -204,11 +200,11 @@ impl Ellipsoid { /// } /// ``` #[must_use] - pub fn distance(&self, from: &T, to: &T) -> f64 + pub fn distance(&self, from: &G, to: &G) -> f64 where - T: CoordTrait, + G::T: AsPrimitive, { - self.geodesic_inv::(from, to)[2] + self.geodesic_inv::(from, to)[2] } } From 09c2140c265a0ac2ee423036255cf04dfa68757b Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Tue, 2 Apr 2024 22:12:23 +0200 Subject: [PATCH 05/10] Trait impl of ISO-19111 coordinate tuple --- src/coordinate/mod.rs | 203 +++++++++++++++++++------------------ src/ellipsoid/geodesics.rs | 28 +++-- src/inner_op/units.rs | 4 + src/lib.rs | 2 + 4 files changed, 124 insertions(+), 113 deletions(-) diff --git a/src/coordinate/mod.rs b/src/coordinate/mod.rs index 5bbe37f3..caba5153 100644 --- a/src/coordinate/mod.rs +++ b/src/coordinate/mod.rs @@ -82,6 +82,8 @@ pub trait CoordinateMetadata { impl CoordinateMetadata for T where T: ?Sized {} /// CoordinateSet is the fundamental coordinate access interface in ISO-19111. +/// Strictly speaking, it is not a set, but (in abstract terms) rather an +/// indexed list, or (in more concrete terms): An array. /// /// Here it is implemented simply as an accessor trait, that allows us to /// access any user provided data model by iterating over its elements, @@ -89,13 +91,18 @@ impl CoordinateMetadata for T where T: ?Sized {} pub trait CoordinateSet: CoordinateMetadata { /// Number of coordinate tuples in the set fn len(&self) -> usize; + /// Access the `index`th coordinate tuple fn get_coord(&self, index: usize) -> Coor4D; + /// Overwrite the `index`th coordinate tuple fn set_coord(&mut self, index: usize, value: &Coor4D); + + /// Companion to `len()` fn is_empty(&self) -> bool { self.len() == 0 } + /// Set all coordinate tuples in the set to NaN fn stomp(&mut self) { let nanny = Coor4D::nan(); @@ -105,130 +112,130 @@ pub trait CoordinateSet: CoordinateMetadata { } } -//----------------------------------------------------------------------- -// An experiment with an extended version of Kyle Barron's CoordTrait PR -// over at https://github.com/georust/geo/pull/1157 -//----------------------------------------------------------------------- - -// The next 8 lines are mostly copied from georust/geo/geo-types/src/lib.rs -// Although I added ToPrimitive in the first line -use core::fmt::Debug; -use num_traits::{Float, Num, NumCast, ToPrimitive}; -pub trait CoordinateType: Num + Copy + NumCast + PartialOrd + Debug {} -impl CoordinateType for T {} -pub trait CoordNum: CoordinateType + Debug {} -impl CoordNum for T {} -pub trait CoordFloat: CoordNum + Float {} -impl CoordFloat for T {} - -// And here is Kyle Barron's CoordTrait from https://github.com/georust/geo/pull/1157 -// extended with z(), t(), m(), and associated consts DIMENSION and MEASURE -pub trait CoordTrait { - type T: CoordNum; +/// The CoordinateTuple is the ISO-19111 atomic spatial/spatiotemporal +/// referencing element. Loosely speaking, a CoordinateSet consists of +/// CoordinateTuples. +/// +/// Note that (despite the formal name) the underlying data structure +/// need not be a tuple: It can be any item, for which it makes sense +/// to implement the CoordinateTuple trait. +/// +/// The CoordinateTuple trait provides a number of convenience accessors +/// for accessing single coordinate elements or tuples of subsets. +/// These accessors are pragmatically named (x, y, xy, etc.). While these +/// names may be geodetically naive, they are suggestive and practical, and +/// aligns well with the internal coordinate order convention of most +/// Geodesy operators. +/// +/// All accessors have default implementations, except the +/// [`unchecked_nth()`](crate::coordinate::CoordinateTuple::unchecked_nth) function, +/// which must be provided by the implementer. +/// +/// When accessing dimensions outside of the domain of the CoordinateTuple, +/// [NaN](f64::NAN) will be returned. +#[rustfmt::skip] +pub trait CoordinateTuple { const DIMENSION: usize; - const MEASURE: bool; - /// Accessors for the coordinate tuple components - fn x(&self) -> Self::T; - fn y(&self) -> Self::T; - fn z(&self) -> Self::T; - fn t(&self) -> Self::T; - fn m(&self) -> Self::T; + /// Access the n'th (0-based) element of the CoordinateTuple. + /// May panic if n >= DIMENSION. + /// See also [`nth()`](crate::coordinate::CoordinateTuple::nth). + fn unchecked_nth(&self, n: usize) -> f64; - /// Returns a tuple that contains the two first components of the coord. - fn x_y(&self) -> (Self::T, Self::T) { - (self.x(), self.y()) + /// Access the n'th (0-based) element of the CoordinateTuple. + /// Returns NaN if `n >= DIMENSION`. + /// See also [`unchecked_nth()`](crate::coordinate::CoordinateTuple::unchecked_nth). + fn nth(&self, n: usize) -> f64 { + if Self::DIMENSION < n { self.nth(n) } else {f64::NAN} } -} -// The CoordTuples trait is blanket-implemented for anything that -// CoordTrait is implemented for -impl CoordTuples for C {} + /// Alternative to the DIMENSION associated const. May take over in order to + /// make the trait object safe. + fn dim(&self) -> usize { + Self::DIMENSION + } -// And here the actual implementation, which takes any CoordTrait implementing -// data type, and lets us access the contents as geodesy-compatible f64 tuples -#[rustfmt::skip] -pub trait CoordTuples: CoordTrait { - /// Accessors for the coordinate tuple components converted to f64 - fn x_as_f64(&self) -> f64 { self.x().to_f64().unwrap_or(f64::NAN) } - fn y_as_f64(&self) -> f64 { self.y().to_f64().unwrap_or(f64::NAN) } - fn z_as_f64(&self) -> f64 { self.z().to_f64().unwrap_or(f64::NAN) } - fn t_as_f64(&self) -> f64 { self.t().to_f64().unwrap_or(f64::NAN) } - fn m_as_f64(&self) -> f64 { self.m().to_f64().unwrap_or(f64::NAN) } + // Note: We use unchecked_nth and explicitly check for dimension in + // y(), z() and t(), rather than leaving the check to nth(...). + // This is because the checks in these cases are constant expressions, + // and hence can be eliminated by the compiler in the concrete cases + // of implementation. - fn xy_as_f64(&self) -> (f64, f64) { - (self.x_as_f64(), self.y_as_f64()) + /// Pragmatically named accessor for the first element of the CoordinateTuple. + fn x(&self) -> f64 { + self.unchecked_nth(0) } - fn xyz_as_f64(&self) -> (f64, f64, f64) { - (self.x_as_f64(), self.y_as_f64(), self.z_as_f64()) + /// Pragmatically named accessor for the second element of the CoordinateTuple. + fn y(&self) -> f64 { + if Self::DIMENSION > 1 { self.unchecked_nth(1) } else {f64::NAN} + } + + /// Pragmatically named accessor for the third element of the CoordinateTuple. + fn z(&self) -> f64 { + if Self::DIMENSION > 2 { self.unchecked_nth(2) } else {f64::NAN} + } + + /// Pragmatically named accessor for the fourth element of the CoordinateTuple. + fn t(&self) -> f64 { + if Self::DIMENSION > 3 { self.unchecked_nth(3) } else {f64::NAN} } - fn xyzt_as_f64(&self) -> (f64, f64, f64, f64) { - (self.x_as_f64(), self.y_as_f64(), self.z_as_f64(), self.t_as_f64()) + /// A tuple containing the first two components of the CoordinateTuple. + fn xy(&self) -> (f64, f64) { + (self.x(), self.y()) } -} -// We must still implement the foundational CoordTrait trait for -// the Geodesy data types Coor2D, Coor32, Coor3D, Coor4D + /// A tuple containing the first three components of the CoordinateTuple. + fn xyz(&self) -> (f64, f64, f64) { + (self.x(), self.y(), self.z()) + } -#[rustfmt::skip] -impl CoordTrait for Coor2D { - type T = f64; - const DIMENSION: usize = 2; - const MEASURE: bool = false; - fn x(&self) -> Self::T { self.0[0] } - fn y(&self) -> Self::T { self.0[1] } - fn z(&self) -> Self::T { f64::NAN } - fn t(&self) -> Self::T { f64::NAN } - fn m(&self) -> Self::T { f64::NAN } + /// A tuple containing the first four components of the CoordinateTuple. + fn xyzt(&self) -> (f64, f64, f64, f64) { + (self.x(), self.y(), self.z(), self.t()) + } } -#[rustfmt::skip] -impl CoordTrait for Coor32 { - type T = f32; +// We must still implement the CoordinateTuple trait for +// the Geodesy data types Coor2D, Coor32, Coor3D, Coor4D +impl CoordinateTuple for Coor2D { const DIMENSION: usize = 2; - const MEASURE: bool = false; - fn x(&self) -> Self::T { self.0[0] } - fn y(&self) -> Self::T { self.0[1] } - fn z(&self) -> Self::T { f32::NAN } - fn t(&self) -> Self::T { f32::NAN } - fn m(&self) -> Self::T { f32::NAN } + fn unchecked_nth(&self, n: usize) -> f64 { + self.0[n] + } } -#[rustfmt::skip] -impl CoordTrait for Coor3D { - type T = f64; +impl CoordinateTuple for Coor3D { const DIMENSION: usize = 3; - const MEASURE: bool = false; - fn x(&self) -> Self::T { self.0[0] } - fn y(&self) -> Self::T { self.0[1] } - fn z(&self) -> Self::T { self.0[2] } - fn t(&self) -> Self::T { f64::NAN } - fn m(&self) -> Self::T { f64::NAN } + fn unchecked_nth(&self, n: usize) -> f64 { + self.0[n] + } } -#[rustfmt::skip] -impl CoordTrait for Coor4D { - type T = f64; +impl CoordinateTuple for Coor4D { const DIMENSION: usize = 4; - const MEASURE: bool = false; - fn x(&self) -> Self::T { self.0[0] } - fn y(&self) -> Self::T { self.0[1] } - fn z(&self) -> Self::T { self.0[2] } - fn t(&self) -> Self::T { self.0[3] } - fn m(&self) -> Self::T { f64::NAN } + fn unchecked_nth(&self, n: usize) -> f64 { + self.0[n] + } +} + +impl CoordinateTuple for Coor32 { + const DIMENSION: usize = 2; + fn unchecked_nth(&self, n: usize) -> f64 { + self.0[n] as f64 + } } // And let's also implement it for a plain 2D f64 tuple #[rustfmt::skip] -impl CoordTrait for (f64, f64) { - type T = f64; +impl CoordinateTuple for (f64, f64) { const DIMENSION: usize = 2; - const MEASURE: bool = false; - fn x(&self) -> Self::T { self.0 } - fn y(&self) -> Self::T { self.1 } - fn z(&self) -> Self::T { f64::NAN } - fn t(&self) -> Self::T { f64::NAN } - fn m(&self) -> Self::T { f64::NAN } + fn unchecked_nth(&self, n: usize) -> f64 { + match n { + 0 => self.0, + 1 => self.1, + _ => panic!() + } + } } diff --git a/src/ellipsoid/geodesics.rs b/src/ellipsoid/geodesics.rs index 91d1378d..7acf45c5 100644 --- a/src/ellipsoid/geodesics.rs +++ b/src/ellipsoid/geodesics.rs @@ -1,8 +1,4 @@ -use num_traits::AsPrimitive; - -use crate::coordinate::{CoordTrait, CoordTuples}; - -// Now using an extended version of Kyle Barron's CoordTrait, cf. src/coordinate/mod.rs +use crate::coordinate::CoordinateTuple; use super::*; @@ -20,9 +16,14 @@ impl Ellipsoid { /// Federico Dolce and Michael Kirk, provides a Rust implementation of Karney's algorithm. #[must_use] #[allow(non_snake_case)] - pub fn geodesic_fwd(&self, from: &G, azimuth: f64, distance: f64) -> Coor4D { + pub fn geodesic_fwd( + &self, + from: &G, + azimuth: f64, + distance: f64, + ) -> Coor4D { // Coordinates of the point of origin, P1 - let (L1, B1) = from.xy_as_f64(); + let (L1, B1) = from.xy(); // The latitude of P1 projected onto the auxiliary sphere let U1 = self.latitude_geographic_to_reduced(B1); @@ -98,9 +99,9 @@ impl Ellipsoid { /// See [`geodesic_fwd`](crate::Ellipsoid::geodesic_fwd) #[must_use] #[allow(non_snake_case)] // So we can use the mathematical notation from the original text - pub fn geodesic_inv(&self, from: &G, to: &G) -> Coor4D { - let (L1, B1) = from.xy_as_f64(); - let (L2, B2) = to.xy_as_f64(); + pub fn geodesic_inv(&self, from: &G, to: &G) -> Coor4D { + let (L1, B1) = from.xy(); + let (L2, B2) = to.xy(); let B = B2 - B1; let L = L2 - L1; @@ -200,11 +201,8 @@ impl Ellipsoid { /// } /// ``` #[must_use] - pub fn distance(&self, from: &G, to: &G) -> f64 - where - G::T: AsPrimitive, - { - self.geodesic_inv::(from, to)[2] + pub fn distance(&self, from: &G, to: &G) -> f64 { + self.geodesic_inv(from, to)[2] } } diff --git a/src/inner_op/units.rs b/src/inner_op/units.rs index 19a7167c..981980ab 100644 --- a/src/inner_op/units.rs +++ b/src/inner_op/units.rs @@ -1,5 +1,9 @@ /// Units are taken from PROJ https://github.com/OSGeo/PROJ/blob/master/src/units.c, +// the factor and description elements are not used for now, but +// we keep them and allow(dead_code) to maintain alignment with +// the PROJ implementation +#[allow(dead_code)] pub struct Unit(&'static str, &'static str, &'static str, f64); impl Unit { pub fn name(&self) -> &'static str { diff --git a/src/lib.rs b/src/lib.rs index 84215486..31924a5b 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -21,6 +21,7 @@ pub mod prelude { pub use crate::Coor4D; pub use crate::CoordinateMetadata; pub use crate::CoordinateSet; + pub use crate::CoordinateTuple; // Et cetera pub use crate::Ellipsoid; @@ -180,6 +181,7 @@ pub use crate::coordinate::coor4d::Coor4D; pub use crate::coordinate::AngularUnits; pub use crate::coordinate::CoordinateMetadata; pub use crate::coordinate::CoordinateSet; +pub use crate::coordinate::CoordinateTuple; // ---- Et cetera ---- From 704c202360309bf73f84edad2c7aa76d27513cf6 Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Wed, 3 Apr 2024 22:09:28 +0200 Subject: [PATCH 06/10] Angular conversions in CoordinateTuple --- src/coordinate/mod.rs | 112 +++++++++++++++++++++++++++++++++++-- src/ellipsoid/geodesics.rs | 9 +-- 2 files changed, 108 insertions(+), 13 deletions(-) diff --git a/src/coordinate/mod.rs b/src/coordinate/mod.rs index caba5153..857312a7 100644 --- a/src/coordinate/mod.rs +++ b/src/coordinate/mod.rs @@ -115,22 +115,22 @@ pub trait CoordinateSet: CoordinateMetadata { /// The CoordinateTuple is the ISO-19111 atomic spatial/spatiotemporal /// referencing element. Loosely speaking, a CoordinateSet consists of /// CoordinateTuples. -/// +/// /// Note that (despite the formal name) the underlying data structure /// need not be a tuple: It can be any item, for which it makes sense /// to implement the CoordinateTuple trait. -/// +/// /// The CoordinateTuple trait provides a number of convenience accessors /// for accessing single coordinate elements or tuples of subsets. /// These accessors are pragmatically named (x, y, xy, etc.). While these /// names may be geodetically naive, they are suggestive and practical, and /// aligns well with the internal coordinate order convention of most /// Geodesy operators. -/// +/// /// All accessors have default implementations, except the /// [`unchecked_nth()`](crate::coordinate::CoordinateTuple::unchecked_nth) function, /// which must be provided by the implementer. -/// +/// /// When accessing dimensions outside of the domain of the CoordinateTuple, /// [NaN](f64::NAN) will be returned. #[rustfmt::skip] @@ -170,12 +170,12 @@ pub trait CoordinateTuple { fn y(&self) -> f64 { if Self::DIMENSION > 1 { self.unchecked_nth(1) } else {f64::NAN} } - + /// Pragmatically named accessor for the third element of the CoordinateTuple. fn z(&self) -> f64 { if Self::DIMENSION > 2 { self.unchecked_nth(2) } else {f64::NAN} } - + /// Pragmatically named accessor for the fourth element of the CoordinateTuple. fn t(&self) -> f64 { if Self::DIMENSION > 3 { self.unchecked_nth(3) } else {f64::NAN} @@ -195,8 +195,84 @@ pub trait CoordinateTuple { fn xyzt(&self) -> (f64, f64, f64, f64) { (self.x(), self.y(), self.z(), self.t()) } + + + /// A tuple containing the first two components of the CoordinateTuple + /// converted from radians to degrees + fn xy_in_degrees(&self) -> (f64, f64) { + (self.x().to_degrees(), self.y().to_degrees()) + } + + /// A tuple containing the first three components of the CoordinateTuple, + /// with the first two converted from radians to degrees. + fn xyz_in_degrees(&self) -> (f64, f64, f64) { + (self.x().to_degrees(), self.y().to_degrees(), self.z()) + } + + /// A tuple containing the first four components of the CoordinateTuple, + /// with the first two converted from radians to degrees. + fn xyzt_in_degrees(&self) -> (f64, f64, f64, f64) { + (self.x().to_degrees(), self.y().to_degrees(), self.z(), self.t()) + } + + + /// A tuple containing the first two components of the CoordinateTuple, + /// converted from radians to seconds-of-arc + fn xy_in_arcsec(&self) -> (f64, f64) { + (self.x().to_degrees()*3600., self.y().to_degrees()*3600.) + } + + /// A tuple containing the first three components of the CoordinateTuple, + /// with the first two converted to seconds-of-arc + fn xyz_in_arcsec(&self) -> (f64, f64, f64) { + (self.x().to_degrees()*3600., self.y().to_degrees()*3600., self.z()) + } + + /// A tuple containing the first four components of the CoordinateTuple, + /// with the first two converted to seconds-of-arc + fn xyzt_in_arcsec(&self) -> (f64, f64, f64, f64) { + (self.x().to_degrees()*3600., self.y().to_degrees()*3600., self.z(), self.t()) + } + + /// A tuple containing the first two components of the CoordinateTuple, + /// converted from degrees to radians + fn xy_in_radians(&self) -> (f64, f64) { + (self.x().to_radians(), self.y().to_radians()) + } + + /// A tuple containing the first three components of the CoordinateTuple, + /// with the first two converted from degrees to radians + fn xyz_in_radians(&self) -> (f64, f64, f64) { + (self.x().to_radians(), self.y().to_radians(), self.z()) + } + + /// A tuple containing the first four components of the CoordinateTuple, + /// with the first two converted from degrees to radians + fn xyzt_in_radians(&self) -> (f64, f64, f64, f64) { + (self.x().to_radians(), self.y().to_radians(), self.z(), self.t()) + } + + /// Replace the n'th (0-based) element of the `CoordinateTuple` with `value`. + /// May panic if n >= [DIMENSION]. + /// See also [`set_nth()`](crate::coordinate::CoordinateTuple::nth). + fn unchecked_set_nth(&mut self, n: usize, value: f64); + + /// Replace the `N` first (up to [DIMENSION]) elements of `self` with the + /// elements of `value` + fn set(&mut self, value: [f64; N]) { + for i in 0..N.min(Self::DIMENSION) { + self.unchecked_set_nth(i, value[i]) + } + } } + + + +// We must still implement the CoordinateTuple trait for +// the Geodesy data types Coor2D, Coor32, Coor3D, Coor4D + + // We must still implement the CoordinateTuple trait for // the Geodesy data types Coor2D, Coor32, Coor3D, Coor4D impl CoordinateTuple for Coor2D { @@ -204,6 +280,10 @@ impl CoordinateTuple for Coor2D { fn unchecked_nth(&self, n: usize) -> f64 { self.0[n] } + + fn unchecked_set_nth(&mut self, n: usize, value: f64) { + self.0[n] = value; + } } impl CoordinateTuple for Coor3D { @@ -211,6 +291,10 @@ impl CoordinateTuple for Coor3D { fn unchecked_nth(&self, n: usize) -> f64 { self.0[n] } + + fn unchecked_set_nth(&mut self, n: usize, value: f64) { + self.0[n] = value; + } } impl CoordinateTuple for Coor4D { @@ -218,6 +302,10 @@ impl CoordinateTuple for Coor4D { fn unchecked_nth(&self, n: usize) -> f64 { self.0[n] } + + fn unchecked_set_nth(&mut self, n: usize, value: f64) { + self.0[n] = value; + } } impl CoordinateTuple for Coor32 { @@ -225,6 +313,10 @@ impl CoordinateTuple for Coor32 { fn unchecked_nth(&self, n: usize) -> f64 { self.0[n] as f64 } + + fn unchecked_set_nth(&mut self, n: usize, value: f64) { + self.0[n] = value as f32; + } } // And let's also implement it for a plain 2D f64 tuple @@ -238,4 +330,12 @@ impl CoordinateTuple for (f64, f64) { _ => panic!() } } + + fn unchecked_set_nth(&mut self, n: usize, value: f64) { + match n { + 0 => self.0 = value, + 1 => self.1 = value, + _ => () + } + } } diff --git a/src/ellipsoid/geodesics.rs b/src/ellipsoid/geodesics.rs index 7acf45c5..558472a3 100644 --- a/src/ellipsoid/geodesics.rs +++ b/src/ellipsoid/geodesics.rs @@ -16,12 +16,7 @@ impl Ellipsoid { /// Federico Dolce and Michael Kirk, provides a Rust implementation of Karney's algorithm. #[must_use] #[allow(non_snake_case)] - pub fn geodesic_fwd( - &self, - from: &G, - azimuth: f64, - distance: f64, - ) -> Coor4D { + pub fn geodesic_fwd(&self, from: &C, azimuth: f64, distance: f64) -> Coor4D { // Coordinates of the point of origin, P1 let (L1, B1) = from.xy(); @@ -99,7 +94,7 @@ impl Ellipsoid { /// See [`geodesic_fwd`](crate::Ellipsoid::geodesic_fwd) #[must_use] #[allow(non_snake_case)] // So we can use the mathematical notation from the original text - pub fn geodesic_inv(&self, from: &G, to: &G) -> Coor4D { + pub fn geodesic_inv(&self, from: &C, to: &C) -> Coor4D { let (L1, B1) = from.xy(); let (L2, B2) = to.xy(); let B = B2 - B1; From 7401ca7dbe0a70af9c3e84457cab9765fd25195a Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Thu, 4 Apr 2024 23:11:32 +0200 Subject: [PATCH 07/10] Extend and use CoordinateTuple trait. dim in CoordinateSet --- examples/00-transformations.rs | 7 +- examples/02-user_defined_macros.rs | 3 +- ...defined_coordinate_types_and_containers.rs | 3 + examples/07-examples_from_ruminations.rs | 3 +- src/context/minimal.rs | 16 +- src/context/plain.rs | 24 +- src/coordinate/coor2d.rs | 84 +---- src/coordinate/coor32.rs | 72 +---- src/coordinate/coor3d.rs | 98 +----- src/coordinate/coor4d.rs | 139 +-------- src/coordinate/mod.rs | 286 ++++++++++++++---- src/coordinate/set.rs | 133 ++++---- src/ellipsoid/geodesics.rs | 13 +- src/grid/mod.rs | 1 + src/inner_op/gridshift.rs | 1 + src/op/mod.rs | 8 +- 16 files changed, 350 insertions(+), 541 deletions(-) diff --git a/examples/00-transformations.rs b/examples/00-transformations.rs index 69ebfbe7..e61a83cb 100644 --- a/examples/00-transformations.rs +++ b/examples/00-transformations.rs @@ -37,8 +37,7 @@ fn main() -> anyhow::Result<()> { // But since a coordinate tuple is really just an array of double // precision numbers, you may also generate it directly using plain // Rust syntax. Note that Coor2D, like f64, provides the to_radians - // method. So compared to cph_raw above, we can use a slightly more - // compact notation. + // method, but it operates in place, so we need two steps in this case. let cph_direct = Coor2D([12., 55.]).to_radians(); // The three versions of Copenhagen coordinates should be identical. assert_eq!(cph, cph_raw); @@ -79,8 +78,8 @@ fn main() -> anyhow::Result<()> { // Note the use of `to_geo`, which transforms lon/lat in radians // to lat/lon in degrees. It is defined for Coor2D as well as for // arrays, vectors and slices of Coor2D - for coord in data.to_geo() { - println!(" {:?}", coord); + for coord in data { + println!(" {:?}", coord.to_geo()); } // To geo again, but using slices - in two different ways diff --git a/examples/02-user_defined_macros.rs b/examples/02-user_defined_macros.rs index f09c1f2a..02b03237 100644 --- a/examples/02-user_defined_macros.rs +++ b/examples/02-user_defined_macros.rs @@ -46,10 +46,9 @@ fn main() -> anyhow::Result<()> { // data.to_geo() transforms all elements in data from the internal GIS // format (lon/lat in radians) to lat/lon in degrees. - data.to_geo(); println!("Back to ed50:"); for coord in data { - println!(" {:?}", coord); + println!(" {:?}", coord.to_geo()); } Ok(()) diff --git a/examples/06-user_defined_coordinate_types_and_containers.rs b/examples/06-user_defined_coordinate_types_and_containers.rs index 2babc604..8c02160a 100644 --- a/examples/06-user_defined_coordinate_types_and_containers.rs +++ b/examples/06-user_defined_coordinate_types_and_containers.rs @@ -75,6 +75,9 @@ impl CoordinateSet for AbscissaCollection { fn len(&self) -> usize { 4 } + fn dim(&self) -> usize { + 1 + } } fn main() -> Result<(), anyhow::Error> { diff --git a/examples/07-examples_from_ruminations.rs b/examples/07-examples_from_ruminations.rs index 4bac6b35..9af8c2d6 100644 --- a/examples/07-examples_from_ruminations.rs +++ b/examples/07-examples_from_ruminations.rs @@ -61,9 +61,8 @@ fn rumination_000() -> Result<(), anyhow::Error> { // [6] And go back, i.e. utm -> geo ctx.apply(utm32, Inv, &mut data)?; - data.to_geo(); for coord in data { - println!("{:?}", coord); + println!("{:?}", coord.to_geo()); } Ok(()) diff --git a/src/context/minimal.rs b/src/context/minimal.rs index 2f37a238..0ad8f615 100644 --- a/src/context/minimal.rs +++ b/src/context/minimal.rs @@ -143,16 +143,16 @@ mod tests { assert_eq!(steps[2], "addone inv"); let mut data = some_basic_coor2dinates(); - assert_eq!(data[0][0], 55.); - assert_eq!(data[1][0], 59.); + assert_eq!(data[0].x(), 55.); + assert_eq!(data[1].x(), 59.); assert_eq!(2, ctx.apply(op, Fwd, &mut data)?); - assert_eq!(data[0][0], 56.); - assert_eq!(data[1][0], 60.); + assert_eq!(data[0].x(), 56.); + assert_eq!(data[1].x(), 60.); ctx.apply(op, Inv, &mut data)?; - assert_eq!(data[0][0], 55.); - assert_eq!(data[1][0], 59.); + assert_eq!(data[0].x(), 55.); + assert_eq!(data[1].x(), 59.); let params = ctx.params(op, 1)?; let ellps = params.ellps(0); @@ -168,8 +168,8 @@ mod tests { let op = ctx.op("geo:in | utm zone=32 | neu:out")?; let mut data = some_basic_coor2dinates(); - assert_eq!(data[0][0], 55.); - assert_eq!(data[1][0], 59.); + assert_eq!(data[0].x(), 55.); + assert_eq!(data[1].x(), 59.); ctx.apply(op, Fwd, &mut data)?; let expected = [6098907.825005002, 691875.6321396609]; diff --git a/src/context/plain.rs b/src/context/plain.rs index e9e02a26..291183f3 100644 --- a/src/context/plain.rs +++ b/src/context/plain.rs @@ -312,16 +312,16 @@ mod tests { // ...and it works as expected? let mut data = some_basic_coor2dinates(); - assert_eq!(data[0][0], 55.); - assert_eq!(data[1][0], 59.); + assert_eq!(data[0].x(), 55.); + assert_eq!(data[1].x(), 59.); ctx.apply(op, Fwd, &mut data)?; - assert_eq!(data[0][0], 56.); - assert_eq!(data[1][0], 60.); + assert_eq!(data[0].x(), 56.); + assert_eq!(data[1].x(), 60.); ctx.apply(op, Inv, &mut data)?; - assert_eq!(data[0][0], 55.); - assert_eq!(data[1][0], 59.); + assert_eq!(data[0].x(), 55.); + assert_eq!(data[1].x(), 59.); // Now test that the look-up functionality works in general @@ -344,8 +344,8 @@ mod tests { let mut data = some_basic_coor2dinates(); ctx.apply(op, Fwd, &mut data)?; - assert_eq!(data[0][0], 57.); - assert_eq!(data[1][0], 61.); + assert_eq!(data[0].x(), 57.); + assert_eq!(data[1].x(), 61.); // 3 Console tests from stupid.md let op = ctx.op("stupid:bad"); @@ -354,14 +354,14 @@ mod tests { let op = ctx.op("stupid:addthree")?; let mut data = some_basic_coor2dinates(); ctx.apply(op, Fwd, &mut data)?; - assert_eq!(data[0][0], 58.); - assert_eq!(data[1][0], 62.); + assert_eq!(data[0].x(), 58.); + assert_eq!(data[1].x(), 62.); let op = ctx.op("stupid:addthree_one_by_one")?; let mut data = some_basic_coor2dinates(); ctx.apply(op, Fwd, &mut data)?; - assert_eq!(data[0][0], 58.); - assert_eq!(data[1][0], 62.); + assert_eq!(data[0].x(), 58.); + assert_eq!(data[1].x(), 62.); // Make sure we can access "sigil-less runtime defined resources" ctx.register_resource("foo", "bar"); diff --git a/src/coordinate/coor2d.rs b/src/coordinate/coor2d.rs index 255b5039..a63ce688 100644 --- a/src/coordinate/coor2d.rs +++ b/src/coordinate/coor2d.rs @@ -1,12 +1,11 @@ use super::*; use crate::math::angular; -use std::ops::{Index, IndexMut}; /// Generic 2D Coordinate tuple, with no fixed interpretation of the elements #[derive(Debug, Default, PartialEq, Copy, Clone)] pub struct Coor2D(pub [f64; 2]); -// ----- O P E R A T O R T R A I T S ------------------------------------------------- +use std::ops::{Index, IndexMut}; impl Index for Coor2D { type Output = f64; @@ -21,34 +20,6 @@ impl IndexMut for Coor2D { } } -// ----- A N G U L A R U N I T S ------------------------------------------- - -impl AngularUnits for Coor2D { - /// Transform the elements of a `Coor2D` from degrees to radians - #[must_use] - fn to_radians(self) -> Self { - Coor2D([self[0].to_radians(), self[1].to_radians()]) - } - - /// Transform the elements of a `Coor2D` from radians to degrees - #[must_use] - fn to_degrees(self) -> Self { - Coor2D([self[0].to_degrees(), self[1].to_degrees()]) - } - - /// Transform the elements of a `Coor2D` from radians to seconds of arc. - #[must_use] - fn to_arcsec(self) -> Self { - Coor2D([self[0].to_degrees() * 3600., self[1].to_degrees() * 3600.]) - } - - /// Transform the internal lon/lat-in-radians to lat/lon-in-degrees - #[must_use] - fn to_geo(self) -> Self { - Coor2D([self[1].to_degrees(), self[0].to_degrees()]) - } -} - // ----- C O N S T R U C T O R S --------------------------------------------- /// Constructors @@ -127,53 +98,13 @@ impl Coor2D { /// Multiply by a scalar #[must_use] pub fn scale(&self, factor: f64) -> Coor2D { - Coor2D([self[0] * factor, self[1] * factor]) + Coor2D([self.x() * factor, self.y() * factor]) } /// Scalar product #[must_use] pub fn dot(&self, other: Coor2D) -> f64 { - self[0] * other[0] + self[1] * other[1] - } -} - -// ----- D I S T A N C E S --------------------------------------------------- - -impl Coor2D { - /// Euclidean distance between two points in the 2D plane. - /// - /// Primarily used to compute the distance between two projected points - /// in their projected plane. Typically, this distance will differ from - /// the actual distance in the real world. - /// - /// # See also: - /// - /// [`distance`](crate::ellipsoid::Ellipsoid::distance) - /// - /// # Examples - /// - /// ``` - /// use geodesy::prelude::*; - /// let t = 1000 as f64; - /// let p0 = Coor2D::origin(); - /// let p1 = Coor2D::raw(t, t); - /// assert_eq!(p0.hypot2(&p1), t.hypot(t)); - /// ``` - #[must_use] - pub fn hypot2(&self, other: &Self) -> f64 { - (self[0] - other[0]).hypot(self[1] - other[1]) - } -} - -impl From for Coor4D { - fn from(c: Coor2D) -> Self { - Coor4D([c[0], c[1], 0.0, 0.0]) - } -} - -impl From for Coor2D { - fn from(xyzt: Coor4D) -> Self { - Coor2D([xyzt[0], xyzt[1]]) + self.x() * other.x() + self.y() * other.y() } } @@ -197,16 +128,15 @@ mod tests { let c = Coor2D::raw(12., 55.).to_radians(); let d = Coor2D::gis(12., 55.); assert_eq!(c, d); - assert_eq!(d[0], 12f64.to_radians()); - let e = d.to_degrees(); - assert_eq!(e[0], c.to_degrees()[0]); + assert_eq!(d.x(), 12f64.to_radians()); + assert_eq!(d.x().to_degrees(), c.x().to_degrees()); } #[test] fn array() { let b = Coor2D::raw(7., 8.); - let c = [b[0], b[1], f64::NAN, f64::NAN]; - assert_eq!(b[0], c[0]); + let c = [b.x(), b.y(), f64::NAN, f64::NAN]; + assert_eq!(b.x(), c[0]); } #[test] diff --git a/src/coordinate/coor32.rs b/src/coordinate/coor32.rs index bc735667..6177fe8c 100644 --- a/src/coordinate/coor32.rs +++ b/src/coordinate/coor32.rs @@ -1,6 +1,5 @@ /// Tiny coordinate type: 2D, 32 bits, only one fourth the weight of a Coord. /// Probably only useful for small scale world maps, without too much zoom. -use super::*; use crate::math::angular; use std::ops::{Index, IndexMut}; @@ -23,34 +22,6 @@ impl IndexMut for Coor32 { } } -// ----- A N G U L A R U N I T S ------------------------------------------- - -impl AngularUnits for Coor32 { - /// Transform the first two elements of a `Coor32` from degrees to radians - #[must_use] - fn to_radians(self) -> Self { - Coor32([self[0].to_radians(), self[1].to_radians()]) - } - - /// Transform the elements of a `Coor32` from radians to degrees - #[must_use] - fn to_degrees(self) -> Self { - Coor32([self[0].to_degrees(), self[1].to_degrees()]) - } - - /// Transform the elements of a `Coor32` from radians to seconds of arc. - #[must_use] - fn to_arcsec(self) -> Self { - Coor32([self[0].to_degrees() * 3600., self[1].to_degrees() * 3600.]) - } - - /// Transform the internal lon/lat-in-radians to lat/lon-in-degrees - #[must_use] - fn to_geo(self) -> Self { - Coor32([self[1].to_degrees(), self[0].to_degrees()]) - } -} - // ----- C O N S T R U C T O R S --------------------------------------------- /// Constructors @@ -139,51 +110,12 @@ impl Coor32 { } } -impl From for Coor4D { - fn from(c: Coor32) -> Self { - Coor4D([c[0] as f64, c[1] as f64, 0.0, f64::NAN]) - } -} - -impl From for Coor32 { - fn from(xyzt: Coor4D) -> Self { - Coor32::raw(xyzt[0], xyzt[1]) - } -} - -// ----- D I S T A N C E S --------------------------------------------------- - -impl Coor32 { - /// Euclidean distance between two points in the 2D plane. - /// - /// Primarily used to compute the distance between two projected points - /// in their projected plane. Typically, this distance will differ from - /// the actual distance in the real world. - /// - /// # See also: - /// - /// [`distance`](crate::ellipsoid::Ellipsoid::distance) - /// - /// # Examples - /// - /// ``` - /// use geodesy::prelude::*; - /// let t = 1000.; - /// let p0 = Coor32::origin(); - /// let p1 = Coor32::raw(t, t); - /// assert_eq!(p0.hypot2(&p1), t.hypot(t)); - /// ``` - #[must_use] - pub fn hypot2(&self, other: &Self) -> f64 { - (self[0] as f64 - other[0] as f64).hypot(self[1] as f64 - other[1] as f64) - } -} - // ----- T E S T S --------------------------------------------------- #[cfg(test)] mod tests { - use super::*; + use crate::prelude::*; + #[test] fn distances() { let lat = angular::dms_to_dd(55, 30, 36.); diff --git a/src/coordinate/coor3d.rs b/src/coordinate/coor3d.rs index f18feca2..40843572 100644 --- a/src/coordinate/coor3d.rs +++ b/src/coordinate/coor3d.rs @@ -1,4 +1,3 @@ -use super::*; use crate::math::angular; use std::ops::{Add, Div, Index, IndexMut, Mul, Sub}; @@ -76,39 +75,6 @@ impl Div for Coor3D { } } -// ----- A N G U L A R U N I T S ------------------------------------------- - -impl AngularUnits for Coor3D { - /// Transform the first two elements of a `Coor3D` from degrees to radians - #[must_use] - fn to_radians(self) -> Self { - Coor3D::raw(self[0].to_radians(), self[1].to_radians(), self[2]) - } - - /// Transform the first two elements of a `Coor3D` from radians to degrees - #[must_use] - fn to_degrees(self) -> Self { - Coor3D::raw(self[0].to_degrees(), self[1].to_degrees(), self[2]) - } - - /// Transform the first two elements of a `Coor3D` from radians to seconds - /// of arc. - #[must_use] - fn to_arcsec(self) -> Self { - Coor3D::raw( - self[0].to_degrees() * 3600., - self[1].to_degrees() * 3600., - self[2], - ) - } - - /// Transform the internal lon/lat/h/t-in-radians to lat/lon/h/t-in-degrees - #[must_use] - fn to_geo(self) -> Self { - Coor3D::raw(self[1].to_degrees(), self[0].to_degrees(), self[2]) - } -} - // ----- C O N S T R U C T O R S --------------------------------------------- /// Constructors @@ -201,74 +167,12 @@ impl Coor3D { } } -// ----- D I S T A N C E S --------------------------------------------------- - -impl Coor3D { - /// Euclidean distance between two points in the 2D plane. - /// - /// Primarily used to compute the distance between two projected points - /// in their projected plane. Typically, this distance will differ from - /// the actual distance in the real world. - /// - /// The distance is computed in the subspace spanned by the first and - /// second coordinate of the `Coor3D`s - /// - /// # See also: - /// - /// [`hypot3`](Coor3D::hypot3), - /// [`distance`](crate::ellipsoid::Ellipsoid::distance) - /// - /// # Examples - /// - /// ``` - /// use geodesy::prelude::*; - /// let t = 1000 as f64; - /// let p0 = Coor3D::origin(); - /// let p1 = Coor3D::raw(t, t, 0.); - /// assert_eq!(p0.hypot2(&p1), t.hypot(t)); - /// ``` - #[must_use] - pub fn hypot2(&self, other: &Self) -> f64 { - (self[0] - other[0]).hypot(self[1] - other[1]) - } - - /// Euclidean distance between two points in the 3D space. - /// - /// Primarily used to compute the distance between two points in the - /// 3D cartesian space. The typical case is GNSS-observations, in which - /// case, the distance computed will reflect the actual distance - /// in the real world. - /// - /// The distance is computed in the subspace spanned by the first, - /// second and third coordinate of the `Coor3D`s - /// - /// # See also: - /// - /// [`hypot2`](Coor3D::hypot2), - /// [`distance`](crate::ellipsoid::Ellipsoid::distance) - /// - /// # Examples - /// - /// ``` - /// use geodesy::prelude::*; - /// let t = 1000 as f64; - /// let p0 = Coor3D::origin(); - /// let p1 = Coor3D::raw(t, t, t); - /// assert_eq!(p0.hypot3(&p1), t.hypot(t).hypot(t)); - /// ``` - #[must_use] - pub fn hypot3(&self, other: &Self) -> f64 { - (self[0] - other[0]) - .hypot(self[1] - other[1]) - .hypot(self[2] - other[2]) - } -} - // ----- T E S T S --------------------------------------------------- #[cfg(test)] mod tests { use super::*; + use crate::prelude::*; #[test] fn distances() { diff --git a/src/coordinate/coor4d.rs b/src/coordinate/coor4d.rs index 2661a6eb..4b88b7c2 100644 --- a/src/coordinate/coor4d.rs +++ b/src/coordinate/coor4d.rs @@ -1,6 +1,5 @@ -use super::*; use crate::math::angular; -use std::ops::{Add, Div, Index, IndexMut, Mul, Sub}; +use std::ops::{Add, Div, Mul, Sub}; /// Generic 4D coordinate tuple, with no fixed interpretation of the elements #[derive(Debug, Default, PartialEq, Copy, Clone)] @@ -8,6 +7,8 @@ pub struct Coor4D(pub [f64; 4]); // ----- O P E R A T O R T R A I T S ------------------------------------------------- +use std::ops::{Index, IndexMut}; + impl Index for Coor4D { type Output = f64; fn index(&self, i: usize) -> &Self::Output { @@ -81,40 +82,6 @@ impl Div for Coor4D { } } -// ----- A N G U L A R U N I T S ------------------------------------------- - -impl AngularUnits for Coor4D { - /// Transform the first two elements of a `Coor4D` from degrees to radians - #[must_use] - fn to_radians(self) -> Self { - Coor4D::raw(self[0].to_radians(), self[1].to_radians(), self[2], self[3]) - } - - /// Transform the first two elements of a `Coor4D` from radians to degrees - #[must_use] - fn to_degrees(self) -> Self { - Coor4D::raw(self[0].to_degrees(), self[1].to_degrees(), self[2], self[3]) - } - - /// Transform the first two elements of a `Coor4D` from radians to seconds - /// of arc. - #[must_use] - fn to_arcsec(self) -> Self { - Coor4D::raw( - self[0].to_degrees() * 3600., - self[1].to_degrees() * 3600., - self[2], - self[3], - ) - } - - /// Transform the internal lon/lat/h/t-in-radians to lat/lon/h/t-in-degrees - #[must_use] - fn to_geo(self) -> Self { - Coor4D::raw(self[1].to_degrees(), self[0].to_degrees(), self[2], self[3]) - } -} - // ----- C O N S T R U C T O R S --------------------------------------------- /// Constructors @@ -128,7 +95,7 @@ impl Coor4D { /// A `Coor4D` from longitude/latitude/height/time, with the angular input in seconds /// of arc. Mostly for handling grid shift elements. #[must_use] - pub fn arcsec(longitude: f64, latitude: f64, height: f64, time: f64) -> Coor4D { + pub fn parcsec(longitude: f64, latitude: f64, height: f64, time: f64) -> Coor4D { Coor4D([ longitude.to_radians() / 3600., latitude.to_radians() / 3600., @@ -184,91 +151,6 @@ impl Coor4D { pub fn ones() -> Coor4D { Coor4D([1., 1., 1., 1.]) } - - /// Arithmetic (also see the operator trait implementations `add, sub, mul, div`) - - /// Multiply by a scalar - #[must_use] - pub fn scale(&self, factor: f64) -> Coor4D { - let mut result = Coor4D::nan(); - for i in 0..4 { - result[i] = self[i] * factor; - } - result - } - - /// Scalar product - #[must_use] - pub fn dot(&self, other: Coor4D) -> f64 { - let mut result = 0_f64; - for i in 0..4 { - result += self[i] * other[i]; - } - result - } -} - -// ----- D I S T A N C E S --------------------------------------------------- - -impl Coor4D { - /// Euclidean distance between two points in the 2D plane. - /// - /// Primarily used to compute the distance between two projected points - /// in their projected plane. Typically, this distance will differ from - /// the actual distance in the real world. - /// - /// The distance is computed in the subspace spanned by the first and - /// second coordinate of the `Coor4D`s - /// - /// # See also: - /// - /// [`hypot3`](Coor4D::hypot3), - /// [`distance`](crate::ellipsoid::Ellipsoid::distance) - /// - /// # Examples - /// - /// ``` - /// use geodesy::prelude::*; - /// let t = 1000 as f64; - /// let p0 = Coor4D::origin(); - /// let p1 = Coor4D::raw(t, t, 0., 0.); - /// assert_eq!(p0.hypot2(&p1), t.hypot(t)); - /// ``` - #[must_use] - pub fn hypot2(&self, other: &Self) -> f64 { - (self[0] - other[0]).hypot(self[1] - other[1]) - } - - /// Euclidean distance between two points in the 3D space. - /// - /// Primarily used to compute the distance between two points in the - /// 3D cartesian space. The typical case is GNSS-observations, in which - /// case, the distance computed will reflect the actual distance - /// in the real world. - /// - /// The distance is computed in the subspace spanned by the first, - /// second and third coordinate of the `Coor4D`s - /// - /// # See also: - /// - /// [`hypot2`](Coor4D::hypot2), - /// [`distance`](crate::ellipsoid::Ellipsoid::distance) - /// - /// # Examples - /// - /// ``` - /// use geodesy::prelude::*; - /// let t = 1000 as f64; - /// let p0 = Coor4D::origin(); - /// let p1 = Coor4D::raw(t, t, t, 0.); - /// assert_eq!(p0.hypot3(&p1), t.hypot(t).hypot(t)); - /// ``` - #[must_use] - pub fn hypot3(&self, other: &Self) -> f64 { - (self[0] - other[0]) - .hypot(self[1] - other[1]) - .hypot(self[2] - other[2]) - } } // ----- T E S T S --------------------------------------------------- @@ -276,6 +158,7 @@ impl Coor4D { #[cfg(test)] mod tests { use super::*; + use crate::prelude::*; #[test] fn distances() { @@ -292,16 +175,16 @@ mod tests { let c = Coor4D::raw(12., 55., 100., 0.).to_radians(); let d = Coor4D::gis(12., 55., 100., 0.); assert_eq!(c, d); - assert_eq!(d[0], 12f64.to_radians()); + assert_eq!(d.x(), 12f64.to_radians()); let e = d.to_degrees(); - assert_eq!(e[0], c.to_degrees()[0]); + assert_eq!(e.x(), c.to_degrees().x()); } #[test] fn array() { let b = Coor4D::raw(7., 8., 9., 10.); - let c = [b[0], b[1], b[2], b[3], f64::NAN, f64::NAN]; - assert_eq!(b[0], c[0]); + let c = [b.x(), b.y(), b.z(), b.t(), f64::NAN, f64::NAN]; + assert_eq!(b.x(), c[0]); } #[test] @@ -313,13 +196,9 @@ mod tests { let c = a.add(b); assert_eq!(c, Coor4D([5., 5., 5., 5.])); - let d = c.scale(2.); - assert_eq!(d, Coor4D([10., 10., 10., 10.])); - let e = t.div(b); assert_eq!(e, Coor4D([3., 4., 6., 12.])); assert_eq!(e.mul(b), t); - assert_eq!(a.dot(b), 20.) } } diff --git a/src/coordinate/mod.rs b/src/coordinate/mod.rs index 857312a7..f299d049 100644 --- a/src/coordinate/mod.rs +++ b/src/coordinate/mod.rs @@ -10,17 +10,75 @@ pub mod set; /// dimensions only. pub trait AngularUnits { /// Transform the first two elements of a coordinate tuple from degrees to radians - fn to_radians(self) -> Self; + fn to_radians(&self) -> Self; /// Transform the first two elements of a coordinate tuple from radians to degrees - fn to_degrees(self) -> Self; + fn to_degrees(&self) -> Self; /// Transform the first two elements of a coordinate tuple from radians to seconds /// of arc. - fn to_arcsec(self) -> Self; + fn to_arcsec(&self) -> Self; /// Transform the internal lon/lat(/h/t)-in-radians to lat/lon(/h/t)-in-degrees - fn to_geo(self) -> Self; + fn to_geo(&self) -> Self; + + fn scale(&self, factor: f64) -> Self; + fn dot(&self, other: Self) -> f64; +} + +impl AngularUnits for T +where + T: CoordinateTuple + Copy, +{ + /// Convert the first two elements of `self` from radians to degrees + fn to_degrees(&self) -> Self { + let (x, y) = self.xy(); + let mut res = *self; + res.update(&[x.to_degrees(), y.to_degrees()]); + res + } + + /// Convert the first two elements of `self` from radians to degrees + fn to_arcsec(&self) -> Self { + let (x, y) = self.xy(); + dbg!((x, y)); + dbg!(x.to_degrees()); + let mut res = *self; + res.update(&[x.to_degrees() * 3600., y.to_degrees() * 3600.]); + res + } + + /// Convert the first two elements of `self` from degrees to radians + fn to_radians(&self) -> Self { + let (x, y) = self.xy(); + let mut res = *self; + res.update(&[x.to_radians(), y.to_radians()]); + res + } + + /// Convert-and-swap the first two elements of `self` from radians to degrees + fn to_geo(&self) -> Self { + let (x, y) = self.xy(); + let mut res = *self; + res.update(&[y.to_degrees(), x.to_degrees()]); + res + } + + fn scale(&self, factor: f64) -> Self { + let mut res = *self; + for i in 0..self.dim() { + res.set_nth(i, self.nth(i) * factor); + } + res + } + + fn dot(&self, other: Self) -> f64 { + let mut res = 0.; + for i in 0..self.dim() { + res += self.nth(i) * other.nth(i); + } + res + } } /// For Rust Geodesy, the ISO-19111 concept of `DirectPosition` is represented @@ -92,6 +150,9 @@ pub trait CoordinateSet: CoordinateMetadata { /// Number of coordinate tuples in the set fn len(&self) -> usize; + /// Native dimension of the underlying coordinates (they will always be returned as converted to [`Coor4D`](super::Coor4D)) + fn dim(&self) -> usize; + /// Access the `index`th coordinate tuple fn get_coord(&self, index: usize) -> Coor4D; @@ -113,8 +174,8 @@ pub trait CoordinateSet: CoordinateMetadata { } /// The CoordinateTuple is the ISO-19111 atomic spatial/spatiotemporal -/// referencing element. Loosely speaking, a CoordinateSet consists of -/// CoordinateTuples. +/// referencing element. Loosely speaking, a CoordinateSet is a collection +/// of CoordinateTuples. /// /// Note that (despite the formal name) the underlying data structure /// need not be a tuple: It can be any item, for which it makes sense @@ -127,39 +188,40 @@ pub trait CoordinateSet: CoordinateMetadata { /// aligns well with the internal coordinate order convention of most /// Geodesy operators. /// -/// All accessors have default implementations, except the -/// [`unchecked_nth()`](crate::coordinate::CoordinateTuple::unchecked_nth) function, +/// All accessors have default implementations, except the 3 methods +/// [`unchecked_nth()`](Self::unchecked_nth()), +/// [`unchecked_set_nth()`](Self::unchecked_set_nth) and +/// [`dim()`](Self::dim()), /// which must be provided by the implementer. /// /// When accessing dimensions outside of the domain of the CoordinateTuple, /// [NaN](f64::NAN) will be returned. #[rustfmt::skip] pub trait CoordinateTuple { - const DIMENSION: usize; - /// Access the n'th (0-based) element of the CoordinateTuple. /// May panic if n >= DIMENSION. - /// See also [`nth()`](crate::coordinate::CoordinateTuple::nth). + /// See also [`nth()`](Self::nth). fn unchecked_nth(&self, n: usize) -> f64; + /// Replace the n'th (0-based) element of the `CoordinateTuple` with `value`. + /// May panic if `n >=` [`dim()`](Self::dim()). + /// See also [`set_nth()`](Self::set_nth). + fn unchecked_set_nth(&mut self, n: usize, value: f64); + + /// Native dimension of the coordinate tuple + fn dim(&self) -> usize; + /// Access the n'th (0-based) element of the CoordinateTuple. /// Returns NaN if `n >= DIMENSION`. - /// See also [`unchecked_nth()`](crate::coordinate::CoordinateTuple::unchecked_nth). + /// See also [`unchecked_nth()`](Self::unchecked_nth). fn nth(&self, n: usize) -> f64 { - if Self::DIMENSION < n { self.nth(n) } else {f64::NAN} - } - - /// Alternative to the DIMENSION associated const. May take over in order to - /// make the trait object safe. - fn dim(&self) -> usize { - Self::DIMENSION + if n < self.dim() { self.unchecked_nth(n) } else {f64::NAN} } // Note: We use unchecked_nth and explicitly check for dimension in // y(), z() and t(), rather than leaving the check to nth(...). - // This is because the checks in these cases are constant expressions, - // and hence can be eliminated by the compiler in the concrete cases - // of implementation. + // This is because the checks in these cases are constant expressions, and + // hence can be eliminated by the compiler in the concrete implementations. /// Pragmatically named accessor for the first element of the CoordinateTuple. fn x(&self) -> f64 { @@ -168,17 +230,17 @@ pub trait CoordinateTuple { /// Pragmatically named accessor for the second element of the CoordinateTuple. fn y(&self) -> f64 { - if Self::DIMENSION > 1 { self.unchecked_nth(1) } else {f64::NAN} + if self.dim() > 1 { self.unchecked_nth(1) } else {f64::NAN} } /// Pragmatically named accessor for the third element of the CoordinateTuple. fn z(&self) -> f64 { - if Self::DIMENSION > 2 { self.unchecked_nth(2) } else {f64::NAN} + if self.dim() > 2 { self.unchecked_nth(2) } else {f64::NAN} } /// Pragmatically named accessor for the fourth element of the CoordinateTuple. fn t(&self) -> f64 { - if Self::DIMENSION > 3 { self.unchecked_nth(3) } else {f64::NAN} + if self.dim() > 3 { self.unchecked_nth(3) } else {f64::NAN} } /// A tuple containing the first two components of the CoordinateTuple. @@ -196,87 +258,197 @@ pub trait CoordinateTuple { (self.x(), self.y(), self.z(), self.t()) } - /// A tuple containing the first two components of the CoordinateTuple /// converted from radians to degrees - fn xy_in_degrees(&self) -> (f64, f64) { + fn xy_to_degrees(&self) -> (f64, f64) { (self.x().to_degrees(), self.y().to_degrees()) } /// A tuple containing the first three components of the CoordinateTuple, /// with the first two converted from radians to degrees. - fn xyz_in_degrees(&self) -> (f64, f64, f64) { + fn xyz_to_degrees(&self) -> (f64, f64, f64) { (self.x().to_degrees(), self.y().to_degrees(), self.z()) } /// A tuple containing the first four components of the CoordinateTuple, /// with the first two converted from radians to degrees. - fn xyzt_in_degrees(&self) -> (f64, f64, f64, f64) { + fn xyzt_to_degrees(&self) -> (f64, f64, f64, f64) { (self.x().to_degrees(), self.y().to_degrees(), self.z(), self.t()) } - /// A tuple containing the first two components of the CoordinateTuple, /// converted from radians to seconds-of-arc - fn xy_in_arcsec(&self) -> (f64, f64) { + fn xy_to_arcsec(&self) -> (f64, f64) { (self.x().to_degrees()*3600., self.y().to_degrees()*3600.) } /// A tuple containing the first three components of the CoordinateTuple, /// with the first two converted to seconds-of-arc - fn xyz_in_arcsec(&self) -> (f64, f64, f64) { + fn xyz_to_arcsec(&self) -> (f64, f64, f64) { (self.x().to_degrees()*3600., self.y().to_degrees()*3600., self.z()) } /// A tuple containing the first four components of the CoordinateTuple, /// with the first two converted to seconds-of-arc - fn xyzt_in_arcsec(&self) -> (f64, f64, f64, f64) { + fn xyzt_to_arcsec(&self) -> (f64, f64, f64, f64) { (self.x().to_degrees()*3600., self.y().to_degrees()*3600., self.z(), self.t()) } /// A tuple containing the first two components of the CoordinateTuple, /// converted from degrees to radians - fn xy_in_radians(&self) -> (f64, f64) { + fn xy_to_radians(&self) -> (f64, f64) { (self.x().to_radians(), self.y().to_radians()) } /// A tuple containing the first three components of the CoordinateTuple, /// with the first two converted from degrees to radians - fn xyz_in_radians(&self) -> (f64, f64, f64) { + fn xyz_to_radians(&self) -> (f64, f64, f64) { (self.x().to_radians(), self.y().to_radians(), self.z()) } /// A tuple containing the first four components of the CoordinateTuple, /// with the first two converted from degrees to radians - fn xyzt_in_radians(&self) -> (f64, f64, f64, f64) { + fn xyzt_to_radians(&self) -> (f64, f64, f64, f64) { (self.x().to_radians(), self.y().to_radians(), self.z(), self.t()) } + /// Fill all elements of `self` with `value` + fn fill(&mut self, value: f64) { + for n in 0..self.dim() { + self.unchecked_set_nth(n, value); + } + } + /// Replace the n'th (0-based) element of the `CoordinateTuple` with `value`. - /// May panic if n >= [DIMENSION]. - /// See also [`set_nth()`](crate::coordinate::CoordinateTuple::nth). - fn unchecked_set_nth(&mut self, n: usize, value: f64); + /// If `n >=` [`dim()`](Self::dim()) fill the coordinate with `f64::NAN`. + /// See also [`unchecked_set_nth()`](Self::unchecked_set_nth). + fn set_nth(&mut self, n: usize, value: f64) { + if n < self.dim() { + self.unchecked_set_nth(n, value) + } else { + self.fill(f64::NAN); + } + } - /// Replace the `N` first (up to [DIMENSION]) elements of `self` with the - /// elements of `value` - fn set(&mut self, value: [f64; N]) { - for i in 0..N.min(Self::DIMENSION) { - self.unchecked_set_nth(i, value[i]) + /// Replace the two first elements of the `CoordinateTuple` with `x` and `y`. + /// If the dimension in less than 2, fill the coordinate with `f64::NAN`. + /// See also [`unchecked_set_nth()`](Self::unchecked_set_nth). + fn set_xy(&mut self, x: f64, y: f64) { + if self.dim() > 1 { + self.unchecked_set_nth(0, x); + self.unchecked_set_nth(1, y); + } else { + self.fill(f64::NAN); } } -} + /// Replace the three first elements of the `CoordinateTuple` with `x`, `y` and `z`. + /// If the dimension is less than 3, fill the coordinate with `f64::NAN`. + fn set_xyz(&mut self, x: f64, y: f64, z: f64) { + if self.dim() > 2 { + self.unchecked_set_nth(0, x); + self.unchecked_set_nth(1, y); + self.unchecked_set_nth(2, z); + } else { + self.fill(f64::NAN); + } + } + /// Replace the four first elements of the `CoordinateTuple` with `x`, `y` `z` and `t`. + /// If the dimension in less than 4, fill the coordinate with `f64::NAN`. + fn set_xyzt(&mut self, x: f64, y: f64, z: f64, t: f64) { + if self.dim() > 3 { + self.unchecked_set_nth(0, x); + self.unchecked_set_nth(1, y); + self.unchecked_set_nth(2, z); + self.unchecked_set_nth(3, t); + } else { + self.fill(f64::NAN); + } +} + /// Replace the `N` first (up to [`dim()`](Self::dim())) elements of `self` with the + /// elements of `value` + #[allow(clippy::needless_range_loop)] + fn update(&mut self, value: &[f64]) { + let n = value.len().min(self.dim()); + for i in 0..n { + self.unchecked_set_nth(i, value[i]) + } + } -// We must still implement the CoordinateTuple trait for -// the Geodesy data types Coor2D, Coor32, Coor3D, Coor4D + /// Euclidean distance between two points in the 2D plane. + /// + /// Primarily used to compute the distance between two projected points + /// in their projected plane. Typically, this distance will differ from + /// the actual distance in the real world. + /// + /// # See also: + /// + /// [`hypot3`](Self::hypot3), + /// [`distance`](crate::ellipsoid::Ellipsoid::distance) + /// + /// # Examples + /// + /// ``` + /// use geodesy::prelude::*; + /// let t = 1000 as f64; + /// let p0 = Coor2D::origin(); + /// let p1 = Coor2D::raw(t, t); + /// assert_eq!(p0.hypot2(&p1), t.hypot(t)); + /// ``` + #[must_use] + fn hypot2(&self, other: &Self) -> f64 + where Self: Sized { + let (u, v) = self.xy(); + let (x, y) = other.xy(); + (u - x).hypot(v - y) + } + + /// Euclidean distance between two points in the 3D space. + /// + /// Primarily used to compute the distance between two points in the + /// 3D cartesian space. The typical case is GNSS-observations, in which + /// case, the distance computed will reflect the actual distance + /// in the real world. + /// + /// The distance is computed in the subspace spanned by the first, + /// second and third coordinate of the `Coor3D`s + /// + /// # See also: + /// + /// [`hypot2()`](Self::hypot2), + /// [`distance`](crate::ellipsoid::Ellipsoid::distance) + /// + /// # Examples + /// + /// ``` + /// use geodesy::prelude::*; + /// let t = 1000 as f64; + /// let p0 = Coor3D::origin(); + /// let p1 = Coor3D::raw(t, t, t); + /// assert_eq!(p0.hypot3(&p1), t.hypot(t).hypot(t)); + /// ``` + #[must_use] + fn hypot3(&self, other: &Self) -> f64 + where Self: Sized { + if self.dim() < 3 { + return f64::NAN; + } + let (u, v, w) = self.xyz(); + let (x, y, z) = other.xyz(); + (u - x).hypot(v - y).hypot(w - z) + } +} // We must still implement the CoordinateTuple trait for // the Geodesy data types Coor2D, Coor32, Coor3D, Coor4D impl CoordinateTuple for Coor2D { - const DIMENSION: usize = 2; + fn dim(&self) -> usize { + 2 + } + fn unchecked_nth(&self, n: usize) -> f64 { self.0[n] } @@ -287,7 +459,10 @@ impl CoordinateTuple for Coor2D { } impl CoordinateTuple for Coor3D { - const DIMENSION: usize = 3; + fn dim(&self) -> usize { + 3 + } + fn unchecked_nth(&self, n: usize) -> f64 { self.0[n] } @@ -298,7 +473,10 @@ impl CoordinateTuple for Coor3D { } impl CoordinateTuple for Coor4D { - const DIMENSION: usize = 4; + fn dim(&self) -> usize { + 4 + } + fn unchecked_nth(&self, n: usize) -> f64 { self.0[n] } @@ -309,7 +487,10 @@ impl CoordinateTuple for Coor4D { } impl CoordinateTuple for Coor32 { - const DIMENSION: usize = 2; + fn dim(&self) -> usize { + 2 + } + fn unchecked_nth(&self, n: usize) -> f64 { self.0[n] as f64 } @@ -322,7 +503,8 @@ impl CoordinateTuple for Coor32 { // And let's also implement it for a plain 2D f64 tuple #[rustfmt::skip] impl CoordinateTuple for (f64, f64) { - const DIMENSION: usize = 2; + fn dim(&self) -> usize { 2 } + fn unchecked_nth(&self, n: usize) -> f64 { match n { 0 => self.0, diff --git a/src/coordinate/set.rs b/src/coordinate/set.rs index 69d99f7d..1ac08111 100644 --- a/src/coordinate/set.rs +++ b/src/coordinate/set.rs @@ -1,52 +1,13 @@ use super::*; -impl AngularUnits for &mut T -where - T: CoordinateSet, -{ - /// Transform the first two elements of all elements in a - /// coordinate set from degrees to radians - fn to_radians(self) -> Self { - for i in 0..self.len() { - self.set_coord(i, &self.get_coord(i).to_radians()); - } - self - } - - /// Transform the first two elements of a all elements in a - /// coordinate set from radians to degrees - #[must_use] - fn to_degrees(self) -> Self { - for i in 0..self.len() { - self.set_coord(i, &self.get_coord(i).to_degrees()); - } - self - } - - /// Transform the first two elements of a all elements in a - /// coordinate set from radians to seconds of arc. - #[must_use] - fn to_arcsec(self) -> Self { - for i in 0..self.len() { - self.set_coord(i, &self.get_coord(i).to_arcsec()); - } - self - } - - /// Transform all elements in a coordinate set from the internal - /// lon/lat/h/t-in-radians to lat/lon/h/t-in-degrees - fn to_geo(self) -> Self { - for i in 0..self.len() { - self.set_coord(i, &self.get_coord(i).to_geo()); - } - self - } -} - // ----- CoordinateSet implementations for some Coor4D containers ------------ -impl CoordinateSet for [Coor4D; N] { + +impl CoordinateSet for &mut [Coor4D] { fn len(&self) -> usize { - N + (**self).len() + } + fn dim(&self) -> usize { + 4 } fn get_coord(&self, index: usize) -> Coor4D { self[index] @@ -56,9 +17,12 @@ impl CoordinateSet for [Coor4D; N] { } } -impl CoordinateSet for &mut [Coor4D] { +impl CoordinateSet for [Coor4D; N] { fn len(&self) -> usize { - (**self).len() + N + } + fn dim(&self) -> usize { + 4 } fn get_coord(&self, index: usize) -> Coor4D { self[index] @@ -72,6 +36,9 @@ impl CoordinateSet for Vec { fn len(&self) -> usize { self.len() } + fn dim(&self) -> usize { + 4 + } fn get_coord(&self, index: usize) -> Coor4D { self[index] } @@ -81,10 +48,14 @@ impl CoordinateSet for Vec { } // ----- CoordinateSet implementations for some Coor3D containers ------------ + impl CoordinateSet for [Coor3D; N] { fn len(&self) -> usize { N } + fn dim(&self) -> usize { + 3 + } fn get_coord(&self, index: usize) -> Coor4D { Coor4D([self[index][0], self[index][1], self[index][2], f64::NAN]) } @@ -97,6 +68,9 @@ impl CoordinateSet for &mut [Coor3D] { fn len(&self) -> usize { (**self).len() } + fn dim(&self) -> usize { + 3 + } fn get_coord(&self, index: usize) -> Coor4D { Coor4D([self[index][0], self[index][1], self[index][2], f64::NAN]) } @@ -109,6 +83,9 @@ impl CoordinateSet for Vec { fn len(&self) -> usize { self.len() } + fn dim(&self) -> usize { + 3 + } fn get_coord(&self, index: usize) -> Coor4D { Coor4D([self[index][0], self[index][1], self[index][2], f64::NAN]) } @@ -141,6 +118,9 @@ impl CoordinateSet for [Coor2D; N] { fn len(&self) -> usize { N } + fn dim(&self) -> usize { + 2 + } fn get_coord(&self, index: usize) -> Coor4D { Coor4D([self[index][0], self[index][1], 0., f64::NAN]) } @@ -153,6 +133,9 @@ impl CoordinateSet for &mut [Coor2D] { fn len(&self) -> usize { (**self).len() } + fn dim(&self) -> usize { + 2 + } fn get_coord(&self, index: usize) -> Coor4D { Coor4D([self[index][0], self[index][1], 0.0, f64::NAN]) } @@ -165,6 +148,9 @@ impl CoordinateSet for Vec { fn len(&self) -> usize { self.len() } + fn dim(&self) -> usize { + 2 + } fn get_coord(&self, index: usize) -> Coor4D { Coor4D([self[index][0], self[index][1], 0., f64::NAN]) } @@ -183,6 +169,9 @@ where fn len(&self) -> usize { self.0.len() } + fn dim(&self) -> usize { + 4 + } fn get_coord(&self, index: usize) -> Coor4D { let c = self.0.get_coord(index); Coor4D([c[0], c[1], self.1, self.2]) @@ -202,6 +191,9 @@ where fn len(&self) -> usize { self.0.len() } + fn dim(&self) -> usize { + 4 + } fn get_coord(&self, index: usize) -> Coor4D { let c = self.0.get_coord(index); Coor4D([c[0], c[1], c[2], self.1]) @@ -213,36 +205,15 @@ where // ----- CoordinateSet implementations for some Coor32 containers ------------ -impl CoordinateSet for [Coor32; N] { - fn len(&self) -> usize { - N - } - fn get_coord(&self, index: usize) -> Coor4D { - Coor4D([self[index][0] as f64, self[index][1] as f64, 0., f64::NAN]) - } - fn set_coord(&mut self, index: usize, value: &Coor4D) { - self[index] = Coor32::raw(value[0], value[1]); - } -} - impl CoordinateSet for &mut [Coor32] { fn len(&self) -> usize { (**self).len() } - fn get_coord(&self, index: usize) -> Coor4D { - Coor4D([self[index][0] as f64, self[index][1] as f64, 0.0, f64::NAN]) - } - fn set_coord(&mut self, index: usize, value: &Coor4D) { - self[index] = Coor32::raw(value[0], value[1]); - } -} - -impl CoordinateSet for Vec { - fn len(&self) -> usize { - self.len() + fn dim(&self) -> usize { + 2 } fn get_coord(&self, index: usize) -> Coor4D { - Coor4D([self[index][0] as f64, self[index][1] as f64, 0., f64::NAN]) + Coor4D([self[index][0] as f64, self[index][1] as f64, 0.0, f64::NAN]) } fn set_coord(&mut self, index: usize, value: &Coor4D) { self[index] = Coor32::raw(value[0], value[1]); @@ -321,20 +292,24 @@ mod tests { // Test the "AngularUnits" conversion trait #[test] fn angular() { - let mut operands = some_basic_coor2dinates(); + let operands = some_basic_coor2dinates(); let cph = operands.get_coord(0); // Note the different usage patterns when using the AngularUnits trait with // a Coor4D and a CoordinateSet: For the latter, the blanket implementation // is for an `&mut T where T: CoordinateSet`, and we just mutate the contents // in situ. For the former, we return a newly computed `Coor4D`. - operands.to_radians(); let cph = cph.to_radians(); - assert_eq!(cph[0], operands.get_coord(0)[0]); - assert_eq!(cph[1], operands.get_coord(0)[1]); - - operands.to_arcsec(); - assert_eq!(cph[0].to_degrees() * 3600., operands.get_coord(0)[0]); - assert_eq!(cph[1].to_degrees() * 3600., operands.get_coord(0)[1]); + assert_eq!(cph[0], operands.get_coord(0).to_radians()[0]); + assert_eq!(cph[1], operands.get_coord(0).to_radians()[1]); + + assert_eq!( + cph[0].to_degrees() * 3600., + operands.get_coord(0).to_radians().to_arcsec()[0] + ); + assert_eq!( + cph[1].to_degrees() * 3600., + operands.get_coord(0).to_radians().to_arcsec()[1] + ); } } diff --git a/src/ellipsoid/geodesics.rs b/src/ellipsoid/geodesics.rs index 558472a3..e0b3b1f7 100644 --- a/src/ellipsoid/geodesics.rs +++ b/src/ellipsoid/geodesics.rs @@ -16,7 +16,12 @@ impl Ellipsoid { /// Federico Dolce and Michael Kirk, provides a Rust implementation of Karney's algorithm. #[must_use] #[allow(non_snake_case)] - pub fn geodesic_fwd(&self, from: &C, azimuth: f64, distance: f64) -> Coor4D { + pub fn geodesic_fwd( + &self, + from: &C, + azimuth: f64, + distance: f64, + ) -> Coor4D { // Coordinates of the point of origin, P1 let (L1, B1) = from.xy(); @@ -238,9 +243,9 @@ mod tests { assert!((d[2] - 2365723.367715).abs() < 1e-4); // And the other way round... - let b = ellps.geodesic_fwd(&p1, d[0], d[2]).to_degrees(); - assert!((b[0] - p2[0].to_degrees()).abs() < 1e-9); - assert!((b[1] - p2[1].to_degrees()).abs() < 1e-9); + let b = ellps.geodesic_fwd(&p1, d[0], d[2]); + assert!((b[0].to_degrees() - p2[0].to_degrees()).abs() < 1e-9); + assert!((b[1].to_degrees() - p2[1].to_degrees()).abs() < 1e-9); Ok(()) } } diff --git a/src/grid/mod.rs b/src/grid/mod.rs index 728c3bbc..fea9756c 100644 --- a/src/grid/mod.rs +++ b/src/grid/mod.rs @@ -340,6 +340,7 @@ pub fn grids_at(grids: &[Arc], coord: &Coor4D, use_null_grid: bool) -> #[cfg(test)] mod tests { use super::*; + use crate::coordinate::AngularUnits; // lat_n, lat_s, lon_w, lon_e, dlat, dlon const HEADER: [f64; 6] = [58., 54., 8., 16., -1., 1.]; diff --git a/src/inner_op/gridshift.rs b/src/inner_op/gridshift.rs index 9f67d2cc..ba8efc71 100644 --- a/src/inner_op/gridshift.rs +++ b/src/inner_op/gridshift.rs @@ -148,6 +148,7 @@ pub fn new(parameters: &RawParameters, ctx: &dyn Context) -> Result { #[cfg(test)] mod tests { use super::*; + use crate::coordinate::AngularUnits; #[test] fn gridshift() -> Result<(), Error> { diff --git a/src/op/mod.rs b/src/op/mod.rs index 5dfbe9ba..b08af4d4 100644 --- a/src/op/mod.rs +++ b/src/op/mod.rs @@ -381,11 +381,11 @@ mod tests { // Now instantiating the macro with ham = 2 let op = ctx.op("helmert:ham ham=2")?; ctx.apply(op, Fwd, &mut data)?; - assert_eq!(data[0][0], 57.); - assert_eq!(data[1][0], 61.); + assert_eq!(data[0].x(), 57.); + assert_eq!(data[1].x(), 61.); ctx.apply(op, Inv, &mut data)?; - assert_eq!(data[0][0], 55.); - assert_eq!(data[1][0], 59.); + assert_eq!(data[0].x(), 55.); + assert_eq!(data[1].x(), 59.); Ok(()) } From 57485b30f9fe0a7ca417696fc097eed045fca364 Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Thu, 4 Apr 2024 23:20:28 +0200 Subject: [PATCH 08/10] parcsec->arcsec in Coor4D --- src/coordinate/coor4d.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coordinate/coor4d.rs b/src/coordinate/coor4d.rs index 4b88b7c2..94a11548 100644 --- a/src/coordinate/coor4d.rs +++ b/src/coordinate/coor4d.rs @@ -95,7 +95,7 @@ impl Coor4D { /// A `Coor4D` from longitude/latitude/height/time, with the angular input in seconds /// of arc. Mostly for handling grid shift elements. #[must_use] - pub fn parcsec(longitude: f64, latitude: f64, height: f64, time: f64) -> Coor4D { + pub fn arcsec(longitude: f64, latitude: f64, height: f64, time: f64) -> Coor4D { Coor4D([ longitude.to_radians() / 3600., latitude.to_radians() / 3600., From 94e6fea79cc02c487a080d24e1957398429d015c Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Fri, 5 Apr 2024 09:11:28 +0200 Subject: [PATCH 09/10] Prefix unchecked to postfix --- src/coordinate/mod.rs | 70 +++++++++++++++++++++---------------------- 1 file changed, 35 insertions(+), 35 deletions(-) diff --git a/src/coordinate/mod.rs b/src/coordinate/mod.rs index f299d049..5e403e58 100644 --- a/src/coordinate/mod.rs +++ b/src/coordinate/mod.rs @@ -189,8 +189,8 @@ pub trait CoordinateSet: CoordinateMetadata { /// Geodesy operators. /// /// All accessors have default implementations, except the 3 methods -/// [`unchecked_nth()`](Self::unchecked_nth()), -/// [`unchecked_set_nth()`](Self::unchecked_set_nth) and +/// [`nth_unchecked()`](Self::nth_unchecked()), +/// [`set_nth_unchecked()`](Self::set_nth_unchecked) and /// [`dim()`](Self::dim()), /// which must be provided by the implementer. /// @@ -201,46 +201,46 @@ pub trait CoordinateTuple { /// Access the n'th (0-based) element of the CoordinateTuple. /// May panic if n >= DIMENSION. /// See also [`nth()`](Self::nth). - fn unchecked_nth(&self, n: usize) -> f64; + fn nth_unchecked(&self, n: usize) -> f64; /// Replace the n'th (0-based) element of the `CoordinateTuple` with `value`. /// May panic if `n >=` [`dim()`](Self::dim()). /// See also [`set_nth()`](Self::set_nth). - fn unchecked_set_nth(&mut self, n: usize, value: f64); + fn set_nth_unchecked(&mut self, n: usize, value: f64); /// Native dimension of the coordinate tuple fn dim(&self) -> usize; /// Access the n'th (0-based) element of the CoordinateTuple. /// Returns NaN if `n >= DIMENSION`. - /// See also [`unchecked_nth()`](Self::unchecked_nth). + /// See also [`nth()`](Self::nth_unchecked). fn nth(&self, n: usize) -> f64 { - if n < self.dim() { self.unchecked_nth(n) } else {f64::NAN} + if n < self.dim() { self.nth_unchecked(n) } else {f64::NAN} } - // Note: We use unchecked_nth and explicitly check for dimension in + // Note: We use nth_unchecked and explicitly check for dimension in // y(), z() and t(), rather than leaving the check to nth(...). // This is because the checks in these cases are constant expressions, and // hence can be eliminated by the compiler in the concrete implementations. /// Pragmatically named accessor for the first element of the CoordinateTuple. fn x(&self) -> f64 { - self.unchecked_nth(0) + self.nth_unchecked(0) } /// Pragmatically named accessor for the second element of the CoordinateTuple. fn y(&self) -> f64 { - if self.dim() > 1 { self.unchecked_nth(1) } else {f64::NAN} + if self.dim() > 1 { self.nth_unchecked(1) } else {f64::NAN} } /// Pragmatically named accessor for the third element of the CoordinateTuple. fn z(&self) -> f64 { - if self.dim() > 2 { self.unchecked_nth(2) } else {f64::NAN} + if self.dim() > 2 { self.nth_unchecked(2) } else {f64::NAN} } /// Pragmatically named accessor for the fourth element of the CoordinateTuple. fn t(&self) -> f64 { - if self.dim() > 3 { self.unchecked_nth(3) } else {f64::NAN} + if self.dim() > 3 { self.nth_unchecked(3) } else {f64::NAN} } /// A tuple containing the first two components of the CoordinateTuple. @@ -315,16 +315,16 @@ pub trait CoordinateTuple { /// Fill all elements of `self` with `value` fn fill(&mut self, value: f64) { for n in 0..self.dim() { - self.unchecked_set_nth(n, value); + self.set_nth_unchecked(n, value); } } /// Replace the n'th (0-based) element of the `CoordinateTuple` with `value`. /// If `n >=` [`dim()`](Self::dim()) fill the coordinate with `f64::NAN`. - /// See also [`unchecked_set_nth()`](Self::unchecked_set_nth). + /// See also [`set_nth_unchecked()`](Self::set_nth_unchecked). fn set_nth(&mut self, n: usize, value: f64) { if n < self.dim() { - self.unchecked_set_nth(n, value) + self.set_nth_unchecked(n, value) } else { self.fill(f64::NAN); } @@ -332,11 +332,11 @@ pub trait CoordinateTuple { /// Replace the two first elements of the `CoordinateTuple` with `x` and `y`. /// If the dimension in less than 2, fill the coordinate with `f64::NAN`. - /// See also [`unchecked_set_nth()`](Self::unchecked_set_nth). + /// See also [`set_nth_unchecked()`](Self::set_nth_unchecked). fn set_xy(&mut self, x: f64, y: f64) { if self.dim() > 1 { - self.unchecked_set_nth(0, x); - self.unchecked_set_nth(1, y); + self.set_nth_unchecked(0, x); + self.set_nth_unchecked(1, y); } else { self.fill(f64::NAN); } @@ -346,9 +346,9 @@ pub trait CoordinateTuple { /// If the dimension is less than 3, fill the coordinate with `f64::NAN`. fn set_xyz(&mut self, x: f64, y: f64, z: f64) { if self.dim() > 2 { - self.unchecked_set_nth(0, x); - self.unchecked_set_nth(1, y); - self.unchecked_set_nth(2, z); + self.set_nth_unchecked(0, x); + self.set_nth_unchecked(1, y); + self.set_nth_unchecked(2, z); } else { self.fill(f64::NAN); } @@ -358,10 +358,10 @@ pub trait CoordinateTuple { /// If the dimension in less than 4, fill the coordinate with `f64::NAN`. fn set_xyzt(&mut self, x: f64, y: f64, z: f64, t: f64) { if self.dim() > 3 { - self.unchecked_set_nth(0, x); - self.unchecked_set_nth(1, y); - self.unchecked_set_nth(2, z); - self.unchecked_set_nth(3, t); + self.set_nth_unchecked(0, x); + self.set_nth_unchecked(1, y); + self.set_nth_unchecked(2, z); + self.set_nth_unchecked(3, t); } else { self.fill(f64::NAN); } @@ -373,7 +373,7 @@ pub trait CoordinateTuple { fn update(&mut self, value: &[f64]) { let n = value.len().min(self.dim()); for i in 0..n { - self.unchecked_set_nth(i, value[i]) + self.set_nth_unchecked(i, value[i]) } } @@ -449,11 +449,11 @@ impl CoordinateTuple for Coor2D { 2 } - fn unchecked_nth(&self, n: usize) -> f64 { + fn nth_unchecked(&self, n: usize) -> f64 { self.0[n] } - fn unchecked_set_nth(&mut self, n: usize, value: f64) { + fn set_nth_unchecked(&mut self, n: usize, value: f64) { self.0[n] = value; } } @@ -463,11 +463,11 @@ impl CoordinateTuple for Coor3D { 3 } - fn unchecked_nth(&self, n: usize) -> f64 { + fn nth_unchecked(&self, n: usize) -> f64 { self.0[n] } - fn unchecked_set_nth(&mut self, n: usize, value: f64) { + fn set_nth_unchecked(&mut self, n: usize, value: f64) { self.0[n] = value; } } @@ -477,11 +477,11 @@ impl CoordinateTuple for Coor4D { 4 } - fn unchecked_nth(&self, n: usize) -> f64 { + fn nth_unchecked(&self, n: usize) -> f64 { self.0[n] } - fn unchecked_set_nth(&mut self, n: usize, value: f64) { + fn set_nth_unchecked(&mut self, n: usize, value: f64) { self.0[n] = value; } } @@ -491,11 +491,11 @@ impl CoordinateTuple for Coor32 { 2 } - fn unchecked_nth(&self, n: usize) -> f64 { + fn nth_unchecked(&self, n: usize) -> f64 { self.0[n] as f64 } - fn unchecked_set_nth(&mut self, n: usize, value: f64) { + fn set_nth_unchecked(&mut self, n: usize, value: f64) { self.0[n] = value as f32; } } @@ -505,7 +505,7 @@ impl CoordinateTuple for Coor32 { impl CoordinateTuple for (f64, f64) { fn dim(&self) -> usize { 2 } - fn unchecked_nth(&self, n: usize) -> f64 { + fn nth_unchecked(&self, n: usize) -> f64 { match n { 0 => self.0, 1 => self.1, @@ -513,7 +513,7 @@ impl CoordinateTuple for (f64, f64) { } } - fn unchecked_set_nth(&mut self, n: usize, value: f64) { + fn set_nth_unchecked(&mut self, n: usize, value: f64) { match n { 0 => self.0 = value, 1 => self.1 = value, From ef087ccc4d74cac90041b55ee523d77fe4b083c2 Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Fri, 5 Apr 2024 12:01:38 +0200 Subject: [PATCH 10/10] xy/set_xy for CoordinateSet. Tests via tmerc/merc --- src/coordinate/mod.rs | 21 +++++++++++++++ src/coordinate/set.rs | 60 +++++++++++++++++++++++++++++++++++++++++++ src/inner_op/merc.rs | 26 +++++++++---------- src/inner_op/tmerc.rs | 53 +++++++++++++++++++------------------- 4 files changed, 120 insertions(+), 40 deletions(-) diff --git a/src/coordinate/mod.rs b/src/coordinate/mod.rs index 5e403e58..ae2d2b28 100644 --- a/src/coordinate/mod.rs +++ b/src/coordinate/mod.rs @@ -164,6 +164,27 @@ pub trait CoordinateSet: CoordinateMetadata { self.len() == 0 } + /// Replace the two first elements of the `index`th `CoordinateTuple` with `x` and `y`. + /// If the dimension in less than 2, fill the coordinate with `f64::NAN`. + /// For efficiency, consider implemented a specific implementation, when implementing + /// the CoordinateSet trait for a concrete data type: The provided version is not + /// very efficient. + fn set_xy(&mut self, index: usize, x: f64, y: f64) { + let mut coord = self.get_coord(index); + coord.set_nth_unchecked(0, x); + coord.set_nth_unchecked(1, y); + self.set_coord(index, &coord); + } + + /// Replace the two first elements of the `index`th `CoordinateTuple` with `x` and `y`. + /// If the dimension in less than 2, fill the coordinate with `f64::NAN`. + /// For efficiency, consider implemented a specific implementation, when implementing + /// the CoordinateSet trait for a concrete data type: The provided version is not + /// very efficient. + fn xy(&self, index: usize) -> (f64, f64) { + self.get_coord(index).xy() + } + /// Set all coordinate tuples in the set to NaN fn stomp(&mut self) { let nanny = Coor4D::nan(); diff --git a/src/coordinate/set.rs b/src/coordinate/set.rs index 1ac08111..b4e56560 100644 --- a/src/coordinate/set.rs +++ b/src/coordinate/set.rs @@ -15,6 +15,12 @@ impl CoordinateSet for &mut [Coor4D] { fn set_coord(&mut self, index: usize, value: &Coor4D) { self[index] = *value; } + fn xy(&self, index: usize) -> (f64, f64) { + self[index].xy() + } + fn set_xy(&mut self, index: usize, x: f64, y: f64) { + self[index].set_xy(x, y); + } } impl CoordinateSet for [Coor4D; N] { @@ -30,6 +36,12 @@ impl CoordinateSet for [Coor4D; N] { fn set_coord(&mut self, index: usize, value: &Coor4D) { self[index] = *value; } + fn xy(&self, index: usize) -> (f64, f64) { + self[index].xy() + } + fn set_xy(&mut self, index: usize, x: f64, y: f64) { + self[index].set_xy(x, y); + } } impl CoordinateSet for Vec { @@ -45,6 +57,12 @@ impl CoordinateSet for Vec { fn set_coord(&mut self, index: usize, value: &Coor4D) { self[index] = *value; } + fn xy(&self, index: usize) -> (f64, f64) { + self[index].xy() + } + fn set_xy(&mut self, index: usize, x: f64, y: f64) { + self[index].set_xy(x, y); + } } // ----- CoordinateSet implementations for some Coor3D containers ------------ @@ -62,6 +80,12 @@ impl CoordinateSet for [Coor3D; N] { fn set_coord(&mut self, index: usize, value: &Coor4D) { self[index] = Coor3D([value[0], value[1], value[2]]); } + fn xy(&self, index: usize) -> (f64, f64) { + self[index].xy() + } + fn set_xy(&mut self, index: usize, x: f64, y: f64) { + self[index].set_xy(x, y); + } } impl CoordinateSet for &mut [Coor3D] { @@ -77,6 +101,12 @@ impl CoordinateSet for &mut [Coor3D] { fn set_coord(&mut self, index: usize, value: &Coor4D) { self[index] = Coor3D([value[0], value[1], value[2]]); } + fn xy(&self, index: usize) -> (f64, f64) { + self[index].xy() + } + fn set_xy(&mut self, index: usize, x: f64, y: f64) { + self[index].set_xy(x, y); + } } impl CoordinateSet for Vec { @@ -92,6 +122,12 @@ impl CoordinateSet for Vec { fn set_coord(&mut self, index: usize, value: &Coor4D) { self[index] = Coor3D([value[0], value[1], value[2]]); } + fn xy(&self, index: usize) -> (f64, f64) { + self[index].xy() + } + fn set_xy(&mut self, index: usize, x: f64, y: f64) { + self[index].set_xy(x, y); + } } // ----- CoordinateSet implementations for some Coor2D containers ------------ @@ -127,6 +163,12 @@ impl CoordinateSet for [Coor2D; N] { fn set_coord(&mut self, index: usize, value: &Coor4D) { self[index] = Coor2D([value[0], value[1]]); } + fn xy(&self, index: usize) -> (f64, f64) { + self[index].xy() + } + fn set_xy(&mut self, index: usize, x: f64, y: f64) { + self[index].set_xy(x, y); + } } impl CoordinateSet for &mut [Coor2D] { @@ -142,6 +184,12 @@ impl CoordinateSet for &mut [Coor2D] { fn set_coord(&mut self, index: usize, value: &Coor4D) { self[index] = Coor2D([value[0], value[1]]); } + fn xy(&self, index: usize) -> (f64, f64) { + self[index].xy() + } + fn set_xy(&mut self, index: usize, x: f64, y: f64) { + self[index].set_xy(x, y); + } } impl CoordinateSet for Vec { @@ -157,6 +205,12 @@ impl CoordinateSet for Vec { fn set_coord(&mut self, index: usize, value: &Coor4D) { self[index] = Coor2D([value[0], value[1]]); } + fn xy(&self, index: usize) -> (f64, f64) { + self[index].xy() + } + fn set_xy(&mut self, index: usize, x: f64, y: f64) { + self[index].set_xy(x, y); + } } /// User defined values for third and fourth coordinate dimension. @@ -218,6 +272,12 @@ impl CoordinateSet for &mut [Coor32] { fn set_coord(&mut self, index: usize, value: &Coor4D) { self[index] = Coor32::raw(value[0], value[1]); } + fn xy(&self, index: usize) -> (f64, f64) { + self[index].xy() + } + fn set_xy(&mut self, index: usize, x: f64, y: f64) { + self[index].set_xy(x, y); + } } // ----- Implementations: Coordinate Metadata --------------------------------- diff --git a/src/inner_op/merc.rs b/src/inner_op/merc.rs index a9a35fe6..289e7a40 100644 --- a/src/inner_op/merc.rs +++ b/src/inner_op/merc.rs @@ -15,15 +15,13 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { let mut successes = 0_usize; let length = operands.len(); for i in 0..length { - let mut coord = operands.get_coord(i); - // Easting - coord[0] = (coord[0] - lon_0) * k_0 * a - x_0; + let (lon, lat) = operands.xy(i); - // Northing - let lat = coord[1] + lat_0; - coord[1] = a * k_0 * ellps.latitude_geographic_to_isometric(lat) - y_0; + let easting = (lon - lon_0) * k_0 * a - x_0; + let isometric = ellps.latitude_geographic_to_isometric(lat + lat_0); + let northing = a * k_0 * isometric - y_0; - operands.set_coord(i, &coord); + operands.set_xy(i, easting, northing); successes += 1; } @@ -44,17 +42,17 @@ fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { let mut successes = 0_usize; let length = operands.len(); for i in 0..length { - let mut coord = operands.get_coord(i); + let (mut x, mut y) = operands.xy(i); // Easting -> Longitude - let x = coord[0] + x_0; - coord[0] = x / (a * k_0) - lon_0; + x += x_0; + let lon = x / (a * k_0) - lon_0; // Northing -> Latitude - let y = coord[1] + y_0; + y += y_0; let psi = y / (a * k_0); - coord[1] = ellps.latitude_isometric_to_geographic(psi) - lat_0; - operands.set_coord(i, &coord); + let lat = ellps.latitude_isometric_to_geographic(psi) - lat_0; + operands.set_xy(i, lon, lat); successes += 1; } @@ -160,9 +158,9 @@ mod tests { let definition = "merc lat_ts=56"; let op = ctx.op(definition)?; - // Validation value from PROJ: echo 12 55 0 0 | cct -d18 +proj=merc +lat_ts=56 let geo = [Coor4D::geo(55., 12., 0., 0.)]; + // Validation value from PROJ: echo 12 55 0 0 | cct -d18 +proj=merc +lat_ts=56 let projected = [Coor4D::raw( 748_713.257_925_886_8, 4_106_573.862_841_270_4, diff --git a/src/inner_op/tmerc.rs b/src/inner_op/tmerc.rs index b163b6b3..ac4f79ad 100644 --- a/src/inner_op/tmerc.rs +++ b/src/inner_op/tmerc.rs @@ -29,14 +29,15 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { let range = 0..operands.len(); let mut successes = 0_usize; for i in range { - let mut coord = operands.get_coord(i); + //let mut coord = operands.get_coord(i); + let (lon, lat) = operands.xy(i); // --- 1. Geographical -> Conformal latitude, rotated longitude // The conformal latitude - let lat = ellps.latitude_geographic_to_conformal(coord[1], conformal); + let lat = ellps.latitude_geographic_to_conformal(lat, conformal); // The longitude as reckoned from the central meridian - let lon = coord[0] - lon_0; + let lon = lon - lon_0; // --- 2. Conformal LAT, LNG -> complex spherical LAT @@ -48,7 +49,7 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { // --- 3. Complex spherical N, E -> ellipsoidal normalized N, E // Some numerical optimizations from PROJ modifications by Even Rouault, - let inv_denom_tan_lon = 1. / sin_lat.hypot(cos_lat_lon); + let inv_denom_tan_lon = sin_lat.hypot(cos_lat_lon).recip(); let tan_lon = sin_lon * cos_lat * inv_denom_tan_lon; // Inverse Gudermannian, using the precomputed tan(lon) let mut lon = tan_lon.asinh(); @@ -74,17 +75,18 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { // Don't wanna play if we're too far from the center meridian if lon.abs() > 2.623395162778 { - coord[0] = f64::NAN; - coord[1] = f64::NAN; + operands.set_xy(i, f64::NAN, f64::NAN); continue; } // --- 4. ellipsoidal normalized N, E -> metric N, E - coord[0] = qs * lon + x_0; // Easting - coord[1] = qs * lat + zb; // Northing + let easting = qs * lon + x_0; // Easting + let northing = qs * lat + zb; // Northing + + // Done! + operands.set_xy(i, easting, northing); successes += 1; - operands.set_coord(i, &coord); } successes @@ -118,17 +120,16 @@ fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { let range = 0..operands.len(); let mut successes = 0_usize; for i in range { - let mut coord = operands.get_coord(i); + let (x, y) = operands.xy(i); // --- 1. Normalize N, E - let mut lon = (coord[0] - x_0) / qs; - let mut lat = (coord[1] - zb) / qs; + let mut lon = (x - x_0) / qs; + let mut lat = (y - zb) / qs; // Don't wanna play if we're too far from the center meridian if lon.abs() > 2.623395162778 { - coord[0] = f64::NAN; - coord[1] = f64::NAN; + operands.set_xy(i, f64::NAN, f64::NAN); continue; } @@ -151,10 +152,10 @@ fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { let lon = angular::normalize_symmetric(lon + lon_0); let lat = ellps.latitude_conformal_to_geographic(lat, conformal); - (coord[0], coord[1]) = (lon, lat); + // Done! + operands.set_xy(i, lon, lat); successes += 1; - operands.set_coord(i, &coord); } successes @@ -374,8 +375,8 @@ mod tests { let geo = [ Coor2D::geo( 55., 12.), Coor2D::geo(-55., 12.), - Coor2D::geo( 55., -6.), - Coor2D::geo(-55., -6.) + Coor2D::geo( 55., -6.), + Coor2D::geo(-55., -6.) ]; #[rustfmt::skip] @@ -407,18 +408,18 @@ mod tests { #[rustfmt::skip] let geo = [ - Coor4D::geo( 55., 12., 0., 0.), - Coor4D::geo(-55., 12., 0., 0.), - Coor4D::geo( 55., -6., 0., 0.), - Coor4D::geo(-55., -6., 0., 0.) + Coor2D::geo( 55., 12.), + Coor2D::geo(-55., 12.), + Coor2D::geo( 55., -6.,), + Coor2D::geo(-55., -6.,) ]; #[rustfmt::skip] let projected = [ - Coor4D::raw( 691_875.632_139_661, 1e7+6_098_907.825_005_012, 0., 0.), - Coor4D::raw( 691_875.632_139_661, 1e7-6_098_907.825_005_012, 0., 0.), - Coor4D::raw(-455_673.814_189_040, 1e7+6_198_246.671_090_279, 0., 0.), - Coor4D::raw(-455_673.814_189_040, 1e7-6_198_246.671_090_279, 0., 0.) + Coor2D::raw( 691_875.632_139_661, 1e7+6_098_907.825_005_012), + Coor2D::raw( 691_875.632_139_661, 1e7-6_098_907.825_005_012), + Coor2D::raw(-455_673.814_189_040, 1e7+6_198_246.671_090_279), + Coor2D::raw(-455_673.814_189_040, 1e7-6_198_246.671_090_279) ]; let mut operands = geo;