Add CrossTrackDistance calculation, should be independently useful, b…

…ut is also a starting point for the implementation of #834. (#961)

* Add CrossTrackDistance calculation.

* Remove confusing duplicate test.

---------

Co-authored-by: Michael Kirk <michael.code@endoftheworl.de>

geo/CHANGES.md

@@ -2,6 +2,9 @@
 
 ## UNRELEASED
 
+- Added `CrossTrackDistance` trait to calculate the distance from a point
+  to the nearest point on a line
+  - <https://github.com/georust/geo/pull/961>
 - BREAKING: Remove deprecated methods
   - <https://github.com/georust/geo/pull/1012>
   - Instead of `map_coords_inplace` use `map_coords_in_place`
@@ -89,7 +92,7 @@
   - <https://github.com/georust/geo/pull/880>
 - Removed deprecated `ToGeo` trait. Use `std::convert::TryFrom<$geometry>`
   instead.
-  - <https://github.com/georust/geo/pull/892>
+  * <https://github.com/georust/geo/pull/892>
 
 ## 0.22.1
 

geo/src/algorithm/cross_track_distance.rs

@@ -0,0 +1,116 @@
+use num_traits::FromPrimitive;
+use geo_types::{CoordFloat, Point};
+use crate::{Bearing, HaversineDistance, MEAN_EARTH_RADIUS};
+
+/// Determine the cross track distance (also known as the cross track error) which is the shortest
+/// distance between a point and a continuous line.
+pub trait CrossTrackDistance<T, Rhs = Self> {
+    /// Determine the cross track distance between this point and a line
+    /// which passes through line_point_a and line_point_b
+    ///
+    /// # Units
+    ///
+    /// - return value: meters
+    ///
+    /// # Example
+    ///
+    /// ```rust
+    /// use geo::prelude::*;
+    /// use geo::point;
+    ///
+    /// // New York City
+    /// let p1 = point!(x: -74.006f64, y: 40.7128f64);
+    ///
+    /// // Miami
+    /// let line_point_a = point!(x: -80.1918f64, y: 25.7617f64);
+    ///
+    /// // Washington
+    /// let line_point_b = point!(x: -120.7401, y: 47.7511f64);
+    ///
+    /// let distance = p1.cross_track_distance(&line_point_a, &line_point_b);
+    ///
+    /// assert_eq!(
+    ///     1_547_104., // meters
+    ///     distance.round()
+    /// );
+    /// ```
+    fn cross_track_distance(&self, line_point_a: &Rhs, line_point_b: &Rhs) -> T;
+}
+
+impl<T> CrossTrackDistance<T, Point<T>> for Point<T>
+    where
+        T: CoordFloat + FromPrimitive,
+{
+    fn cross_track_distance(&self, line_point_a: &Point<T>, line_point_b: &Point<T>) -> T {
+        let mean_earth_radius = T::from(MEAN_EARTH_RADIUS).unwrap();
+        let l_delta_13: T = line_point_a.haversine_distance(self) / mean_earth_radius;
+        let theta_13: T = line_point_a.bearing(*self).to_radians();
+        let theta_12: T = line_point_a.bearing(*line_point_b).to_radians();
+        let l_delta_xt: T = (l_delta_13.sin() * (theta_12 - theta_13).sin()).asin();
+        mean_earth_radius * l_delta_xt.abs()
+    }
+}
+
+#[cfg(test)]
+mod test {
+    use crate::CrossTrackDistance;
+    use crate::HaversineDistance;
+    use crate::Point;
+
+    #[test]
+    fn distance1_test() {
+        let p = Point::new(-0.7972, 53.2611);
+        let line_point_a = Point::new(-1.7297, 53.3206);
+        let line_point_b = Point::new(0.1334, 53.1887);
+        assert_relative_eq!(
+            p.cross_track_distance(&line_point_a, &line_point_b),
+            307.549995,
+            epsilon = 1.0e-6
+        );
+    }
+
+    #[test]
+    fn cross_track_distance_to_line_passing_through_point() {
+        let p = Point::new(0., 0.);
+        let line_point_a = Point::new(1., 0.);
+        let line_point_b = Point::new(2., 0.);
+
+        assert_relative_eq!(
+            p.cross_track_distance(&line_point_a, &line_point_b),
+            0.,
+            epsilon = 1.0e-6
+        );
+    }
+
+    #[test]
+    fn cross_track_distance_to_line_orthogonal_to_point() {
+        let p = Point::new(0., 0.);
+        let line_point_a = Point::new(1., -1.);
+        let line_point_b = Point::new(1., 1.);
+
+        assert_relative_eq!(
+            p.cross_track_distance(&line_point_a, &line_point_b),
+            p.haversine_distance(&Point::new(1., 0.)),
+            epsilon = 1.0e-6
+        );
+
+        assert_relative_eq!(
+            p.cross_track_distance(&line_point_b, &line_point_a),
+            p.haversine_distance(&Point::new(1., 0.)),
+            epsilon = 1.0e-6
+        );
+    }
+
+    #[test]
+    fn new_york_to_line_between_miami_and_washington() {
+        let p1 = Point::new(-74.006f64, 40.7128f64);
+        let line_point_a = Point::new(-80.1918f64, 25.7617f64);
+        let line_point_b = Point::new(-120.7401f64, 47.7511f64);
+
+        assert_relative_eq!(
+            p1.cross_track_distance(&line_point_a, &line_point_b),
+            1_547_104.,
+            epsilon = 1.0
+        );
+    }
+}

geo/src/algorithm/mod.rs

@@ -60,6 +60,10 @@ pub use convert::{Convert, TryConvert};
 pub mod convex_hull;
 pub use convex_hull::ConvexHull;
 
+/// Cross track distance
+pub mod cross_track_distance;
+pub use cross_track_distance::CrossTrackDistance;
+
 /// Determine whether a `Coord` lies inside, outside, or on the boundary of a geometry.
 pub mod coordinate_position;
 pub use coordinate_position::CoordinatePosition;
@@ -243,4 +247,6 @@ pub mod sweep;
 
 /// Detect outliers in a group of points using [LOF](https://en.wikipedia.org/wiki/Local_outlier_factor)
 pub mod outlier_detection;
+
+
 pub use outlier_detection::OutlierDetection;