gistools/tools/
orthodrome.rs

1use libm::{atan2, cos, sin, sqrt};
2use s2json::{GetXY, NewXY};
3
4/// # Orthodrome
5///
6/// ## Description
7/// Represents an orthodrome, which is the shortest path between two points on a sphere.
8///
9/// [Learn more here](http://www.movable-type.co.uk/scripts/latlong.html)
10///
11/// ## Usage
12///
13/// The methods you have access to:
14/// - [`Orthodrome::new`]: Create a new Orthodrome
15/// - [`Orthodrome::from_points`]: Create an orthodrome from two points. The points must implement the [`GetXY`] trait
16/// - [`Orthodrome::intermediate_point`]: input t 0->1. Find a point along the orthodrome.
17/// - [`Orthodrome::bearing`]: returns the bearing in degrees between the two points
18/// - [`Orthodrome::distance_to`]: Finds the distance between the two points in kilometers projected normalized (0->1)
19///
20/// The properties you have access to:
21/// - [`Orthodrome::lon1`]: start longitude in radians
22/// - [`Orthodrome::lat1`]: start latitude in radians
23/// - [`Orthodrome::lon2`]: end longitude in radians
24/// - [`Orthodrome::lat2`]: end latitude in radians
25/// - [`Orthodrome::a`]: distance property
26/// - [`Orthodrome::dist`]: distance property
27///
28/// ```rust
29/// use gistools::{geometry::LonLat, tools::Orthodrome};
30///
31/// let ortho = Orthodrome::new(-60., -40., 20., 10.);
32/// assert_eq!(ortho.intermediate_point::<LonLat>(0.), LonLat::new(-59.99999999999999, -40., None));
33/// assert_eq!(
34///     ortho.intermediate_point::<LonLat>(0.2),
35///     LonLat::new(-39.13793657428956, -33.728521975616516, None)
36/// );
37/// ```
38///
39/// ## NOTE
40/// There is no reason to use this outside verbosity. You can create an S1Angle or use the utility functions in LonLat
41///
42/// ## Links
43/// - <http://www.movable-type.co.uk/scripts/latlong.html>
44#[derive(Debug, Clone, Default)]
45pub struct Orthodrome {
46    /// start longitude in radians
47    pub lon1: f64,
48    /// start latitude in radians
49    pub lat1: f64,
50    /// end longitude in radians
51    pub lon2: f64,
52    /// end latitude in radians
53    pub lat2: f64,
54    /// distance property
55    pub a: f64,
56    /// distance property
57    pub dist: f64,
58}
59impl Orthodrome {
60    /// Create an orthodrome
61    pub fn new(start_lon: f64, start_lat: f64, end_lon: f64, end_lat: f64) -> Orthodrome {
62        let lon1 = start_lon.to_radians();
63        let lat1 = start_lat.to_radians();
64        let lon2 = end_lon.to_radians();
65        let lat2 = end_lat.to_radians();
66        let d_lat = lat2 - lat1;
67        let d_lon = lon2 - lon1;
68        let a = sin(d_lat / 2.) * sin(d_lat / 2.)
69            + cos(lat1) * cos(lat2) * sin(d_lon / 2.) * sin(d_lon / 2.);
70        let dist = 2. * atan2(sqrt(a), sqrt(1. - a));
71
72        Orthodrome { lon1, lat1, lon2, lat2, a, dist }
73    }
74
75    /// Create an orthodrome from two points
76    ///
77    /// The points must implement the [`GetXY`] trait
78    pub fn from_points<P1: GetXY, P2: GetXY>(p1: &P1, p2: &P2) -> Orthodrome {
79        Orthodrome::new(p1.x(), p1.y(), p2.x(), p2.y())
80    }
81
82    /// input t 0->1. Find a point along the orthodrome.
83    ///
84    /// ## Parameters
85    /// - `t`: distance along the orthodrome to find
86    ///
87    /// ## Returns
88    /// Returns any Point type that implements the [`NewXY`].
89    /// X is the longitude in degrees.
90    /// Y is the latitude in degrees.
91    pub fn intermediate_point<P: NewXY>(&self, t: f64) -> P {
92        let Self { lon1, lon2, lat1, lat2, dist, .. } = self;
93
94        // check corner cases first
95        if t == 0. {
96            return P::new_xy(lon1.to_degrees(), lat1.to_degrees());
97        } else if t == 1. {
98            return P::new_xy(lon2.to_degrees(), lat2.to_degrees());
99        }
100        // check if points are equal
101        else if lon1 == lon2 && lat1 == lat2 {
102            return P::new_xy(lon1.to_degrees(), lat1.to_degrees());
103        }
104
105        let a = sin((1. - t) * dist) / sin(*dist);
106        let b = sin(t * dist) / sin(*dist);
107
108        let x = a * cos(*lat1) * cos(*lon1) + b * cos(*lat2) * cos(*lon2);
109        let y = a * cos(*lat1) * sin(*lon1) + b * cos(*lat2) * sin(*lon2);
110        let z = a * sin(*lat1) + b * sin(*lat2);
111
112        let lat = atan2(z, sqrt(x * x + y * y));
113        let lon = atan2(y, x);
114
115        P::new_xy(lon.to_degrees(), lat.to_degrees())
116    }
117
118    /// returns the bearing in degrees between the two points
119    pub fn bearing(&self) -> f64 {
120        let Self { lon1, lat1, lon2, lat2, .. } = self;
121
122        let y = sin(lon2 - lon1) * cos(*lat2);
123        let x = cos(*lat1) * sin(*lat2) - sin(*lat1) * cos(*lat2) * cos(lon2 - lon1);
124        let angle_rad = atan2(y, x);
125
126        (angle_rad.to_degrees() + 360.) % 360. // in degrees
127    }
128
129    /// Finds the distance between the two points in kilometers projected normalized (0->1)
130    ///
131    /// ## Returns
132    /// The total distance between the two points
133    pub fn distance_to(&self) -> f64 {
134        2. * atan2(sqrt(self.a), sqrt(1. - self.a))
135    }
136}
gistools/tools/orthodrome.rs

gistools/tools/
orthodrome.rs
