gistools/geometry/ll/
mod.rs

1use super::{S1Angle, S2CellId, S2Point};
2use core::{
3    cmp::Ordering,
4    ops::{Add, Div, Mul, Neg, Sub},
5};
6use libm::{asin, atan2, cos, fmin, sin, sqrt};
7use s2json::{GetM, GetXY, MValue, NewXY, VectorPoint};
8
9/// # Longitude-Latitude Point container
10///
11/// ## Description
12/// This class represents a point on the unit sphere as a pair
13/// of latitude-longitude coordinates.  Like the rest of the "geometry"
14/// package, the intent is to represent spherical geometry as a mathematical
15/// abstraction, so functions that are specifically related to the Earth's
16/// geometry (e.g. easting/northing conversions) should be put elsewhere.
17///
18/// This struct implements the [`GetXY`] and [`GetM`] traits.
19///
20/// ## Usage
21///
22/// Methods that are available:
23/// - [`LonLat::new`]: Create a new LonLat
24/// - [`LonLat::take`]: Take ownership of the underlying VectorPoint
25/// - [`LonLat::lon`]: Return the longitude in degrees
26/// - [`LonLat::lat`]: Return the latitude in degrees
27/// - [`LonLat::from_s2cellid`]: Convert a [`S2CellId`] to an LonLat
28/// - [`LonLat::from_s2_point`]: Convert a [`S2Point`] to an LonLat
29/// - [`LonLat::normalize`]: Normalize the coordinates to the range [-180, 180] and [-90, 90] deg.
30/// - [`LonLat::coords`]: Return the latitude or longitude coordinates in degrees.
31/// - [`LonLat::to_angles`]: Return the latitude and longitude coordinates in radians.
32/// - [`LonLat::is_valid`]: Return true if the latitude is between -90 and 90 degrees inclusive and the longitude is between -180 and 180 degrees inclusive.
33/// - [`LonLat::to_point`]: Converts an LonLat to the equivalent unit-length vector.
34/// - [`LonLat::to_point_gl`]: An alternative to [`LonLat::to_point`] that returns a GPU compatible vector.
35/// - [`LonLat::get_distance`]:Returns the distance (measured along the surface of the sphere) to the given LonLat.
36/// - [`LonLat::get_bearing`]: Returns the bearing from the first point to the second point.
37#[derive(Clone, PartialEq, Debug)]
38pub struct LonLat<M: Clone + Default = MValue>(pub VectorPoint<M>);
39impl GetXY for LonLat {
40    fn x(&self) -> f64 {
41        self.0.x
42    }
43    fn y(&self) -> f64 {
44        self.0.y
45    }
46}
47impl<M: Clone + Default> GetM<M> for LonLat<M> {
48    fn m(&self) -> Option<&M> {
49        self.0.m.as_ref()
50    }
51}
52impl NewXY for LonLat {
53    fn new_xy(x: f64, y: f64) -> LonLat {
54        LonLat(VectorPoint::new(x, y, None, None))
55    }
56}
57impl<M: Clone + Default> LonLat<M> {
58    /// Build a new LonLat
59    pub fn new(lon: f64, lat: f64, m: Option<M>) -> LonLat<M> {
60        LonLat(VectorPoint::new(lon, lat, None, m))
61    }
62
63    /// Take the underlying VectorPoint
64    pub fn take(&mut self) -> VectorPoint<M> {
65        core::mem::take(&mut self.0)
66    }
67
68    /// Return the longitude in degrees
69    #[inline]
70    pub fn lon(&self) -> f64 {
71        self.0.x
72    }
73
74    /// Return the latitude in degrees
75    #[inline]
76    pub fn lat(&self) -> f64 {
77        self.0.y
78    }
79
80    /// Convert a S2CellId to an LonLat
81    pub fn from_s2cellid(cellid: S2CellId) -> Self {
82        let p = cellid.to_point();
83        Self::from_s2_point(&p)
84    }
85
86    /// Convert a direction vector (not necessarily unit length) to an LonLat.
87    pub fn from_s2_point(p: &S2Point) -> Self {
88        let lon_rad = atan2(p.y + 0.0, p.x + 0.0);
89        let lat_rad = atan2(p.z, sqrt(p.x * p.x + p.y * p.y));
90        let ll = LonLat::new((lon_rad).to_degrees(), (lat_rad).to_degrees(), None);
91        debug_assert!(ll.is_valid(), "from_s2_point() called on invalid LonLat");
92
93        ll
94    }
95
96    /// Normalize the coordinates to the range [-180, 180] and [-90, 90] deg.
97    pub fn normalize(&mut self) {
98        let mut lon = self.lon();
99        let mut lat = self.lat();
100        // Normalize longitude using modulo
101        lon = ((((lon + 180.) % 360.) + 360.) % 360.) - 180.;
102        // Clamp latitude between -90 and 90
103        lat = lat.clamp(-90., 90.);
104
105        self.0.x = lon;
106        self.0.y = lat;
107    }
108
109    /// Return the latitude or longitude coordinates in degrees.
110    pub fn coords(self) -> (f64, f64) {
111        (self.lon(), self.lat())
112    }
113
114    /// Return the latitude and longitude coordinates in radians.
115    pub fn to_angles(&self) -> (S1Angle, S1Angle) {
116        (S1Angle::from_degrees(self.lon()), S1Angle::from_degrees(self.lat()))
117    }
118
119    /// Return true if the latitude is between -90 and 90 degrees inclusive
120    /// and the longitude is between -180 and 180 degrees inclusive.
121    pub fn is_valid(&self) -> bool {
122        // OLD RAD IMPL:
123        // fabs(self.lat()) <= (PI / 2.0) && fabs(self.lon()) <= PI
124        let lat = self.lat();
125        let lon = self.lon();
126
127        (-90.0..=90.0).contains(&lat) && (-180.0..=180.0).contains(&lon)
128    }
129
130    /// Converts an LonLat to the equivalent unit-length vector. Unnormalized
131    /// values (see Normalize()) are wrapped around the sphere as would be expected
132    /// based on their definition as spherical angles.  So for example the
133    /// following pairs yield equivalent points (modulo numerical error):
134    ///     (90.5, 10) =~ (89.5, -170)
135    ///     (a, b) =~ (a + 360 * n, b)
136    /// The maximum error in the result is 1.5 * DBL_EPSILON.  (This does not
137    /// include the error of converting degrees, E5, E6, or E7 to radians.)
138    ///
139    /// Can be used just like an S2Point constructor.  For example:
140    ///   S2Cap cap;
141    ///   cap.AddPoint(S2Point(latlon));
142    pub fn to_point(&self) -> S2Point {
143        debug_assert!(self.is_valid(), "to_point() called on invalid LonLat");
144        let lon: f64 = (self.lon()).to_radians();
145        let lat: f64 = (self.lat()).to_radians();
146        S2Point::new(
147            cos(lat) * cos(lon), // x
148            cos(lat) * sin(lon), // y
149            sin(lat),            // z
150        )
151    }
152
153    /// An alternative to to_point() that returns a GPU compatible vector.
154    pub fn to_point_gl(&self) -> S2Point {
155        debug_assert!(self.is_valid(), "to_point_gl() called on invalid LonLat");
156        let lon: f64 = (self.lon()).to_radians();
157        let lat: f64 = (self.lat()).to_radians();
158        S2Point::new(
159            cos(lat) * sin(lon), // y
160            sin(lat),            // z
161            cos(lat) * cos(lon), // x
162        )
163    }
164
165    /// Returns the distance (measured along the surface of the sphere) to the
166    /// given LonLat, implemented using the Haversine formula.  This is
167    /// equivalent to
168    ///
169    ///   S1Angle(ToPoint(), o.ToPoint())
170    ///
171    /// except that this function is slightly faster, and is also somewhat less
172    /// accurate for distances approaching 180 degrees (see s1angle.h for
173    /// details).  Both LngLats must be normalized.
174    pub fn get_distance<M2: Clone + Default>(&self, b: &LonLat<M2>) -> f64 {
175        // This implements the Haversine formula, which is numerically stable for
176        // small distances but only gets about 8 digits of precision for very large
177        // distances (e.g. antipodal points).  Note that 8 digits is still accurate
178        // to within about 10cm for a sphere the size of the Earth.
179        //
180        // This could be fixed with another sin() and cos() below, but at that point
181        // you might as well just convert both arguments to S2Points and compute the
182        // distance that way (which gives about 15 digits of accuracy for all
183        // distances).
184        debug_assert!(self.is_valid(), "get_bearing() called on invalid LonLat (self)");
185        debug_assert!(b.is_valid(), "get_bearing() called on invalid LonLat (b)");
186
187        let lat1 = self.lat().to_radians();
188        let lat2 = b.lat().to_radians();
189        let lon1 = self.lon().to_radians();
190        let lon2 = b.lon().to_radians();
191        let dlat = sin(0.5 * (lat2 - lat1));
192        let dlon = sin(0.5 * (lon2 - lon1));
193        let x = dlat * dlat + dlon * dlon * cos(lat1) * cos(lat2);
194        2. * asin(sqrt(fmin(1., x)))
195    }
196
197    /// Returns the bearing from the first point to the second point.
198    pub fn get_bearing<M2: Clone + Default>(&self, b: &LonLat<M2>) -> f64 {
199        debug_assert!(self.is_valid(), "get_bearing() called on invalid LonLat (self)");
200        debug_assert!(b.is_valid(), "get_bearing() called on invalid LonLat (self)");
201        let lat1 = self.lat().to_radians();
202        let lat2 = b.lat().to_radians();
203        let lon1 = self.lon().to_radians();
204        let lon2 = b.lon().to_radians();
205        let y = sin(lon2 - lon1) * cos(lat2);
206        let x = cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(lon2 - lon1);
207
208        (atan2(y, x).to_degrees() + 360.) % 360.
209    }
210}
211impl<M: Clone + Default> From<S2CellId> for LonLat<M> {
212    fn from(c: S2CellId) -> Self {
213        LonLat::from_s2cellid(c)
214    }
215}
216impl<M: Clone + Default> From<&S2Point> for LonLat<M> {
217    fn from(p: &S2Point) -> Self {
218        LonLat::from_s2_point(p)
219    }
220}
221impl<M1: Clone + Default, M2: Clone + Default> Add<LonLat<M2>> for LonLat<M1> {
222    type Output = LonLat<M1>;
223    fn add(self, rhs: LonLat<M2>) -> Self::Output {
224        LonLat::new(self.lon() + rhs.lon(), self.lat() + rhs.lat(), self.0.m)
225    }
226}
227impl<M1: Clone + Default, M2: Clone + Default> Sub<LonLat<M2>> for LonLat<M1> {
228    type Output = LonLat<M1>;
229    fn sub(self, rhs: LonLat<M2>) -> Self::Output {
230        LonLat::new(self.lon() - rhs.lon(), self.lat() - rhs.lat(), self.0.m)
231    }
232}
233impl<M1: Clone + Default, M2: Clone + Default> Mul<LonLat<M2>> for LonLat<M1> {
234    type Output = LonLat<M1>;
235    fn mul(self, rhs: LonLat<M2>) -> Self::Output {
236        LonLat::new(self.lon() * rhs.lon(), self.lat() * rhs.lat(), self.0.m)
237    }
238}
239impl<M1: Clone + Default, M2: Clone + Default> Div<LonLat<M2>> for LonLat<M1> {
240    type Output = LonLat<M1>;
241    fn div(self, rhs: LonLat<M2>) -> Self::Output {
242        LonLat::new(self.lon() / rhs.lon(), self.lat() / rhs.lat(), self.0.m)
243    }
244}
245impl Neg for LonLat {
246    type Output = LonLat;
247    fn neg(self) -> Self::Output {
248        LonLat::new(-self.lon(), -self.lat(), self.0.m)
249    }
250}
251impl Eq for LonLat {}
252impl Ord for LonLat {
253    fn cmp(&self, other: &Self) -> Ordering {
254        match self.lon().partial_cmp(&other.lon()) {
255            Some(Ordering::Equal) => {}
256            other => return other.unwrap(), // Handle cases where `lon` comparison is not equal
257        }
258        match self.lat().partial_cmp(&other.lat()) {
259            Some(order) => order,
260            None => Ordering::Equal, // This handles the NaN case safely
261        }
262    }
263}
264impl PartialOrd for LonLat {
265    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
266        Some(self.cmp(other))
267    }
268}
gistools/geometry/ll/mod.rs

gistools/geometry/ll/
mod.rs
