geo/algorithm/rhumb/
mod.rs

1//! This module provides rhumb-line (a.k.a. loxodrome) geometry operations.
2//! The distance, destination, and bearing implementations are adapted in part
3//! from their equivalents in [Turf.js](https://turfjs.org/), which in turn are
4//! adapted from the Movable Type
5//! [spherical geodesy tools](https://www.movable-type.co.uk/scripts/latlong.html).
6//! Turf.js is copyright its authors and the geodesy tools are copyright Chris
7//! Veness; both are available under an MIT license.
8
9use crate::{CoordFloat, Point, point, utils::normalize_longitude};
10use num_traits::FromPrimitive;
11
12mod distance;
13#[allow(deprecated)]
14pub use distance::RhumbDistance;
15
16mod bearing;
17#[allow(deprecated)]
18pub use bearing::RhumbBearing;
19
20mod destination;
21#[allow(deprecated)]
22pub use destination::RhumbDestination;
23
24mod intermediate;
25#[allow(deprecated)]
26pub use intermediate::RhumbIntermediate;
27
28mod length;
29#[allow(deprecated)]
30pub use length::RhumbLength;
31
32pub(crate) struct RhumbCalculations<T: CoordFloat + FromPrimitive> {
33    from: Point<T>,
34    to: Point<T>,
35    phi1: T,
36    delta_lambda: T,
37    delta_phi: T,
38    delta_psi: T,
39}
40
41impl<T: CoordFloat + FromPrimitive> RhumbCalculations<T> {
42    pub(crate) fn new(from: &Point<T>, to: &Point<T>) -> Self {
43        let pi = T::from(std::f64::consts::PI).unwrap();
44        let two = T::one() + T::one();
45        let four = two + two;
46
47        let phi1 = from.y().to_radians();
48        let phi2 = to.y().to_radians();
49        let mut delta_lambda = (to.x() - from.x()).to_radians();
50        // if delta_lambda is over 180° take shorter rhumb line across the anti-meridian:
51        if delta_lambda > pi {
52            delta_lambda = delta_lambda - (two * pi);
53        }
54        if delta_lambda < -pi {
55            delta_lambda = delta_lambda + (two * pi);
56        }
57
58        let delta_psi = ((phi2 / two + pi / four).tan() / (phi1 / two + pi / four).tan()).ln();
59        let delta_phi = phi2 - phi1;
60
61        RhumbCalculations {
62            from: *from,
63            to: *to,
64            phi1,
65            delta_lambda,
66            delta_phi,
67            delta_psi,
68        }
69    }
70
71    pub(crate) fn theta(&self) -> T {
72        self.delta_lambda.atan2(self.delta_psi)
73    }
74
75    pub(crate) fn delta(&self) -> T {
76        let threshold = T::from(10.0e-12).unwrap();
77        let q = if self.delta_psi.abs() > threshold {
78            self.delta_phi / self.delta_psi
79        } else {
80            self.phi1.cos()
81        };
82
83        (self.delta_phi * self.delta_phi + q * q * self.delta_lambda * self.delta_lambda).sqrt()
84    }
85
86    pub(crate) fn intermediate(&self, fraction: T) -> Point<T> {
87        let delta = fraction * self.delta();
88        let theta = self.theta();
89        let lambda1 = self.from.x().to_radians();
90        calculate_destination(delta, lambda1, self.phi1, theta)
91    }
92
93    pub(crate) fn intermediate_fill(&self, max_delta: T, include_ends: bool) -> Vec<Point<T>> {
94        let theta = self.theta();
95        let lambda1 = self.from.x().to_radians();
96
97        let total_delta = self.delta();
98
99        if total_delta <= max_delta {
100            return if include_ends {
101                vec![self.from, self.to]
102            } else {
103                vec![]
104            };
105        }
106
107        let number_of_points = (total_delta / max_delta).ceil();
108        let interval = T::one() / number_of_points;
109
110        let mut current_step = interval;
111        let mut points = if include_ends {
112            vec![self.from]
113        } else {
114            vec![]
115        };
116
117        while current_step < T::one() {
118            let delta = total_delta * current_step;
119            let point = calculate_destination(delta, lambda1, self.phi1, theta);
120            points.push(point);
121            current_step = current_step + interval;
122        }
123
124        if include_ends {
125            points.push(self.to);
126        }
127
128        points
129    }
130}
131
132pub(crate) fn calculate_destination<T: CoordFloat + FromPrimitive>(
133    delta: T,
134    lambda1: T,
135    phi1: T,
136    theta: T,
137) -> Point<T> {
138    let pi = T::from(std::f64::consts::PI).unwrap();
139    let two = T::one() + T::one();
140    let four = two + two;
141    let threshold = T::from(10.0e-12).unwrap();
142
143    let delta_phi = delta * theta.cos();
144    let mut phi2 = phi1 + delta_phi;
145
146    // check for some daft bugger going past the pole, normalise latitude if so
147    if phi2.abs() > pi / two {
148        phi2 = if phi2 > T::zero() {
149            pi - phi2
150        } else {
151            -pi - phi2
152        };
153    }
154
155    let delta_psi = ((phi2 / two + pi / four).tan() / (phi1 / two + pi / four).tan()).ln();
156    // E-W course becomes ill-conditioned with 0/0
157    let q = if delta_psi.abs() > threshold {
158        delta_phi / delta_psi
159    } else {
160        phi1.cos()
161    };
162
163    let delta_lambda = (delta * theta.sin()) / q;
164    let lambda2 = lambda1 + delta_lambda;
165
166    point! {
167        x: normalize_longitude(lambda2.to_degrees()),
168        y: phi2.to_degrees(),
169    }
170}
geo/algorithm/rhumb/mod.rs

geo/algorithm/rhumb/
mod.rs
