gistools/geometry/wm/
lonlat.rs

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

gistools/geometry/wm/
lonlat.rs
