plane_projection/
lib.rs

1// Values that define WGS84 ellipsoid model of the Earth
2const EQUATORIAL_RADIUS: f64 = 6378137.0;
3const FLATTENING: f64 = 1.0 / 298.257223563;
4const SQUARED_ECCENTRICITY: f64 = FLATTENING * (2.0 - FLATTENING);
5
6/// A coordinate in (latitude, longitude) format
7pub type LatLon = (f64, f64);
8
9/// A plane projection, useful for blazingly fast approximate distance calculations.
10/// Based on WGS84 ellipsoid model of the Earth, plane projection provides 0.1% precision
11/// on distances under 500km at latitudes up to the 65°.
12/// See https://blog.mapbox.com/fast-geodesic-approximations-with-cheap-ruler-106f229ad016
13/// for more details about the principle and formulas behind.
14///
15/// ```
16/// use plane_projection::PlaneProjection;
17///
18/// let proj = PlaneProjection::new(55.65);
19/// let distance = proj.distance((55.704141722528554, 13.191304107330561), (55.60330902847681, 13.001973666557435));
20/// assert_eq!(distance as u32, 16373);
21/// ```
22pub struct PlaneProjection {
23    lon_scale: f64,
24    lat_scale: f64,
25}
26
27impl PlaneProjection {
28    /// Creates a plane projection to the Earth at provided latitude
29    pub fn new(latitude: f64) -> Self {
30        // Based on https://en.wikipedia.org/wiki/Earth_radius#Meridional
31        let cos_lat = latitude.to_radians().cos();
32        let w2 = 1.0 / (1.0 - SQUARED_ECCENTRICITY * (1.0 - cos_lat * cos_lat));
33        let w = w2.sqrt();
34
35        // multipliers for converting longitude and latitude degrees into distance
36        let lon_scale = (EQUATORIAL_RADIUS * w * cos_lat).to_radians(); // based on normal radius of curvature
37        let lat_scale = (EQUATORIAL_RADIUS * w * w2 * (1.0 - SQUARED_ECCENTRICITY)).to_radians(); // based on meridonal radius of curvature
38
39        Self {
40            lon_scale,
41            lat_scale,
42        }
43    }
44
45    /// Sqare distance in meters between two points in (lat, lon) format
46    #[inline(always)]
47    pub fn square_distance(&self, a: LatLon, b: LatLon) -> f64 {
48        let lat_dist = (a.0 - b.0) * self.lat_scale;
49        let lon_dist = lon_diff(a.1, b.1) * self.lon_scale;
50        lat_dist * lat_dist + lon_dist * lon_dist
51    }
52
53    /// Distance in meters between two points in (lat, lon) format
54    #[inline(always)]
55    pub fn distance(&self, a: LatLon, b: LatLon) -> f64 {
56        self.square_distance(a, b).sqrt()
57    }
58
59    /// Heading (azimuth) in degrees from point `a` to point `b` in a clockwise direction in range [0.0, 360.0)
60    /// where 0.0 is North, 90.0 is East, 180.0 is South and 270.0 is West.
61    #[inline(always)]
62    pub fn heading(&self, a: LatLon, b: LatLon) -> f32 {
63        // Using f32 calculations for better performance while maintaining sufficient precision
64        let dx = ((a.0 - b.0) * self.lat_scale) as f32;
65        let dy = (lon_diff(b.1, a.1) * self.lon_scale) as f32;
66        // Together with inverted `dx` this converts (-180, 180] `atan2` range into [0, 360) without branching
67        180.0 - dy.atan2(dx).to_degrees()
68    }
69}
70
71/// Returns the difference between two longitudes in range [-180.0, 180.0] degrees
72#[inline(always)]
73fn lon_diff(a: f64, b: f64) -> f64 {
74    let mut lon_diff = a - b;
75    if lon_diff > 180.0 {
76        lon_diff -= 360.0;
77    } else if lon_diff < -180.0 {
78        lon_diff += 360.0;
79    }
80    lon_diff
81}
82
83#[cfg(test)]
84mod tests {
85    use super::*;
86
87    #[test]
88    fn lon_diff_test() {
89        assert_eq!(lon_diff(0.0, 0.0), 0.0);
90        assert_eq!(lon_diff(100.0, 0.0), 100.0);
91        assert_eq!(lon_diff(100.0, -100.0), -160.0);
92        assert_eq!(lon_diff(177.0, -177.0), -6.0);
93        assert_eq!(lon_diff(358.0, 0.0), -2.0);
94        assert_eq!(lon_diff(0.0, 358.0), 2.0);
95        assert_eq!(lon_diff(0.0, -180.0), 180.0);
96        assert_eq!(lon_diff(1.0, -180.0), -179.0);
97        assert_eq!(lon_diff(180.0, 0.0), 180.0);
98        assert_eq!(lon_diff(180.0, -1.0), -179.0);
99        assert_eq!(lon_diff(180.0, -180.0), 0.0);
100    }
101
102    #[test]
103    fn distance_test() {
104        // From Lund C to Malmo C
105        let proj = PlaneProjection::new(55.65);
106        assert_eq!(
107            proj.distance(
108                (55.704141722528554, 13.191304107330561),
109                (55.60330902847681, 13.001973666557435),
110            ) as u32,
111            16373
112        );
113
114        let proj = PlaneProjection::new(51.05);
115        assert_eq!(
116            proj.distance((50.823194, 6.186389), (51.301389, 6.953333)) as u32,
117            75646
118        );
119    }
120
121    #[test]
122    fn heading_test() {
123        let proj = PlaneProjection::new(55.65);
124        assert_eq!(proj.heading((55.70, 13.19), (55.80, 13.19)) as i32, 0);
125        assert_eq!(proj.heading((55.70, 13.19), (55.60, 13.19)) as i32, 180);
126        assert_eq!(proj.heading((55.70, 13.19), (55.70, 13.29)) as i32, 90);
127        assert_eq!(proj.heading((55.70, 13.19), (55.70, 13.09)) as i32, 270);
128
129        // From Malmo C to Lund C
130        assert_eq!(
131            proj.heading(
132                (55.60330902847681, 13.001973666557435),
133                (55.704141722528554, 13.191304107330561),
134            ) as i32,
135            46
136        );
137        assert_eq!(
138            proj.heading(
139                (55.704141722528554, 13.191304107330561),
140                (55.60330902847681, 13.001973666557435),
141            ) as i32,
142            180 + 46
143        );
144    }
145}
plane_projection/lib.rs

plane_projection/
lib.rs
