geometry_strategy/geographic/
rhumb.rs

1//! Constant-bearing rhumb-line measurements on an angular coordinate system.
2//!
3//! Boost.Geometry has no rhumb-line strategy. The closed-form spherical
4//! formulas follow Ed Williams' *Aviation Formulary* and the isometric-latitude
5//! treatment described by Bowring. The strategy is shared by spherical and
6//! geographic coordinate families; geographic use is a spherical mean-radius
7//! approximation rather than an ellipsoidal geodesic.
8
9#[cfg(not(feature = "std"))]
10use geometry_coords::math::Float;
11use geometry_cs::{AngleUnit, CoordinateSystem, GeographicFamily, SphericalFamily};
12use geometry_model::Point2D;
13use geometry_trait::{Linestring, Point};
14
15use crate::azimuth::AzimuthStrategy;
16use crate::destination::DestinationStrategy;
17use crate::distance::DistanceStrategy;
18use crate::length::LengthStrategy;
19use crate::normalise::{HasAngularUnits, lonlat_radians};
20
21/// Coordinate-system families for which a loxodrome is defined.
22#[doc(hidden)]
23pub trait RhumbFamily {}
24
25impl RhumbFamily for SphericalFamily {}
26impl RhumbFamily for GeographicFamily {}
27
28/// Spherical rhumb-line metric with a configurable radius.
29#[derive(Debug, Clone, Copy)]
30pub struct Rhumb {
31    /// Sphere radius in distance units.
32    pub radius: f64,
33}
34
35impl Rhumb {
36    /// IUGG mean Earth radius in metres.
37    pub const EARTH: Self = Self {
38        radius: 6_371_008.8,
39    };
40
41    /// Unit sphere, returning angular distance in radians.
42    pub const UNIT: Self = Self { radius: 1.0 };
43
44    /// Construct a rhumb metric with an application-defined radius.
45    #[must_use]
46    pub const fn with_radius(radius: f64) -> Self {
47        Self { radius }
48    }
49}
50
51impl Default for Rhumb {
52    fn default() -> Self {
53        Self::EARTH
54    }
55}
56
57impl<P1, P2> DistanceStrategy<P1, P2> for Rhumb
58where
59    P1: Point<Scalar = f64>,
60    P2: Point<Scalar = f64, Cs = P1::Cs>,
61    P1::Cs: HasAngularUnits + CoordinateSystem,
62    <P1::Cs as CoordinateSystem>::Family: RhumbFamily,
63{
64    type Out = f64;
65    type Comparable = Self;
66
67    fn distance(&self, first: &P1, second: &P2) -> Self::Out {
68        rhumb_inverse(first, second).0 * self.radius
69    }
70
71    fn comparable(&self) -> Self::Comparable {
72        *self
73    }
74}
75
76impl<P1, P2> AzimuthStrategy<P1, P2> for Rhumb
77where
78    P1: Point<Scalar = f64>,
79    P2: Point<Scalar = f64, Cs = P1::Cs>,
80    P1::Cs: HasAngularUnits + CoordinateSystem,
81    <P1::Cs as CoordinateSystem>::Family: RhumbFamily,
82{
83    type Out = f64;
84
85    fn azimuth(&self, first: &P1, second: &P2) -> Self::Out {
86        rhumb_inverse(first, second).1
87    }
88}
89
90impl<P> DestinationStrategy<P> for Rhumb
91where
92    P: Point<Scalar = f64>,
93    P::Cs: HasAngularUnits + CoordinateSystem,
94    <P::Cs as CoordinateSystem>::Family: RhumbFamily,
95{
96    type Output = Point2D<f64, P::Cs>;
97
98    fn destination(&self, origin: &P, bearing: f64, distance: f64) -> Self::Output {
99        type Units<P> = <<P as Point>::Cs as HasAngularUnits>::Units;
100        let (longitude1, latitude1) = lonlat_radians(origin);
101        let angular_distance = distance / self.radius;
102        let delta_latitude = angular_distance * bearing.cos();
103        let latitude2 = reflect_latitude(latitude1 + delta_latitude);
104        let delta_psi = isometric_latitude(latitude2) - isometric_latitude(latitude1);
105        let q = meridional_scale(delta_latitude, delta_psi, latitude1);
106        let delta_longitude = if q.abs() <= f64::EPSILON {
107            0.0
108        } else {
109            angular_distance * bearing.sin() / q
110        };
111        let longitude2 = normalize_longitude(longitude1 + delta_longitude);
112        Point2D::new(
113            Units::<P>::from_radians(longitude2),
114            Units::<P>::from_radians(latitude2),
115        )
116    }
117}
118
119impl<L> LengthStrategy<L> for Rhumb
120where
121    L: Linestring,
122    L::Point: Point<Scalar = f64>,
123    <L::Point as Point>::Cs: HasAngularUnits + CoordinateSystem,
124    <<L::Point as Point>::Cs as CoordinateSystem>::Family: RhumbFamily,
125{
126    type Out = f64;
127
128    fn length(&self, line: &L) -> Self::Out {
129        let points = line.points();
130        points
131            .clone()
132            .zip(points.skip(1))
133            .map(|(first, second)| {
134                <Self as DistanceStrategy<L::Point, L::Point>>::distance(self, first, second)
135            })
136            .sum()
137    }
138}
139
140fn rhumb_inverse<P1, P2>(first: &P1, second: &P2) -> (f64, f64)
141where
142    P1: Point<Scalar = f64>,
143    P2: Point<Scalar = f64, Cs = P1::Cs>,
144    P1::Cs: HasAngularUnits,
145{
146    let (longitude1, latitude1) = lonlat_radians(first);
147    let (longitude2, latitude2) = lonlat_radians(second);
148    let delta_latitude = latitude2 - latitude1;
149    let delta_longitude = normalize_delta(longitude2 - longitude1);
150    let delta_psi = isometric_latitude(latitude2) - isometric_latitude(latitude1);
151    let q = meridional_scale(delta_latitude, delta_psi, latitude1);
152    let angular_distance = delta_latitude.hypot(q * delta_longitude);
153    let azimuth = delta_longitude
154        .atan2(delta_psi)
155        .rem_euclid(core::f64::consts::TAU);
156    (angular_distance, azimuth)
157}
158
159fn isometric_latitude(latitude: f64) -> f64 {
160    (core::f64::consts::FRAC_PI_4 + latitude / 2.0).tan().ln()
161}
162
163fn meridional_scale(delta_latitude: f64, delta_psi: f64, latitude: f64) -> f64 {
164    if delta_psi.abs() > 1e-12 {
165        delta_latitude / delta_psi
166    } else {
167        latitude.cos()
168    }
169}
170
171fn normalize_delta(delta: f64) -> f64 {
172    (delta + core::f64::consts::PI).rem_euclid(core::f64::consts::TAU) - core::f64::consts::PI
173}
174
175fn normalize_longitude(longitude: f64) -> f64 {
176    normalize_delta(longitude)
177}
178
179fn reflect_latitude(latitude: f64) -> f64 {
180    let latitude = (latitude + core::f64::consts::PI).rem_euclid(core::f64::consts::TAU)
181        - core::f64::consts::PI;
182    if latitude > core::f64::consts::FRAC_PI_2 {
183        core::f64::consts::PI - latitude
184    } else if latitude < -core::f64::consts::FRAC_PI_2 {
185        -core::f64::consts::PI - latitude
186    } else {
187        latitude
188    }
189}
190
191#[cfg(test)]
192mod tests {
193    use geometry_cs::{Degree, Spherical};
194    use geometry_model::{Linestring, Point2D};
195    use geometry_trait::Point as _;
196
197    use super::*;
198
199    #[test]
200    fn equatorial_degree_has_expected_measurements() {
201        type P = Point2D<f64, Spherical<Degree>>;
202        let start = P::new(0.0, 0.0);
203        let east = P::new(1.0, 0.0);
204        let distance = Rhumb::EARTH.distance(&start, &east);
205        assert!((distance - 111_195.080_233_532_9).abs() < 1e-6);
206        assert!((Rhumb::EARTH.azimuth(&start, &east) - core::f64::consts::FRAC_PI_2).abs() < 1e-12);
207        let destination = Rhumb::EARTH.destination(&start, core::f64::consts::FRAC_PI_2, distance);
208        assert!((destination.get::<0>() - 1.0).abs() < 1e-10);
209
210        let line = Linestring::from_vec(alloc::vec![start, east, P::new(2.0, 0.0)]);
211        assert!((Rhumb::EARTH.length(&line) - 2.0 * distance).abs() < 1e-6);
212    }
213}
geometry_strategy/geographic/rhumb.rs

geometry_strategy/geographic/
rhumb.rs
