Skip to content

Commit

Permalink
CrossTrackDistance for other metric spaces
Browse files Browse the repository at this point in the history
  • Loading branch information
michaelkirk committed Oct 31, 2024
1 parent 348f622 commit 36086ab
Show file tree
Hide file tree
Showing 3 changed files with 43 additions and 26 deletions.
1 change: 1 addition & 0 deletions geo/CHANGES.md
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@

## Unreleased

* BREAKING: Make `CrossTrackDistance` generic both to clarify when it is using Haversine, but also implement it for Euclidean.

## 0.29.0 - 2024.10.30

Expand Down
55 changes: 33 additions & 22 deletions geo/src/algorithm/cross_track_distance.rs
Original file line number Diff line number Diff line change
@@ -1,10 +1,13 @@
use crate::{Bearing, Distance, Haversine, MEAN_EARTH_RADIUS};
use geo_types::{CoordFloat, Point};
use crate::Distance;
use geo_types::{CoordFloat, Line, Point};
use num_traits::FromPrimitive;

/// 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> {
pub trait CrossTrackDistance<T, Rhs = Self>
where
T: CoordFloat + FromPrimitive,
{
/// Determine the cross track distance between this point and a line
/// which passes through line_point_a and line_point_b
///
Expand Down Expand Up @@ -34,36 +37,39 @@ pub trait CrossTrackDistance<T, Rhs = Self> {
/// distance.round()
/// );
/// ```
fn cross_track_distance(&self, line_point_a: &Rhs, line_point_b: &Rhs) -> T;
fn cross_track_distance<MetricSpace>(&self, line_point_a: &Rhs, line_point_b: &Rhs) -> T
where
MetricSpace: for<'a> Distance<T, &'a Line<T>, &'a Point<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 = Haversine::distance(*line_point_a, *self) / mean_earth_radius;
let theta_13: T = Haversine::bearing(*line_point_a, *self).to_radians();
let theta_12: T = Haversine::bearing(*line_point_a, *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()
fn cross_track_distance<MetricSpace>(
&self,
line_point_a: &Point<T>,
line_point_b: &Point<T>,
) -> T
where
MetricSpace: for<'a> Distance<T, &'a Line<T>, &'a Point<T>>,
{
MetricSpace::distance(&Line::new(*line_point_a, *line_point_b), self)
}
}

#[cfg(test)]
mod test {
use crate::CrossTrackDistance;
use crate::Point;
use crate::{Distance, Haversine};
use crate::{CrossTrackDistance, Distance, Euclidean, Haversine};

#[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),
p.cross_track_distance::<Haversine>(&line_point_a, &line_point_b),
307.549995,
epsilon = 1.0e-6
);
Expand All @@ -76,7 +82,7 @@ mod test {
let line_point_b = Point::new(2., 0.);

assert_relative_eq!(
p.cross_track_distance(&line_point_a, &line_point_b),
p.cross_track_distance::<Haversine>(&line_point_a, &line_point_b),
0.,
epsilon = 1.0e-6
);
Expand All @@ -89,13 +95,13 @@ mod test {
let line_point_b = Point::new(1., 1.);

assert_relative_eq!(
p.cross_track_distance(&line_point_a, &line_point_b),
p.cross_track_distance::<Haversine>(&line_point_a, &line_point_b),
Haversine::distance(p, Point::new(1., 0.)),
epsilon = 1.0e-6
);

assert_relative_eq!(
p.cross_track_distance(&line_point_b, &line_point_a),
p.cross_track_distance::<Haversine>(&line_point_b, &line_point_a),
Haversine::distance(p, Point::new(1., 0.)),
epsilon = 1.0e-6
);
Expand All @@ -107,10 +113,15 @@ mod test {
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
);
let haversine_distance = p1.cross_track_distance::<Haversine>(&line_point_a, &line_point_b);
assert_relative_eq!(haversine_distance, 1_547_104., epsilon = 1.0);

// Same as above, but projected to EPSG:5070 to test Euclidean
let p1 = Point::new(1826303.9422258963, 2179112.0732980534);
let line_point_a = Point::new(1594077.349564382, 434470.0414052719);
let line_point_b = Point::new(-1847681.2454118263, 2992574.0850278544);

let euclidean_distance = p1.cross_track_distance::<Euclidean>(&line_point_a, &line_point_b);
assert_relative_eq!(euclidean_distance, 1_538_764., epsilon = 1.0);
}
}
13 changes: 9 additions & 4 deletions geo/src/algorithm/line_measures/metric_spaces/haversine.rs
Original file line number Diff line number Diff line change
Expand Up @@ -152,12 +152,17 @@ impl<F: CoordFloat + FromPrimitive> Distance<F, Point<F>, Point<F>> for Haversin
}
}

impl<F: CoordFloat + FromPrimitive> Distance<F, Line<F>, Point<F>> for Haversine {
fn distance(origin: Line<F>, destination: Point<F>) -> F {
todo!()
impl<F: CoordFloat + FromPrimitive> Distance<F, &Line<F>, &Point<F>> for Haversine {
fn distance(line: &Line<F>, point: &Point<F>) -> F {
let mean_earth_radius = F::from(MEAN_EARTH_RADIUS).unwrap();
let l_delta_13 = Haversine::distance(line.start_point(), *point) / mean_earth_radius;
let theta_13 = Haversine::bearing(line.start_point(), *point).to_radians();
let theta_12 = Haversine::bearing(line.start_point(), line.end_point()).to_radians();
let l_delta_xt = (l_delta_13.sin() * (theta_12 - theta_13).sin()).asin();
mean_earth_radius * l_delta_xt.abs()
}
}
symmetric_distance_impl!(Point<F>, Line<F>, for: Haversine, where: CoordFloat + FromPrimitive);
symmetric_distance_impl!(&Point<F>, &Line<F>, for: Haversine, where: CoordFloat + FromPrimitive);

/// Interpolate Point(s) along a [great circle].
///
Expand Down

0 comments on commit 36086ab

Please sign in to comment.