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

plane_projection/
lib.rs
