diff --git a/examples/01-geometric_geodesy.rs b/examples/01-geometric_geodesy.rs index 7b248ae..30c57f5 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 63f9dfd..1f4fe03 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 3f0ed8e..0923c6d 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 cf2398c..bc73566 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 3cbc28f..df7486b 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 f696926..2661a6e 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 25b327c..908c6d2 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 c298bb1..33048ed 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 9f33213..9206a4b 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 3e9f50c..2b250ff 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 ffb12c5..d542a54 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 f926193..978b823 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(())