geo/algorithm/
haversine_destination.rs

1use crate::{CoordFloat, Point, MEAN_EARTH_RADIUS};
2use num_traits::FromPrimitive;
3
4/// Returns a new Point using the distance to the existing Point and a bearing for the direction
5///
6/// *Note*: this implementation uses a mean earth radius of 6371.088 km, based on the [recommendation of
7/// the IUGG](ftp://athena.fsv.cvut.cz/ZFG/grs80-Moritz.pdf)
8pub trait HaversineDestination<T: CoordFloat> {
9    /// Returns a new Point using distance to the existing Point and a bearing for the direction
10    ///
11    /// # Units
12    ///
13    /// - `bearing`: degrees, zero degrees is north
14    /// - `distance`: meters
15    ///
16    /// # Examples
17    ///
18    /// ```rust
19    /// use geo::HaversineDestination;
20    /// use geo::Point;
21    ///
22    /// let p_1 = Point::new(9.177789688110352, 48.776781529534965);
23    /// let p_2 = p_1.haversine_destination(45., 10000.);
24    /// assert_eq!(p_2, Point::new(9.274409949623548, 48.84033274015048))
25    /// ```
26    fn haversine_destination(&self, bearing: T, distance: T) -> Point<T>;
27}
28
29impl<T> HaversineDestination<T> for Point<T>
30where
31    T: CoordFloat + FromPrimitive,
32{
33    fn haversine_destination(&self, bearing: T, distance: T) -> Point<T> {
34        let center_lng = self.x().to_radians();
35        let center_lat = self.y().to_radians();
36        let bearing_rad = bearing.to_radians();
37
38        let rad = distance / T::from(MEAN_EARTH_RADIUS).unwrap();
39
40        let lat =
41            { center_lat.sin() * rad.cos() + center_lat.cos() * rad.sin() * bearing_rad.cos() }
42                .asin();
43        let lng = { bearing_rad.sin() * rad.sin() * center_lat.cos() }
44            .atan2(rad.cos() - center_lat.sin() * lat.sin())
45            + center_lng;
46
47        Point::new(lng.to_degrees(), lat.to_degrees())
48    }
49}
50
51#[cfg(test)]
52mod test {
53    use super::*;
54    use crate::HaversineDistance;
55    use num_traits::pow;
56
57    #[test]
58    fn returns_a_new_point() {
59        let p_1 = Point::new(9.177789688110352, 48.776781529534965);
60        let p_2 = p_1.haversine_destination(45., 10000.);
61        assert_eq!(p_2, Point::new(9.274409949623548, 48.84033274015048));
62        let distance = p_1.haversine_distance(&p_2);
63        assert_relative_eq!(distance, 10000., epsilon = 1.0e-6)
64    }
65
66    #[test]
67    fn direct_and_indirect_destinations_are_close() {
68        let p_1 = Point::new(9.177789688110352, 48.776781529534965);
69        let p_2 = p_1.haversine_destination(45., 10000.);
70        let square_edge = { pow(10000., 2) / 2f64 }.sqrt();
71        let p_3 = p_1.haversine_destination(0., square_edge);
72        let p_4 = p_3.haversine_destination(90., square_edge);
73        assert_relative_eq!(p_4, p_2, epsilon = 1.0e-6);
74    }
75
76    #[test]
77    fn bearing_zero_is_north() {
78        let p_1 = Point::new(9.177789688110352, 48.776781529534965);
79        let p_2 = p_1.haversine_destination(0., 1000.);
80        assert_relative_eq!(p_1.x(), p_2.x(), epsilon = 1.0e-6);
81        assert!(p_2.y() > p_1.y())
82    }
83}
geo/algorithm/haversine_destination.rs

geo/algorithm/
haversine_destination.rs
