Rhumb no longer goes past pole

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.
   * <https://github.com/georust/geo/pull/1134>
+* BREAKING: `RhumbBearing` now returns an `Option<Point<T>>` rather than a raw
+  `Point<T>` 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

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);

geo/src/algorithm/rhumb/destination.rs

@@ -12,6 +12,9 @@ pub trait RhumbDestination<T: CoordFloat> {
     /// 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<T: CoordFloat> {
     /// 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<T>;
+    fn rhumb_destination(&self, bearing: T, distance: T) -> Option<Point<T>>;
 }
 
 impl<T> RhumbDestination<T> for Point<T>
 where
     T: CoordFloat + FromPrimitive,
 {
-    fn rhumb_destination(&self, bearing: T, distance: T) -> Point<T> {
+    fn rhumb_destination(&self, bearing: T, distance: T) -> Option<Point<T>> {
         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());
     }
 }

geo/src/algorithm/rhumb/intermediate.rs

@@ -8,6 +8,8 @@ use super::RhumbCalculations;
 pub trait RhumbIntermediate<T: CoordFloat> {
     /// 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]);
+    }
 }

geo/src/algorithm/rhumb/mod.rs

@@ -35,12 +35,29 @@ struct RhumbCalculations<T: CoordFloat + FromPrimitive> {
 
 impl<T: CoordFloat + FromPrimitive> RhumbCalculations<T> {
     fn new(from: &Point<T>, to: &Point<T>) -> 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 {
@@ -53,9 +70,9 @@ impl<T: CoordFloat + FromPrimitive> RhumbCalculations<T> {
         let delta_psi = ((phi2 / two + pi / four).tan() / (phi1 / two + pi / four).tan()).ln();
         let delta_phi = phi2 - phi1;
 
-        RhumbCalculations {
-            from: *from,
-            to: *to,
+        Self {
+            from: from,
+            to: to,
             phi1,
             delta_lambda,
             delta_phi,
@@ -82,7 +99,7 @@ impl<T: CoordFloat + FromPrimitive> RhumbCalculations<T> {
         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<Point<T>> {
@@ -112,7 +129,7 @@ impl<T: CoordFloat + FromPrimitive> RhumbCalculations<T> {
         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 +146,18 @@ fn calculate_destination<T: CoordFloat + FromPrimitive>(
     lambda1: T,
     phi1: T,
     theta: T,
-) -> Point<T> {
+) -> Option<Point<T>> {
     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 +171,8 @@ fn calculate_destination<T: CoordFloat + FromPrimitive>(
     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(),
-    }
+    })
 }

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;