geo/algorithm/
haversine_intermediate.rs

1use crate::{CoordFloat, Point, MEAN_EARTH_RADIUS};
2use num_traits::FromPrimitive;
3
4/// Returns a new Point along a great circle route between two existing points
5
6pub trait HaversineIntermediate<T: CoordFloat> {
7    /// Returns a new Point along a great circle route between two existing points.
8    ///
9    /// # Examples
10    ///
11    /// ```
12    /// # #[macro_use] extern crate approx;
13    /// #
14    /// use geo::HaversineIntermediate;
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.haversine_intermediate(&p2, 0.2);
20    /// let i50 = p1.haversine_intermediate(&p2, 0.5);
21    /// let i80 = p1.haversine_intermediate(&p2, 0.8);
22    /// let i20_should = Point::new(29.8, 29.9);
23    /// let i50_should = Point::new(65.8, 37.7);
24    /// let i80_should = Point::new(103.5, 33.5);
25    /// assert_relative_eq!(i20.x(), i20_should.x(), epsilon = 0.2);
26    /// assert_relative_eq!(i20.y(), i20_should.y(), epsilon = 0.2);
27    /// assert_relative_eq!(i50.x(), i50_should.x(), epsilon = 0.2);
28    /// assert_relative_eq!(i50.y(), i50_should.y(), epsilon = 0.2);
29    /// assert_relative_eq!(i80.x(), i80_should.x(), epsilon = 0.2);
30    /// assert_relative_eq!(i80.y(), i80_should.y(), epsilon = 0.2);
31    /// ```
32
33    fn haversine_intermediate(&self, other: &Point<T>, f: T) -> Point<T>;
34    fn haversine_intermediate_fill(
35        &self,
36        other: &Point<T>,
37        max_dist: T,
38        include_ends: bool,
39    ) -> Vec<Point<T>>;
40}
41
42impl<T> HaversineIntermediate<T> for Point<T>
43where
44    T: CoordFloat + FromPrimitive,
45{
46    fn haversine_intermediate(&self, other: &Point<T>, f: T) -> Point<T> {
47        let params = get_params(self, other);
48        get_point(&params, f)
49    }
50
51    fn haversine_intermediate_fill(
52        &self,
53        other: &Point<T>,
54        max_dist: T,
55        include_ends: bool,
56    ) -> Vec<Point<T>> {
57        let params = get_params(self, other);
58        let HaversineParams { d, .. } = params;
59
60        let total_distance = d * T::from(MEAN_EARTH_RADIUS).unwrap();
61
62        if total_distance <= max_dist {
63            return if include_ends {
64                vec![*self, *other]
65            } else {
66                vec![]
67            };
68        }
69
70        let number_of_points = (total_distance / max_dist).ceil();
71        let interval = T::one() / number_of_points;
72
73        let mut current_step = interval;
74        let mut points = if include_ends { vec![*self] } else { vec![] };
75
76        while current_step < T::one() {
77            let point = get_point(&params, current_step);
78            points.push(point);
79            current_step = current_step + interval;
80        }
81
82        if include_ends {
83            points.push(*other);
84        }
85
86        points
87    }
88}
89
90#[allow(clippy::many_single_char_names)]
91struct HaversineParams<T: num_traits::Float + FromPrimitive> {
92    d: T,
93    n: T,
94    o: T,
95    p: T,
96    q: T,
97    r: T,
98    s: T,
99}
100
101#[allow(clippy::many_single_char_names)]
102fn get_point<T: CoordFloat + FromPrimitive>(params: &HaversineParams<T>, f: T) -> Point<T> {
103    let one = T::one();
104
105    let HaversineParams {
106        d,
107        n,
108        o,
109        p,
110        q,
111        r,
112        s,
113    } = *params;
114
115    let a = ((one - f) * d).sin() / d.sin();
116    let b = (f * d).sin() / d.sin();
117
118    let x = a * n + b * o;
119    let y = a * p + b * q;
120    let z = a * r + b * s;
121
122    let lat = z.atan2(x.hypot(y));
123    let lon = y.atan2(x);
124
125    Point::new(lon.to_degrees(), lat.to_degrees())
126}
127
128#[allow(clippy::many_single_char_names)]
129fn get_params<T: CoordFloat + FromPrimitive>(p1: &Point<T>, p2: &Point<T>) -> HaversineParams<T> {
130    let one = T::one();
131    let two = one + one;
132
133    let lat1 = p1.y().to_radians();
134    let lon1 = p1.x().to_radians();
135    let lat2 = p2.y().to_radians();
136    let lon2 = p2.x().to_radians();
137
138    let (lat1_sin, lat1_cos) = lat1.sin_cos();
139    let (lat2_sin, lat2_cos) = lat2.sin_cos();
140    let (lon1_sin, lon1_cos) = lon1.sin_cos();
141    let (lon2_sin, lon2_cos) = lon2.sin_cos();
142
143    let m = lat1_cos * lat2_cos;
144
145    let n = lat1_cos * lon1_cos;
146    let o = lat2_cos * lon2_cos;
147    let p = lat1_cos * lon1_sin;
148    let q = lat2_cos * lon2_sin;
149
150    let k = (((lat1 - lat2) / two).sin().powi(2) + m * ((lon1 - lon2) / two).sin().powi(2)).sqrt();
151
152    let d = two * k.asin();
153
154    HaversineParams {
155        d,
156        n,
157        o,
158        p,
159        q,
160        r: lat1_sin,
161        s: lat2_sin,
162    }
163}
164
165#[cfg(test)]
166mod test {
167    use super::*;
168    use crate::HaversineIntermediate;
169
170    #[test]
171    fn f_is_zero_or_one_test() {
172        let p1 = Point::new(10.0, 20.0);
173        let p2 = Point::new(15.0, 25.0);
174        let i0 = p1.haversine_intermediate(&p2, 0.0);
175        let i100 = p1.haversine_intermediate(&p2, 1.0);
176        assert_relative_eq!(i0.x(), p1.x(), epsilon = 1.0e-6);
177        assert_relative_eq!(i0.y(), p1.y(), epsilon = 1.0e-6);
178        assert_relative_eq!(i100.x(), p2.x(), epsilon = 1.0e-6);
179        assert_relative_eq!(i100.y(), p2.y(), epsilon = 1.0e-6);
180    }
181
182    #[test]
183    fn various_f_values_test() {
184        let p1 = Point::new(10.0, 20.0);
185        let p2 = Point::new(125.0, 25.0);
186        let i20 = p1.haversine_intermediate(&p2, 0.2);
187        let i50 = p1.haversine_intermediate(&p2, 0.5);
188        let i80 = p1.haversine_intermediate(&p2, 0.8);
189        let i20_should = Point::new(29.83519, 29.94841);
190        let i50_should = Point::new(65.87471, 37.72201);
191        let i80_should = Point::new(103.56036, 33.50518);
192        assert_relative_eq!(i20.x(), i20_should.x(), epsilon = 0.2);
193        assert_relative_eq!(i20.y(), i20_should.y(), epsilon = 0.2);
194        assert_relative_eq!(i50.x(), i50_should.x(), epsilon = 0.2);
195        assert_relative_eq!(i50.y(), i50_should.y(), epsilon = 0.2);
196        assert_relative_eq!(i80.x(), i80_should.x(), epsilon = 0.2);
197        assert_relative_eq!(i80.y(), i80_should.y(), epsilon = 0.2);
198    }
199
200    #[test]
201    fn should_be_north_pole_test() {
202        let p1 = Point::new(0.0, 10.0);
203        let p2 = Point::new(180.0, 10.0);
204        let i50 = p1.haversine_intermediate(&p2, 0.5);
205        let i50_should = Point::new(90.0, 90.0);
206        assert_relative_eq!(i50.x(), i50_should.x(), epsilon = 1.0e-6);
207        assert_relative_eq!(i50.y(), i50_should.y(), epsilon = 1.0e-6);
208    }
209
210    #[test]
211    fn should_be_start_end_test() {
212        let p1 = Point::new(30.0, 40.0);
213        let p2 = Point::new(40.0, 50.0);
214        let max_dist = 1500000.0; // meters
215        let include_ends = true;
216        let route = p1.haversine_intermediate_fill(&p2, max_dist, include_ends);
217        assert_eq!(route, vec![p1, p2]);
218    }
219
220    #[test]
221    fn should_add_i50_test() {
222        let p1 = Point::new(30.0, 40.0);
223        let p2 = Point::new(40.0, 50.0);
224        let max_dist = 1000000.0; // meters
225        let include_ends = true;
226        let i50 = p1.clone().haversine_intermediate(&p2, 0.5);
227        let route = p1.haversine_intermediate_fill(&p2, max_dist, include_ends);
228        assert_eq!(route, vec![p1, i50, p2]);
229    }
230
231    #[test]
232    fn should_add_i25_i50_i75_test() {
233        let p1 = Point::new(30.0, 40.0);
234        let p2 = Point::new(40.0, 50.0);
235        let max_dist = 400000.0; // meters
236        let include_ends = true;
237        let i25 = p1.clone().haversine_intermediate(&p2, 0.25);
238        let i50 = p1.clone().haversine_intermediate(&p2, 0.5);
239        let i75 = p1.clone().haversine_intermediate(&p2, 0.75);
240        let route = p1.haversine_intermediate_fill(&p2, max_dist, include_ends);
241        assert_eq!(route, vec![p1, i25, i50, i75, p2]);
242    }
243}
geo/algorithm/haversine_intermediate.rs

geo/algorithm/
haversine_intermediate.rs
