From c9fa18fefa7aef83cccbd5d677eb40846a0e4236 Mon Sep 17 00:00:00 2001 From: David Morgan-Brown Date: Tue, 16 May 2023 21:12:33 +0100 Subject: [PATCH] Add CrossTrackDistance calculation, should be independently useful, but is also a starting point for the implementation of #834. (#961) * Add CrossTrackDistance calculation. * Remove confusing duplicate test. --------- Co-authored-by: Michael Kirk --- geo/CHANGES.md | 5 +- geo/src/algorithm/cross_track_distance.rs | 116 ++++++++++++++++++++++ geo/src/algorithm/mod.rs | 6 ++ 3 files changed, 126 insertions(+), 1 deletion(-) create mode 100644 geo/src/algorithm/cross_track_distance.rs diff --git a/geo/CHANGES.md b/geo/CHANGES.md index fb046d83dc..e3b1ad0336 100644 --- a/geo/CHANGES.md +++ b/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 + - - BREAKING: Remove deprecated methods - - Instead of `map_coords_inplace` use `map_coords_in_place` @@ -89,7 +92,7 @@ - - Removed deprecated `ToGeo` trait. Use `std::convert::TryFrom<$geometry>` instead. - - + * ## 0.22.1 diff --git a/geo/src/algorithm/cross_track_distance.rs b/geo/src/algorithm/cross_track_distance.rs new file mode 100644 index 0000000000..cbd8829092 --- /dev/null +++ b/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 { + /// 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 CrossTrackDistance> for Point + where + T: CoordFloat + FromPrimitive, +{ + fn cross_track_distance(&self, line_point_a: &Point, line_point_b: &Point) -> 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 + ); + } +} \ No newline at end of file diff --git a/geo/src/algorithm/mod.rs b/geo/src/algorithm/mod.rs index 49f97048e0..d7e57b590f 100644 --- a/geo/src/algorithm/mod.rs +++ b/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;