This commit is deliberately borked to show why clamping is impractical

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 {
@@ -54,8 +71,8 @@ impl<T: CoordFloat + FromPrimitive> RhumbCalculations<T> {
         let delta_phi = phi2 - phi1;
 
         Self {
-            from: *from,
-            to: *to,
+            from: from,
+            to: to,
             phi1,
             delta_lambda,
             delta_phi,