From 92c7903d7cc0f9e15b6471b2a5583ef5fd086a35 Mon Sep 17 00:00:00 2001 From: Daniel Pettersson Date: Tue, 19 Jan 2021 18:07:17 +0100 Subject: [PATCH 1/5] Add GeodesicIntermediate algorithm --- geo/src/algorithm/geodesic_intermediate.rs | 138 +++++++++++++++++++++ geo/src/algorithm/mod.rs | 2 + geo/src/lib.rs | 1 + 3 files changed, 141 insertions(+) create mode 100644 geo/src/algorithm/geodesic_intermediate.rs diff --git a/geo/src/algorithm/geodesic_intermediate.rs b/geo/src/algorithm/geodesic_intermediate.rs new file mode 100644 index 000000000..af1c3b3f7 --- /dev/null +++ b/geo/src/algorithm/geodesic_intermediate.rs @@ -0,0 +1,138 @@ +use crate::{CoordFloat, Point}; +use geographiclib_rs::{DirectGeodesic, Geodesic, InverseGeodesic}; + +/// Returns a new Point along a route between two existing points on an ellipsoidal model of the earth + +pub trait GeodesicIntermediate { + /// Returns a new Point along a route between two existing points on an ellipsoidal model of the earth + /// + /// # Examples + /// + /// ``` + /// # #[macro_use] extern crate approx; + /// # + /// use geo::algorithm::geodesic_intermediate::GeodesicIntermediate; + /// use geo::Point; + /// + /// let p1 = Point::::new(10.0, 20.0); + /// let p2 = Point::::new(125.0, 25.0); + /// let i20 = p1.geodesic_intermediate(&p2, 0.2); + /// let i50 = p1.geodesic_intermediate(&p2, 0.5); + /// let i80 = p1.geodesic_intermediate(&p2, 0.8); + /// let i20_should = Point::new(29.842907, 29.951445); + /// let i50_should = Point::new(65.879360, 37.722253); + /// let i80_should = Point::new(103.556796, 33.506196); + /// assert_relative_eq!(i20.x(), i20_should.x(), epsilon = 1.0e-6); + /// assert_relative_eq!(i20.y(), i20_should.y(), epsilon = 1.0e-6); + /// assert_relative_eq!(i50.x(), i50_should.x(), epsilon = 1.0e-6); + /// assert_relative_eq!(i50.y(), i50_should.y(), epsilon = 1.0e-6); + /// assert_relative_eq!(i80.x(), i80_should.x(), epsilon = 1.0e-6); + /// assert_relative_eq!(i80.y(), i80_should.y(), epsilon = 1.0e-6); + /// ``` + + fn geodesic_intermediate(&self, other: &Point, f: T) -> Point; + fn geodesic_intermediate_fill( + &self, + other: &Point, + max_dist: T, + include_ends: bool, + ) -> Vec>; +} + +impl GeodesicIntermediate for Point { + fn geodesic_intermediate(&self, other: &Point, f: f64) -> Point { + let g = Geodesic::wgs84(); + let (total_distance, azi1, _azi2, _a12) = + g.inverse(self.lat(), self.lng(), other.lat(), other.lng()); + let distance = total_distance * f; + let (lat2, lon2) = g.direct(self.lat(), self.lng(), azi1, distance); + + Point::new(lon2, lat2) + } + + fn geodesic_intermediate_fill( + &self, + other: &Point, + max_dist: f64, + include_ends: bool, + ) -> Vec> { + let g = Geodesic::wgs84(); + let (total_distance, azi1, _azi2, _a12) = + g.inverse(self.lat(), self.lng(), other.lat(), other.lng()); + + if total_distance <= max_dist { + if include_ends { + return vec![*self, *other]; + } else { + return vec![]; + } + } + + let number_of_points = (total_distance / max_dist).ceil(); + let interval = 1.0 / number_of_points; + + let mut current_step = interval; + let mut points = if include_ends { vec![*self] } else { vec![] }; + + while current_step < 1.0 { + let (lat2, lon2) = + g.direct(self.lat(), self.lng(), azi1, total_distance * current_step); + let point = Point::new(lon2, lat2); + points.push(point); + current_step = current_step + interval; + } + + if include_ends { + points.push(*other); + } + + points + } +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + #[test] + fn f_is_zero_or_one_test() { + let p1 = Point::::new(10.0, 20.0); + let p2 = Point::::new(15.0, 25.0); + let i0 = p1.geodesic_intermediate(&p2, 0.0); + let i100 = p1.geodesic_intermediate(&p2, 1.0); + assert_relative_eq!(i0.x(), p1.x(), epsilon = 1.0e-6); + assert_relative_eq!(i0.y(), p1.y(), epsilon = 1.0e-6); + assert_relative_eq!(i100.x(), p2.x(), epsilon = 1.0e-6); + assert_relative_eq!(i100.y(), p2.y(), epsilon = 1.0e-6); + } + + #[test] + fn various_f_values_test() { + let p1 = Point::::new(10.0, 20.0); + let p2 = Point::::new(125.0, 25.0); + let i20 = p1.geodesic_intermediate(&p2, 0.2); + let i50 = p1.geodesic_intermediate(&p2, 0.5); + let i80 = p1.geodesic_intermediate(&p2, 0.8); + let i20_should = Point::new(29.842907, 29.951445); + let i50_should = Point::new(65.879360, 37.722253); + let i80_should = Point::new(103.556796, 33.506196); + assert_relative_eq!(i20.x(), i20_should.x(), epsilon = 1.0e-6); + assert_relative_eq!(i20.y(), i20_should.y(), epsilon = 1.0e-6); + assert_relative_eq!(i50.x(), i50_should.x(), epsilon = 1.0e-6); + assert_relative_eq!(i50.y(), i50_should.y(), epsilon = 1.0e-6); + assert_relative_eq!(i80.x(), i80_should.x(), epsilon = 1.0e-6); + assert_relative_eq!(i80.y(), i80_should.y(), epsilon = 1.0e-6); + } + + #[test] + fn should_add_i50_test() { + let p1 = Point::::new(30.0, 40.0); + let p2 = Point::::new(40.0, 50.0); + let max_dist = 1000000.0; // meters + let include_ends = true; + let i50 = p1.clone().geodesic_intermediate(&p2, 0.5); + let route = p1.geodesic_intermediate_fill(&p2, max_dist, include_ends); + assert_eq!(route, vec![p1, i50, p2]); + } +} diff --git a/geo/src/algorithm/mod.rs b/geo/src/algorithm/mod.rs index 6612a5dc4..0f326004a 100644 --- a/geo/src/algorithm/mod.rs +++ b/geo/src/algorithm/mod.rs @@ -37,6 +37,8 @@ pub mod frechet_distance; pub mod geodesic_distance; /// Calculate the Geodesic length of a line. pub mod geodesic_length; +/// Calculate a new `Point` lying on a Geodesic arc between two `Point`s. +pub mod geodesic_intermediate; /// Calculate a destination `Point`, given a distance and a bearing. pub mod haversine_destination; /// Calculate the Haversine distance between two `Geometries`. diff --git a/geo/src/lib.rs b/geo/src/lib.rs index 69c198c37..9c52b25b2 100644 --- a/geo/src/lib.rs +++ b/geo/src/lib.rs @@ -105,6 +105,7 @@ pub mod prelude { pub use crate::algorithm::frechet_distance::FrechetDistance; pub use crate::algorithm::geodesic_distance::GeodesicDistance; pub use crate::algorithm::geodesic_length::GeodesicLength; + pub use crate::algorithm::geodesic_intermediate::GeodesicIntermediate; pub use crate::algorithm::haversine_destination::HaversineDestination; pub use crate::algorithm::haversine_distance::HaversineDistance; pub use crate::algorithm::haversine_intermediate::HaversineIntermediate; From 7199f110582d65c76e95284cb1ea48e414b3c357 Mon Sep 17 00:00:00 2001 From: Daniel Pettersson Date: Tue, 19 Jan 2021 21:04:27 +0100 Subject: [PATCH 2/5] Remove clone since Points are Copy --- geo/src/algorithm/geodesic_intermediate.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geo/src/algorithm/geodesic_intermediate.rs b/geo/src/algorithm/geodesic_intermediate.rs index af1c3b3f7..ce914706a 100644 --- a/geo/src/algorithm/geodesic_intermediate.rs +++ b/geo/src/algorithm/geodesic_intermediate.rs @@ -131,7 +131,7 @@ mod tests { let p2 = Point::::new(40.0, 50.0); let max_dist = 1000000.0; // meters let include_ends = true; - let i50 = p1.clone().geodesic_intermediate(&p2, 0.5); + let i50 = p1.geodesic_intermediate(&p2, 0.5); let route = p1.geodesic_intermediate_fill(&p2, max_dist, include_ends); assert_eq!(route, vec![p1, i50, p2]); } From dbbf6032adcd0b604cabe9727265daef6f64304a Mon Sep 17 00:00:00 2001 From: Daniel Pettersson Date: Tue, 19 Jan 2021 21:04:53 +0100 Subject: [PATCH 3/5] Compare Points directly --- geo/src/algorithm/geodesic_intermediate.rs | 24 ++++++++-------------- 1 file changed, 8 insertions(+), 16 deletions(-) diff --git a/geo/src/algorithm/geodesic_intermediate.rs b/geo/src/algorithm/geodesic_intermediate.rs index ce914706a..38a6120fa 100644 --- a/geo/src/algorithm/geodesic_intermediate.rs +++ b/geo/src/algorithm/geodesic_intermediate.rs @@ -22,12 +22,9 @@ pub trait GeodesicIntermediate { /// let i20_should = Point::new(29.842907, 29.951445); /// let i50_should = Point::new(65.879360, 37.722253); /// let i80_should = Point::new(103.556796, 33.506196); - /// assert_relative_eq!(i20.x(), i20_should.x(), epsilon = 1.0e-6); - /// assert_relative_eq!(i20.y(), i20_should.y(), epsilon = 1.0e-6); - /// assert_relative_eq!(i50.x(), i50_should.x(), epsilon = 1.0e-6); - /// assert_relative_eq!(i50.y(), i50_should.y(), epsilon = 1.0e-6); - /// assert_relative_eq!(i80.x(), i80_should.x(), epsilon = 1.0e-6); - /// assert_relative_eq!(i80.y(), i80_should.y(), epsilon = 1.0e-6); + /// assert_relative_eq!(i20, i20_should, epsilon = 1.0e-6); + /// assert_relative_eq!(i50, i50_should, epsilon = 1.0e-6); + /// assert_relative_eq!(i80, i80_should, epsilon = 1.0e-6); /// ``` fn geodesic_intermediate(&self, other: &Point, f: T) -> Point; @@ -101,10 +98,8 @@ mod tests { let p2 = Point::::new(15.0, 25.0); let i0 = p1.geodesic_intermediate(&p2, 0.0); let i100 = p1.geodesic_intermediate(&p2, 1.0); - assert_relative_eq!(i0.x(), p1.x(), epsilon = 1.0e-6); - assert_relative_eq!(i0.y(), p1.y(), epsilon = 1.0e-6); - assert_relative_eq!(i100.x(), p2.x(), epsilon = 1.0e-6); - assert_relative_eq!(i100.y(), p2.y(), epsilon = 1.0e-6); + assert_relative_eq!(i0, p1, epsilon = 1.0e-6); + assert_relative_eq!(i100, p2, epsilon = 1.0e-6); } #[test] @@ -117,12 +112,9 @@ mod tests { let i20_should = Point::new(29.842907, 29.951445); let i50_should = Point::new(65.879360, 37.722253); let i80_should = Point::new(103.556796, 33.506196); - assert_relative_eq!(i20.x(), i20_should.x(), epsilon = 1.0e-6); - assert_relative_eq!(i20.y(), i20_should.y(), epsilon = 1.0e-6); - assert_relative_eq!(i50.x(), i50_should.x(), epsilon = 1.0e-6); - assert_relative_eq!(i50.y(), i50_should.y(), epsilon = 1.0e-6); - assert_relative_eq!(i80.x(), i80_should.x(), epsilon = 1.0e-6); - assert_relative_eq!(i80.y(), i80_should.y(), epsilon = 1.0e-6); + assert_relative_eq!(i20, i20_should, epsilon = 1.0e-6); + assert_relative_eq!(i50, i50_should, epsilon = 1.0e-6); + assert_relative_eq!(i80, i80_should, epsilon = 1.0e-6); } #[test] From 43afd8d4eb3e2ff2c33c510fd4a89f9f65309cb4 Mon Sep 17 00:00:00 2001 From: Daniel Pettersson Date: Tue, 19 Jan 2021 22:06:20 +0100 Subject: [PATCH 4/5] Update changelog --- geo/CHANGES.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geo/CHANGES.md b/geo/CHANGES.md index fe153b2a5..04d7691fc 100644 --- a/geo/CHANGES.md +++ b/geo/CHANGES.md @@ -2,7 +2,7 @@ ## Unreleased -* Add Changes Here +* Add `GeodesicIntermediate` algorithm ## 0.17.0 From a099b3b9a9d2decc028e65bc4a6cfd0401caddd2 Mon Sep 17 00:00:00 2001 From: Daniel Pettersson Date: Sat, 23 Jan 2021 20:00:34 +0100 Subject: [PATCH 5/5] Use add assign operator --- geo/src/algorithm/geodesic_intermediate.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geo/src/algorithm/geodesic_intermediate.rs b/geo/src/algorithm/geodesic_intermediate.rs index 38a6120fa..57a286bea 100644 --- a/geo/src/algorithm/geodesic_intermediate.rs +++ b/geo/src/algorithm/geodesic_intermediate.rs @@ -76,7 +76,7 @@ impl GeodesicIntermediate for Point { g.direct(self.lat(), self.lng(), azi1, total_distance * current_step); let point = Point::new(lon2, lat2); points.push(point); - current_step = current_step + interval; + current_step += interval; } if include_ends {