Merge pull request #1090 from apendleton/rhumb-ops

Add rhumb-line operations

Author: urschrei

geo/CHANGES.md

@@ -18,6 +18,8 @@
   * <https://github.com/georust/geo/pull/908>
 * Add `ToDegrees` and `ToRadians` traits.
   * <https://github.com/georust/geo/pull/1070>
+* Add rhumb-line operations analogous to several current haversine operations: `RhumbBearing`, `RhumbDestination`, `RhumbDistance`, `RhumbIntermediate`, `RhumbLength`.
+  * <https://github.com/georust/geo/pull/1090>
 
 ## 0.26.0
 

geo/src/algorithm/mod.rs

@@ -283,3 +283,7 @@ pub use outlier_detection::OutlierDetection;
 /// Monotonic polygon subdivision
 pub mod monotone;
 pub use monotone::{monotone_subdivision, MonoPoly, MonotonicPolygons};
+
+/// Rhumb-line-related algorithms and utils
+pub mod rhumb;
+pub use rhumb::{RhumbBearing, RhumbDestination, RhumbDistance, RhumbIntermediate, RhumbLength};

geo/src/algorithm/rhumb/bearing.rs

@@ -0,0 +1,91 @@
+use num_traits::FromPrimitive;
+
+use crate::{CoordFloat, Point};
+
+use super::RhumbCalculations;
+
+/// Returns the bearing to another Point in degrees.
+///
+/// Bullock, R.: Great Circle Distances and Bearings Between Two Locations, 2007.
+/// (<https://dtcenter.org/met/users/docs/write_ups/gc_simple.pdf>)
+
+pub trait RhumbBearing<T: CoordFloat + FromPrimitive> {
+    /// Returns the bearing to another Point in degrees along a [rhumb line], where North is 0° and East is 90°.
+    ///
+    /// # Examples
+    ///
+    /// ```
+    /// # use approx::assert_relative_eq;
+    /// use geo::RhumbBearing;
+    /// use geo::Point;
+    ///
+    /// let p_1 = Point::new(9.177789688110352, 48.776781529534965);
+    /// let p_2 = Point::new(9.274348757829898, 48.84037308229984);
+    /// let bearing = p_1.rhumb_bearing(p_2);
+    /// assert_relative_eq!(bearing, 45., epsilon = 1.0e-6);
+    /// ```
+    /// [rhumb line]: https://en.wikipedia.org/wiki/Rhumb_line
+
+    fn rhumb_bearing(&self, point: Point<T>) -> T;
+}
+
+impl<T> RhumbBearing<T> for Point<T>
+where
+    T: CoordFloat + FromPrimitive,
+{
+    fn rhumb_bearing(&self, point: Point<T>) -> T {
+        let three_sixty = T::from(360.0f64).unwrap();
+
+        let calculations = RhumbCalculations::new(self, &point);
+        (calculations.theta().to_degrees() + three_sixty) % three_sixty
+    }
+}
+
+#[cfg(test)]
+mod test {
+    use crate::point;
+    use crate::RhumbBearing;
+    use crate::RhumbDestination;
+
+    #[test]
+    fn north_bearing() {
+        let p_1 = point!(x: 9., y: 47.);
+        let p_2 = point!(x: 9., y: 48.);
+        let bearing = p_1.rhumb_bearing(p_2);
+        assert_relative_eq!(bearing, 0.);
+    }
+
+    #[test]
+    fn equatorial_east_bearing() {
+        let p_1 = point!(x: 9., y: 0.);
+        let p_2 = point!(x: 10., y: 0.);
+        let bearing = p_1.rhumb_bearing(p_2);
+        assert_relative_eq!(bearing, 90.);
+    }
+
+    #[test]
+    fn east_bearing() {
+        let p_1 = point!(x: 9., y: 10.);
+        let p_2 = point!(x: 18.131938299366652, y: 10.);
+
+        let bearing = p_1.rhumb_bearing(p_2);
+        assert_relative_eq!(bearing, 90.);
+    }
+
+    #[test]
+    fn northeast_bearing() {
+        let p_1 = point!(x: 9.177789688110352f64, y: 48.776781529534965);
+        let p_2 = point!(x: 9.274348757829898, y: 48.84037308229984);
+        let bearing = p_1.rhumb_bearing(p_2);
+        assert_relative_eq!(bearing, 45., epsilon = 1.0e-6);
+    }
+
+    #[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 b_1 = p_1.rhumb_bearing(p_2);
+        assert_relative_eq!(b_1, 45., epsilon = 1.0e-6);
+    }
+}

geo/src/algorithm/rhumb/destination.rs

@@ -0,0 +1,80 @@
+use crate::{CoordFloat, Point, MEAN_EARTH_RADIUS};
+use num_traits::FromPrimitive;
+
+use super::calculate_destination;
+
+/// Returns the destination Point having travelled the given distance along a [rhumb line]
+/// from the origin geometry with the given bearing
+///
+/// *Note*: this implementation uses a mean earth radius of 6371.088 km, based on the [recommendation of
+/// the IUGG](ftp://athena.fsv.cvut.cz/ZFG/grs80-Moritz.pdf)
+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
+    ///
+    /// # Units
+    ///
+    /// - `bearing`: degrees, zero degrees is north
+    /// - `distance`: meters
+    ///
+    /// # Examples
+    ///
+    /// ```
+    /// use geo::RhumbDestination;
+    /// use geo::Point;
+    ///
+    /// let p_1 = Point::new(9.177789688110352, 48.776781529534965);
+    /// let p_2 = p_1.rhumb_destination(45., 10000.);
+    /// 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>;
+}
+
+impl<T> RhumbDestination<T> for Point<T>
+where
+    T: CoordFloat + FromPrimitive,
+{
+    fn rhumb_destination(&self, bearing: T, distance: T) -> 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();
+        let theta = bearing.to_radians();
+
+        calculate_destination(delta, lambda1, phi1, theta)
+    }
+}
+
+#[cfg(test)]
+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.);
+        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)
+    }
+
+    #[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);
+        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.);
+        assert_relative_eq!(p_1.x(), p_2.x(), epsilon = 1.0e-6);
+        assert!(p_2.y() > p_1.y())
+    }
+}

geo/src/algorithm/rhumb/distance.rs

@@ -0,0 +1,99 @@
+use crate::{CoordFloat, Point, MEAN_EARTH_RADIUS};
+use num_traits::FromPrimitive;
+
+use super::RhumbCalculations;
+
+/// Determine the distance between two geometries along a [rhumb line].
+///
+/// [rhumb line]: https://en.wikipedia.org/wiki/Rhumb_line
+///
+/// *Note*: this implementation uses a mean earth radius of 6371.088 km, based on the [recommendation of
+/// the IUGG](ftp://athena.fsv.cvut.cz/ZFG/grs80-Moritz.pdf)
+pub trait RhumbDistance<T, Rhs = Self> {
+    /// Determine the distance between along the [rhumb line] between two geometries.
+    ///
+    /// # Units
+    ///
+    /// - return value: meters
+    ///
+    /// # Examples
+    ///
+    /// ```rust
+    /// use geo::prelude::*;
+    /// use geo::point;
+    ///
+    /// // New York City
+    /// let p1 = point!(x: -74.006f64, y: 40.7128f64);
+    ///
+    /// // London
+    /// let p2 = point!(x: -0.1278f64, y: 51.5074f64);
+    ///
+    /// let distance = p1.rhumb_distance(&p2);
+    ///
+    /// assert_eq!(
+    ///     5_794_129., // meters
+    ///     distance.round()
+    /// );
+    /// ```
+    ///
+    /// [rhumb line]: https://en.wikipedia.org/wiki/Rhumb_line
+    fn rhumb_distance(&self, rhs: &Rhs) -> T;
+}
+
+impl<T> RhumbDistance<T, Point<T>> for Point<T>
+where
+    T: CoordFloat + FromPrimitive,
+{
+    fn rhumb_distance(&self, rhs: &Point<T>) -> T {
+        let calculations = RhumbCalculations::new(self, rhs);
+        calculations.delta() * T::from(MEAN_EARTH_RADIUS).unwrap()
+    }
+}
+
+#[cfg(test)]
+mod test {
+    use crate::Point;
+    use crate::RhumbDistance;
+
+    #[test]
+    fn distance1_test() {
+        let a = Point::new(0., 0.);
+        let b = Point::new(1., 0.);
+        assert_relative_eq!(
+            a.rhumb_distance(&b),
+            111195.0802335329_f64,
+            epsilon = 1.0e-6
+        );
+    }
+
+    #[test]
+    fn distance2_test() {
+        let a = Point::new(-72.1235, 42.3521);
+        let b = Point::new(72.1260, 70.612);
+        assert_relative_eq!(
+            a.rhumb_distance(&b),
+            8903668.508603323_f64,
+            epsilon = 1.0e-6
+        );
+    }
+
+    #[test]
+    fn distance3_test() {
+        // this input comes from issue #100
+        let a = Point::new(-77.036585, 38.897448);
+        let b = Point::new(-77.009080, 38.889825);
+        assert_relative_eq!(
+            a.rhumb_distance(&b),
+            2526.7031699343006_f64,
+            epsilon = 1.0e-6
+        );
+    }
+
+    #[test]
+    fn distance3_test_f32() {
+        // this input comes from issue #100
+        let a = Point::<f32>::new(-77.03658, 38.89745);
+        let b = Point::<f32>::new(-77.00908, 38.889825);
+        assert_relative_eq!(a.rhumb_distance(&b), 2526.7273_f32, epsilon = 1.0e-6);
+    }
+}