geo/algorithm/
geodesic_destination.rs

1use crate::Point;
2use geo_types::CoordNum;
3use geographiclib_rs::{DirectGeodesic, Geodesic};
4
5/// Returns a new Point using the distance to the existing Point and a bearing for the direction on a geodesic.
6///
7/// This uses the geodesic methods given by [Karney (2013)].
8///
9/// [Karney (2013)]:  https://arxiv.org/pdf/1109.4448.pdf
10pub trait GeodesicDestination<T: CoordNum> {
11    /// Returns a new Point using distance to the existing Point and a bearing for the direction.
12    ///
13    /// # Units
14    ///
15    /// - `bearing`: degrees, zero degrees is north
16    /// - `distance`: meters
17    ///
18    /// # Examples
19    ///
20    /// ```rust
21    /// use geo::GeodesicDestination;
22    /// use geo::Point;
23    ///
24    /// // Determine the point 10000 km NE of JFK.
25    /// let jfk = Point::new(-73.78, 40.64);
26    /// let northeast_bearing = 45.0;
27    /// let distance = 10e6;
28    ///
29    /// let p_1 = jfk.geodesic_destination(northeast_bearing, distance);
30    /// use approx::assert_relative_eq;
31    /// assert_relative_eq!(p_1.x(), 49.052487092959836);
32    /// assert_relative_eq!(p_1.y(), 32.621100463725796);
33    /// ```
34    fn geodesic_destination(&self, bearing: T, distance: T) -> Point<T>;
35}
36
37impl GeodesicDestination<f64> for Point<f64> {
38    fn geodesic_destination(&self, bearing: f64, distance: f64) -> Point<f64> {
39        let (lat, lon) = Geodesic::wgs84().direct(self.y(), self.x(), bearing, distance);
40        Point::new(lon, lat)
41    }
42}
43
44#[cfg(test)]
45mod test {
46    use super::*;
47    use crate::GeodesicDistance;
48
49    #[test]
50    fn returns_a_new_point() {
51        let p_1 = Point::new(9.177789688110352, 48.776781529534965);
52        let p_2 = p_1.geodesic_destination(45., 10000.);
53
54        assert_relative_eq!(
55            p_2,
56            Point::new(9.27411867078536, 48.8403266058781),
57            epsilon = 1.0e-6
58        );
59
60        let distance = p_1.geodesic_distance(&p_2);
61        assert_relative_eq!(distance, 10000., epsilon = 1.0e-6)
62    }
63
64    #[test]
65    fn bearing_zero_is_north() {
66        let p_1 = Point::new(9.177789688110352, 48.776781529534965);
67        let p_2 = p_1.geodesic_destination(0., 1000.);
68        assert_relative_eq!(p_1.x(), p_2.x(), epsilon = 1.0e-6);
69        assert!(p_2.y() > p_1.y())
70    }
71
72    #[test]
73    fn bearing_90_is_east() {
74        let p_1 = Point::new(9.177789688110352, 48.776781529534965);
75        let p_2 = p_1.geodesic_destination(90., 1000.);
76        assert_relative_eq!(p_1.y(), p_2.y(), epsilon = 1.0e-6);
77        assert!(p_2.x() > p_1.x())
78    }
79}