From b4eeaa098fc3b3a1711fa6388d3fb52755b64880 Mon Sep 17 00:00:00 2001 From: Joseph Fox-Rabinovitz Date: Wed, 21 Feb 2024 20:34:53 -0600 Subject: [PATCH 1/3] Rhumb no longer goes past pole --- geo/src/algorithm/rhumb/bearing.rs | 2 +- geo/src/algorithm/rhumb/destination.rs | 33 ++++++++++++++++---------- geo/src/algorithm/rhumb/mod.rs | 20 +++++++--------- geo/src/lib.rs | 2 +- 4 files changed, 31 insertions(+), 26 deletions(-) diff --git a/geo/src/algorithm/rhumb/bearing.rs b/geo/src/algorithm/rhumb/bearing.rs index c1c6e3c85..e6d18b7a5 100644 --- a/geo/src/algorithm/rhumb/bearing.rs +++ b/geo/src/algorithm/rhumb/bearing.rs @@ -83,7 +83,7 @@ mod test { #[test] fn consistent_with_destination() { let p_1 = point!(x: 9.177789688110352f64, y: 48.776781529534965); - let p_2 = p_1.rhumb_destination(45., 10000.); + let p_2 = p_1.rhumb_destination(45., 10000.).unwrap(); let b_1 = p_1.rhumb_bearing(p_2); assert_relative_eq!(b_1, 45., epsilon = 1.0e-6); diff --git a/geo/src/algorithm/rhumb/destination.rs b/geo/src/algorithm/rhumb/destination.rs index a9b9be85e..774c62b62 100644 --- a/geo/src/algorithm/rhumb/destination.rs +++ b/geo/src/algorithm/rhumb/destination.rs @@ -12,6 +12,9 @@ pub trait RhumbDestination { /// Returns the destination Point having travelled the given distance along a [rhumb line] /// from the origin Point with the given bearing /// + /// A rhumb line has a finite length, so if the path distance exceeds the location of the pole, + /// the result is None. + /// /// # Units /// /// - `bearing`: degrees, zero degrees is north @@ -24,18 +27,18 @@ pub trait RhumbDestination { /// use geo::Point; /// /// let p_1 = Point::new(9.177789688110352, 48.776781529534965); - /// let p_2 = p_1.rhumb_destination(45., 10000.); + /// let p_2 = p_1.rhumb_destination(45., 10000.).unwrap(); /// assert_eq!(p_2, Point::new(9.274348757829898, 48.84037308229984)) /// ``` /// [rhumb line]: https://en.wikipedia.org/wiki/Rhumb_line - fn rhumb_destination(&self, bearing: T, distance: T) -> Point; + fn rhumb_destination(&self, bearing: T, distance: T) -> Option>; } impl RhumbDestination for Point where T: CoordFloat + FromPrimitive, { - fn rhumb_destination(&self, bearing: T, distance: T) -> Point { + fn rhumb_destination(&self, bearing: T, distance: T) -> Option> { let delta = distance / T::from(MEAN_EARTH_RADIUS).unwrap(); // angular distance in radians let lambda1 = self.x().to_radians(); let phi1 = self.y().to_radians(); @@ -49,32 +52,38 @@ where mod test { use super::*; use crate::RhumbDistance; - use num_traits::pow; #[test] fn returns_a_new_point() { let p_1 = Point::new(9.177789688110352, 48.776781529534965); - let p_2 = p_1.rhumb_destination(45., 10000.); + let p_2 = p_1.rhumb_destination(45., 10000.).unwrap(); assert_eq!(p_2, Point::new(9.274348757829898, 48.84037308229984)); let distance = p_1.rhumb_distance(&p_2); - assert_relative_eq!(distance, 10000., epsilon = 1.0e-6) + assert_relative_eq!(distance, 10000., epsilon = 1.0e-6); } #[test] fn direct_and_indirect_destinations_are_close() { let p_1 = Point::new(9.177789688110352, 48.776781529534965); - let p_2 = p_1.rhumb_destination(45., 10000.); - let square_edge = { pow(10000., 2) / 2f64 }.sqrt(); - let p_3 = p_1.rhumb_destination(0., square_edge); - let p_4 = p_3.rhumb_destination(90., square_edge); + let p_2 = p_1.rhumb_destination(45., 10000.).unwrap(); + let square_edge = 10000.0 * 0.5f64.sqrt(); + let p_3 = p_1.rhumb_destination(0., square_edge).unwrap(); + let p_4 = p_3.rhumb_destination(90., square_edge).unwrap(); assert_relative_eq!(p_4, p_2, epsilon = 1.0e-3); } #[test] fn bearing_zero_is_north() { let p_1 = Point::new(9.177789688110352, 48.776781529534965); - let p_2 = p_1.rhumb_destination(0., 1000.); + let p_2 = p_1.rhumb_destination(0., 1000.).unwrap(); assert_relative_eq!(p_1.x(), p_2.x(), epsilon = 1.0e-6); - assert!(p_2.y() > p_1.y()) + assert!(p_2.y() > p_1.y()); + } + + #[test] + fn past_the_pole() { + let p = Point::new(9.177789688110352, 48.776781529534965); + let result = p.rhumb_destination(0., 1e8); + assert!(result.is_none()); } } diff --git a/geo/src/algorithm/rhumb/mod.rs b/geo/src/algorithm/rhumb/mod.rs index 7ffe9f8b7..7ea42c575 100644 --- a/geo/src/algorithm/rhumb/mod.rs +++ b/geo/src/algorithm/rhumb/mod.rs @@ -53,7 +53,7 @@ impl RhumbCalculations { let delta_psi = ((phi2 / two + pi / four).tan() / (phi1 / two + pi / four).tan()).ln(); let delta_phi = phi2 - phi1; - RhumbCalculations { + Self { from: *from, to: *to, phi1, @@ -82,7 +82,7 @@ impl RhumbCalculations { let delta = fraction * self.delta(); let theta = self.theta(); let lambda1 = self.from.x().to_radians(); - calculate_destination(delta, lambda1, self.phi1, theta) + calculate_destination(delta, lambda1, self.phi1, theta).unwrap() } fn intermediate_fill(&self, max_delta: T, include_ends: bool) -> Vec> { @@ -112,7 +112,7 @@ impl RhumbCalculations { while current_step < T::one() { let delta = total_delta * current_step; let point = calculate_destination(delta, lambda1, self.phi1, theta); - points.push(point); + points.push(point.unwrap()); current_step = current_step + interval; } @@ -129,22 +129,18 @@ fn calculate_destination( lambda1: T, phi1: T, theta: T, -) -> Point { +) -> Option> { let pi = T::from(std::f64::consts::PI).unwrap(); let two = T::one() + T::one(); let four = two + two; let threshold = T::from(10.0e-12).unwrap(); let delta_phi = delta * theta.cos(); - let mut phi2 = phi1 + delta_phi; + let phi2 = phi1 + delta_phi; // check for some daft bugger going past the pole, normalise latitude if so if phi2.abs() > pi / two { - phi2 = if phi2 > T::zero() { - pi - phi2 - } else { - -pi - phi2 - }; + return None; } let delta_psi = ((phi2 / two + pi / four).tan() / (phi1 / two + pi / four).tan()).ln(); @@ -158,8 +154,8 @@ fn calculate_destination( let delta_lambda = (delta * theta.sin()) / q; let lambda2 = lambda1 + delta_lambda; - point! { + Some(point! { x: normalize_longitude(lambda2.to_degrees()), y: phi2.to_degrees(), - } + }) } diff --git a/geo/src/lib.rs b/geo/src/lib.rs index e5e291667..771fb36b1 100644 --- a/geo/src/lib.rs +++ b/geo/src/lib.rs @@ -242,7 +242,7 @@ extern crate log; /// https://link.springer.com/article/10.1007%2Fs001900050278 /// https://sci-hub.se/https://doi.org/10.1007/s001900050278 /// https://en.wikipedia.org/wiki/Earth_radius#Mean_radius -const MEAN_EARTH_RADIUS: f64 = 6371008.8; +const MEAN_EARTH_RADIUS: f64 = 6_371_008.8; // Radius of Earth at the equator in meters (derived from the WGS-84 ellipsoid) const EQUATORIAL_EARTH_RADIUS: f64 = 6_378_137.0; From 17b1a00a3aa7c43a5aae80d97e34c7e2f32d6b4a Mon Sep 17 00:00:00 2001 From: Joseph Fox-Rabinovitz Date: Wed, 21 Feb 2024 20:41:44 -0600 Subject: [PATCH 2/3] Updated CHANGES --- geo/CHANGES.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/geo/CHANGES.md b/geo/CHANGES.md index 2a47e9908..11a87d42a 100644 --- a/geo/CHANGES.md +++ b/geo/CHANGES.md @@ -16,6 +16,9 @@ shouldn't break for any common numeric types, but if you are using something exotic you'll need to manually implement `GeoNum` for your numeric type. * +* BREAKING: `RhumbBearing` now returns an `Option>` rather than a raw + `Point` to reflect the fact that rhumb lines have a fixed length and are + undefined at and past the poles. * POSSIBLY BREAKING: `SimplifyVwPreserve` trait implementation moved from `geo_types::CoordNum` to `geo::GeoNum` as a consequence of introducing the `GeoNum::total_cmp`. This shouldn't break anything for common numeric From af50f59e7fbb6a07a5dd7975597713f8305931c5 Mon Sep 17 00:00:00 2001 From: Joseph Fox-Rabinovitz Date: Fri, 23 Feb 2024 10:03:53 -0600 Subject: [PATCH 3/3] This commit is deliberately borked to show why clamping is impractical --- geo/src/algorithm/rhumb/intermediate.rs | 15 +++++++++++++++ geo/src/algorithm/rhumb/mod.rs | 25 +++++++++++++++++++++---- 2 files changed, 36 insertions(+), 4 deletions(-) diff --git a/geo/src/algorithm/rhumb/intermediate.rs b/geo/src/algorithm/rhumb/intermediate.rs index 54651c6ec..9bb885129 100644 --- a/geo/src/algorithm/rhumb/intermediate.rs +++ b/geo/src/algorithm/rhumb/intermediate.rs @@ -8,6 +8,8 @@ use super::RhumbCalculations; pub trait RhumbIntermediate { /// Returns a new Point along a [rhumb line] between two existing points. /// + /// Start and end points are clamped to the range [-90.0, +90.0]. + /// /// # Examples /// /// ```rust @@ -140,4 +142,17 @@ mod test { let route = p1.rhumb_intermediate_fill(&p2, max_dist, include_ends); assert_eq!(route, vec![p1, i25, i50, i75, p2]); } + + #[test] + fn clamp_on_both_ends() { + let p1 = Point::new(30.0, -100.0); + let p2 = Point::new(40.0, 100.0); + let max_dist = 12000000.0; // meters + let include_ends = true; + let i_start = Point::new(30.0, -90.0); + let i_half = Point::new(35.0, 0.0); // Vertical symmetry across equator + let i_end = Point::new(40.0, 90.0); + let route = p1.rhumb_intermediate_fill(&p2, max_dist, include_ends); + assert_eq!(route, vec![i_start, i_half, i_end]); + } } diff --git a/geo/src/algorithm/rhumb/mod.rs b/geo/src/algorithm/rhumb/mod.rs index 7ea42c575..ded9a368a 100644 --- a/geo/src/algorithm/rhumb/mod.rs +++ b/geo/src/algorithm/rhumb/mod.rs @@ -35,12 +35,29 @@ struct RhumbCalculations { impl RhumbCalculations { fn new(from: &Point, to: &Point) -> Self { + let ninety = T::from(90.0).unwrap(); let pi = T::from(std::f64::consts::PI).unwrap(); let two = T::one() + T::one(); let four = two + two; - let phi1 = from.y().to_radians(); - let phi2 = to.y().to_radians(); + let phi1 = if from.y() < -ninety { + -ninety + } else if from.y() > ninety { + ninety + } else { + from.y() + }; + let phi2 = if to.y() < -ninety { + -ninety + } else if to.y() > ninety { + ninety + } else { + to.y() + }; + let from = Point::new(from.x(), phi1); + let to = Point::new(to.x(), phi2); + let phi1 = phi1.to_radians(); + let phi2 = phi2.to_radians(); let mut delta_lambda = (to.x() - from.x()).to_radians(); // if delta_lambda is over 180° take shorter rhumb line across the anti-meridian: if delta_lambda > pi { @@ -54,8 +71,8 @@ impl RhumbCalculations { let delta_phi = phi2 - phi1; Self { - from: *from, - to: *to, + from: from, + to: to, phi1, delta_lambda, delta_phi,