plane_projection/
lib.rs

1// Values that define WGS84 ellipsoid model of the Earth in meters.
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///
22/// let heading = proj.heading((55.704141722528554, 13.191304107330561), (55.60330902847681, 13.001973666557435));
23/// assert_eq!(heading as u32, 226);
24/// ```
25#[derive(Clone)]
26pub struct PlaneProjection {
27    /// Meters per degree of longitude.
28    lon_scale: f64,
29    /// Meters per degree of latitude.
30    lat_scale: f64,
31}
32
33impl PlaneProjection {
34    /// Creates a plane projection to the Earth at provided latitude.
35    pub fn new(latitude: f64) -> Self {
36        // `cosf32` gives sufficient precision (adds approx. 0.0001 meter error) with much better performance
37        let cos_lat = (latitude as f32).to_radians().cos() as f64;
38
39        // Based on https://en.wikipedia.org/wiki/Earth_radius#Meridional
40        let w2 = 1.0 / (1.0 - SQUARED_ECCENTRICITY * (1.0 - cos_lat * cos_lat));
41        let w = w2.sqrt();
42        let lon_scale = (EQUATORIAL_RADIUS * w * cos_lat).to_radians(); // based on normal radius of curvature
43        let lat_scale = (EQUATORIAL_RADIUS * w * w2 * (1.0 - SQUARED_ECCENTRICITY)).to_radians(); // based on meridonal radius of curvature
44
45        Self {
46            lon_scale,
47            lat_scale,
48        }
49    }
50
51    /// Square distance in meters between two points in (lat, lon) format.
52    #[inline(always)]
53    pub fn square_distance(&self, a: LatLon, b: LatLon) -> f64 {
54        let lat_dist = (a.0 - b.0) * self.lat_scale;
55        let lon_dist = lon_diff(a.1, b.1) * self.lon_scale;
56        lat_dist * lat_dist + lon_dist * lon_dist
57    }
58
59    /// Distance in meters between two points in (lat, lon) format.
60    #[inline(always)]
61    pub fn distance(&self, a: LatLon, b: LatLon) -> f64 {
62        self.square_distance(a, b).sqrt()
63    }
64
65    /// Square distance in meters from point to the segment.
66    pub fn square_distance_to_segment(&self, point: LatLon, segment: (LatLon, LatLon)) -> f64 {
67        // Convert point and segment to projected space with origin at segment start
68        let mut point = (
69            (point.0 - segment.0.0) * self.lat_scale,
70            lon_diff(point.1, segment.0.1) * self.lon_scale,
71        );
72        let segment = (
73            (segment.1.0 - segment.0.0) * self.lat_scale,
74            lon_diff(segment.1.1, segment.0.1) * self.lon_scale,
75        );
76        if segment.0 != 0.0 || segment.1 != 0.0 {
77            let projection = (point.0 * segment.0 + point.1 * segment.1)
78                / (segment.0 * segment.0 + segment.1 * segment.1);
79            if projection > 1.0 {
80                point.0 -= segment.0;
81                point.1 -= segment.1;
82            } else if projection > 0.0 {
83                point.0 -= segment.0 * projection;
84                point.1 -= segment.1 * projection;
85            }
86        }
87        point.0 * point.0 + point.1 * point.1
88    }
89
90    /// Distance in meters from point to the segment.
91    #[inline(always)]
92    pub fn distance_to_segment(&self, point: LatLon, segment: (LatLon, LatLon)) -> f64 {
93        self.square_distance_to_segment(point, segment).sqrt()
94    }
95
96    /// Heading (azimuth) in degrees from point `a` to point `b` in the range [0.0, 360.0) degrees,
97    /// measured clockwise from North: 0.0 is North, 90.0 is East, 180.0 is South and 270.0 is West.
98    #[inline(always)]
99    pub fn heading(&self, a: LatLon, b: LatLon) -> f32 {
100        // Convert to f32 for better `atan2` performance while maintaining sufficient precision
101        let dx = ((a.0 - b.0) * self.lat_scale) as f32;
102        let dy = (lon_diff(b.1, a.1) * self.lon_scale) as f32;
103
104        // Together with inverted `dx` this converts (-180, 180] `atan2` range into [0, 360) without branching
105        180.0 - dy.atan2(dx).to_degrees()
106    }
107}
108
109/// Returns the difference between two longitudes in range [-180.0, 180.0] degrees.
110#[inline(always)]
111fn lon_diff(a: f64, b: f64) -> f64 {
112    let mut lon_diff = a - b;
113    if lon_diff > 180.0 {
114        lon_diff -= 360.0;
115    } else if lon_diff < -180.0 {
116        lon_diff += 360.0;
117    }
118    lon_diff
119}
120
121#[cfg(test)]
122mod tests {
123    use super::*;
124
125    const MALMO_C: LatLon = (55.60330902847681, 13.001973666557435);
126    const LUND_C: LatLon = (55.704141722528554, 13.191304107330561);
127    const STOCKHOLM_C: LatLon = (59.33036105663399, 18.058682977850953);
128
129    #[test]
130    fn lon_diff_test() {
131        assert_eq!(lon_diff(0.0, 0.0), 0.0);
132        assert_eq!(lon_diff(100.0, 0.0), 100.0);
133        assert_eq!(lon_diff(100.0, -100.0), -160.0);
134        assert_eq!(lon_diff(177.0, -177.0), -6.0);
135        assert_eq!(lon_diff(358.0, 0.0), -2.0);
136        assert_eq!(lon_diff(0.0, 358.0), 2.0);
137        assert_eq!(lon_diff(0.0, -180.0), 180.0);
138        assert_eq!(lon_diff(1.0, -180.0), -179.0);
139        assert_eq!(lon_diff(180.0, 0.0), 180.0);
140        assert_eq!(lon_diff(180.0, -1.0), -179.0);
141        assert_eq!(lon_diff(180.0, -180.0), 0.0);
142    }
143
144    #[test]
145    fn distance_test() {
146        let proj = PlaneProjection::new(55.65);
147        assert_eq!(proj.distance(MALMO_C, LUND_C).round() as u32, 16374);
148
149        // Geodesic distance is between Malmo and Stockholm is 513_861m and the best precision from
150        // the plane projection is when halfway latitude is used.
151        let proj = PlaneProjection::new((MALMO_C.0 + STOCKHOLM_C.0) * 0.5);
152        assert_eq!(proj.distance(MALMO_C, STOCKHOLM_C).round() as u32, 514_168); // 0.06% error
153        let proj = PlaneProjection::new(MALMO_C.0);
154        assert_eq!(proj.distance(MALMO_C, STOCKHOLM_C).round() as u32, 523_230); // 1.8% error
155        let proj = PlaneProjection::new(STOCKHOLM_C.0);
156        assert_eq!(proj.distance(MALMO_C, STOCKHOLM_C).round() as u32, 505_217); // 1.7% error
157    }
158
159    #[test]
160    fn distance_to_segment_test() {
161        let proj = PlaneProjection::new(0.0);
162        assert_eq!(
163            proj.distance_to_segment((0.0, 0.0), ((0.0, 0.0), (0.0, 1.0))),
164            0.0
165        );
166        assert_eq!(
167            proj.distance_to_segment((0.0, 0.5), ((0.0, 0.0), (0.0, 1.0))),
168            0.0
169        );
170        assert_eq!(
171            proj.distance_to_segment((0.0, 1.0), ((0.0, 0.0), (0.0, 1.0))),
172            0.0
173        );
174        assert_eq!(
175            proj.distance_to_segment((0.0, 0.0), ((0.0, 1.0), (0.0, 2.0))),
176            proj.distance((0.0, 0.0), (0.0, 1.0))
177        );
178        assert_eq!(
179            proj.distance_to_segment((0.0, 4.5), ((0.0, 1.0), (0.0, 3.0))),
180            proj.distance((0.0, 4.5), (0.0, 3.0))
181        );
182
183        // zero-length segment
184        assert_eq!(
185            proj.distance_to_segment((1.0, 2.0), ((3.0, 3.0), (3.0, 3.0))),
186            proj.distance((1.0, 2.0), (3.0, 3.0))
187        );
188
189        let proj = PlaneProjection::new(55.65);
190        assert_eq!(
191            proj.distance_to_segment(MALMO_C, (LUND_C, LUND_C)) as u32,
192            16373
193        );
194
195        // This is why it's tricky - point get projected not to a straight line in the lat/lon space,
196        // but to a segment in the projected space.
197        assert_eq!(
198            proj.distance_to_segment((55.67817981392954, 13.058789566271836), (MALMO_C, LUND_C))
199                as u32,
200            3615
201        );
202    }
203
204    #[test]
205    fn heading_test() {
206        let proj = PlaneProjection::new(55.65);
207        assert_eq!(proj.heading((55.70, 13.19), (55.80, 13.19)) as i32, 0);
208        assert_eq!(proj.heading((55.70, 13.19), (55.60, 13.19)) as i32, 180);
209        assert_eq!(proj.heading((55.70, 13.19), (55.70, 13.29)) as i32, 90);
210        assert_eq!(proj.heading((55.70, 13.19), (55.70, 13.09)) as i32, 270);
211
212        assert_eq!(proj.heading(MALMO_C, LUND_C,) as i32, 46);
213        assert_eq!(proj.heading(LUND_C, MALMO_C,) as i32, 180 + 46);
214    }
215}
plane_projection/lib.rs

plane_projection/
lib.rs
