geo/algorithm/
geodesic_intermediate.rs

1use crate::{CoordFloat, Point};
2use geographiclib_rs::{DirectGeodesic, Geodesic, InverseGeodesic};
3
4/// Returns a new Point along a route between two existing points on an ellipsoidal model of the earth
5
6pub trait GeodesicIntermediate<T: CoordFloat> {
7    /// Returns a new Point along a route between two existing points on an ellipsoidal model of the earth
8    ///
9    /// # Examples
10    ///
11    /// ```
12    /// # #[macro_use] extern crate approx;
13    /// #
14    /// use geo::GeodesicIntermediate;
15    /// use geo::Point;
16    ///
17    /// let p1 = Point::new(10.0, 20.0);
18    /// let p2 = Point::new(125.0, 25.0);
19    /// let i20 = p1.geodesic_intermediate(&p2, 0.2);
20    /// let i50 = p1.geodesic_intermediate(&p2, 0.5);
21    /// let i80 = p1.geodesic_intermediate(&p2, 0.8);
22    /// let i20_should = Point::new(29.842907, 29.951445);
23    /// let i50_should = Point::new(65.879360, 37.722253);
24    /// let i80_should = Point::new(103.556796, 33.506196);
25    /// assert_relative_eq!(i20, i20_should, epsilon = 1.0e-6);
26    /// assert_relative_eq!(i50, i50_should, epsilon = 1.0e-6);
27    /// assert_relative_eq!(i80, i80_should, epsilon = 1.0e-6);
28    /// ```
29
30    fn geodesic_intermediate(&self, other: &Point<T>, f: T) -> Point<T>;
31    fn geodesic_intermediate_fill(
32        &self,
33        other: &Point<T>,
34        max_dist: T,
35        include_ends: bool,
36    ) -> Vec<Point<T>>;
37}
38
39impl GeodesicIntermediate<f64> for Point {
40    fn geodesic_intermediate(&self, other: &Point, f: f64) -> Point {
41        let g = Geodesic::wgs84();
42        let (total_distance, azi1, _azi2, _a12) =
43            g.inverse(self.y(), self.x(), other.y(), other.x());
44        let distance = total_distance * f;
45        let (lat2, lon2) = g.direct(self.y(), self.x(), azi1, distance);
46
47        Point::new(lon2, lat2)
48    }
49
50    fn geodesic_intermediate_fill(
51        &self,
52        other: &Point,
53        max_dist: f64,
54        include_ends: bool,
55    ) -> Vec<Point> {
56        let g = Geodesic::wgs84();
57        let (total_distance, azi1, _azi2, _a12) =
58            g.inverse(self.y(), self.x(), other.y(), other.x());
59
60        if total_distance <= max_dist {
61            return if include_ends {
62                vec![*self, *other]
63            } else {
64                vec![]
65            };
66        }
67
68        let number_of_points = (total_distance / max_dist).ceil();
69        let interval = 1.0 / number_of_points;
70
71        let mut current_step = interval;
72        let mut points = if include_ends { vec![*self] } else { vec![] };
73
74        while current_step < 1.0 {
75            let (lat2, lon2) = g.direct(self.y(), self.x(), azi1, total_distance * current_step);
76            let point = Point::new(lon2, lat2);
77            points.push(point);
78            current_step += interval;
79        }
80
81        if include_ends {
82            points.push(*other);
83        }
84
85        points
86    }
87}
88
89#[cfg(test)]
90mod tests {
91    use super::*;
92    use approx::assert_relative_eq;
93
94    #[test]
95    fn f_is_zero_or_one_test() {
96        let p1 = Point::new(10.0, 20.0);
97        let p2 = Point::new(15.0, 25.0);
98        let i0 = p1.geodesic_intermediate(&p2, 0.0);
99        let i100 = p1.geodesic_intermediate(&p2, 1.0);
100        assert_relative_eq!(i0, p1, epsilon = 1.0e-6);
101        assert_relative_eq!(i100, p2, epsilon = 1.0e-6);
102    }
103
104    #[test]
105    fn various_f_values_test() {
106        let p1 = Point::new(10.0, 20.0);
107        let p2 = Point::new(125.0, 25.0);
108        let i20 = p1.geodesic_intermediate(&p2, 0.2);
109        let i50 = p1.geodesic_intermediate(&p2, 0.5);
110        let i80 = p1.geodesic_intermediate(&p2, 0.8);
111        let i20_should = Point::new(29.842907, 29.951445);
112        let i50_should = Point::new(65.879360, 37.722253);
113        let i80_should = Point::new(103.556796, 33.506196);
114        assert_relative_eq!(i20, i20_should, epsilon = 1.0e-6);
115        assert_relative_eq!(i50, i50_should, epsilon = 1.0e-6);
116        assert_relative_eq!(i80, i80_should, epsilon = 1.0e-6);
117    }
118
119    #[test]
120    fn should_add_i50_test() {
121        let p1 = Point::new(30.0, 40.0);
122        let p2 = Point::new(40.0, 50.0);
123        let max_dist = 1000000.0; // meters
124        let include_ends = true;
125        let i50 = p1.geodesic_intermediate(&p2, 0.5);
126        let route = p1.geodesic_intermediate_fill(&p2, max_dist, include_ends);
127        assert_eq!(route, vec![p1, i50, p2]);
128    }
129}
geo/algorithm/geodesic_intermediate.rs

geo/algorithm/
geodesic_intermediate.rs
