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    /// Projects a coordinate from (latitude, longitude) to the plane projection space.
52    ///
53    /// This function is intended for low-level coordinate manipulation (like vector math) in the projection space
54    /// and should not be used unless the built-in methods like [`PlaneProjection::distance()`] and
55    /// [`PlaneProjection::distance_to_segment()`] are insufficient for your use case.
56    #[inline(always)]
57    pub fn project(&self, ll: LatLon) -> (f64, f64) {
58        (ll.0 * self.lat_scale, ll.1 * self.lon_scale)
59    }
60
61    /// Square distance in meters between two points in (lat, lon) format.
62    #[inline(always)]
63    pub fn square_distance(&self, a: LatLon, b: LatLon) -> f64 {
64        let lat_dist = (a.0 - b.0) * self.lat_scale;
65        let lon_dist = lon_diff(a.1, b.1) * self.lon_scale;
66        lat_dist * lat_dist + lon_dist * lon_dist
67    }
68
69    /// Distance in meters between two points in (lat, lon) format.
70    #[inline(always)]
71    pub fn distance(&self, a: LatLon, b: LatLon) -> f64 {
72        self.square_distance(a, b).sqrt()
73    }
74
75    /// Square distance in meters from point to the segment.
76    pub fn square_distance_to_segment(&self, point: LatLon, segment: (LatLon, LatLon)) -> f64 {
77        // Convert point and segment to projected space with origin at segment start
78        let mut point = (
79            (point.0 - segment.0.0) * self.lat_scale,
80            lon_diff(point.1, segment.0.1) * self.lon_scale,
81        );
82        let segment = (
83            (segment.1.0 - segment.0.0) * self.lat_scale,
84            lon_diff(segment.1.1, segment.0.1) * self.lon_scale,
85        );
86        if segment.0 != 0.0 || segment.1 != 0.0 {
87            let projection = (point.0 * segment.0 + point.1 * segment.1)
88                / (segment.0 * segment.0 + segment.1 * segment.1);
89            if projection > 1.0 {
90                point.0 -= segment.0;
91                point.1 -= segment.1;
92            } else if projection > 0.0 {
93                point.0 -= segment.0 * projection;
94                point.1 -= segment.1 * projection;
95            }
96        }
97        point.0 * point.0 + point.1 * point.1
98    }
99
100    /// Distance in meters from point to the segment.
101    #[inline(always)]
102    pub fn distance_to_segment(&self, point: LatLon, segment: (LatLon, LatLon)) -> f64 {
103        self.square_distance_to_segment(point, segment).sqrt()
104    }
105
106    /// Heading (azimuth) in degrees from point `a` to point `b` in the range [0.0, 360.0) degrees,
107    /// measured clockwise from North: 0.0 is North, 90.0 is East, 180.0 is South and 270.0 is West.
108    #[inline(always)]
109    pub fn heading(&self, a: LatLon, b: LatLon) -> f32 {
110        // Convert to f32 for better `atan2` performance while maintaining sufficient precision
111        let dx = ((a.0 - b.0) * self.lat_scale) as f32;
112        let dy = (lon_diff(b.1, a.1) * self.lon_scale) as f32;
113
114        // Together with inverted `dx` this converts (-180, 180] `atan2` range into [0, 360) without branching
115        180.0 - dy.atan2(dx).to_degrees()
116    }
117}
118
119/// Returns the difference between two longitudes in range [-180.0, 180.0] degrees.
120#[inline(always)]
121fn lon_diff(a: f64, b: f64) -> f64 {
122    let mut lon_diff = a - b;
123    if lon_diff > 180.0 {
124        lon_diff -= 360.0;
125    } else if lon_diff < -180.0 {
126        lon_diff += 360.0;
127    }
128    lon_diff
129}
130
131#[cfg(test)]
132mod tests {
133    use super::*;
134
135    const MALMO_C: LatLon = (55.60330902847681, 13.001973666557435);
136    const LUND_C: LatLon = (55.704141722528554, 13.191304107330561);
137    const STOCKHOLM_C: LatLon = (59.33036105663399, 18.058682977850953);
138
139    #[test]
140    fn lon_diff_test() {
141        assert_eq!(lon_diff(0.0, 0.0), 0.0);
142        assert_eq!(lon_diff(100.0, 0.0), 100.0);
143        assert_eq!(lon_diff(100.0, -100.0), -160.0);
144        assert_eq!(lon_diff(177.0, -177.0), -6.0);
145        assert_eq!(lon_diff(358.0, 0.0), -2.0);
146        assert_eq!(lon_diff(0.0, 358.0), 2.0);
147        assert_eq!(lon_diff(0.0, -180.0), 180.0);
148        assert_eq!(lon_diff(1.0, -180.0), -179.0);
149        assert_eq!(lon_diff(180.0, 0.0), 180.0);
150        assert_eq!(lon_diff(180.0, -1.0), -179.0);
151        assert_eq!(lon_diff(180.0, -180.0), 0.0);
152    }
153
154    #[test]
155    fn distance_test() {
156        let proj = PlaneProjection::new(55.65);
157        assert_eq!(proj.distance(MALMO_C, LUND_C).round() as u32, 16374);
158
159        // Geodesic distance is between Malmo and Stockholm is 513_861m and the best precision from
160        // the plane projection is when halfway latitude is used.
161        let proj = PlaneProjection::new((MALMO_C.0 + STOCKHOLM_C.0) * 0.5);
162        assert_eq!(proj.distance(MALMO_C, STOCKHOLM_C).round() as u32, 514_168); // 0.06% error
163        let proj = PlaneProjection::new(MALMO_C.0);
164        assert_eq!(proj.distance(MALMO_C, STOCKHOLM_C).round() as u32, 523_230); // 1.8% error
165        let proj = PlaneProjection::new(STOCKHOLM_C.0);
166        assert_eq!(proj.distance(MALMO_C, STOCKHOLM_C).round() as u32, 505_217); // 1.7% error
167    }
168
169    #[test]
170    fn distance_to_segment_test() {
171        let proj = PlaneProjection::new(0.0);
172        assert_eq!(
173            proj.distance_to_segment((0.0, 0.0), ((0.0, 0.0), (0.0, 1.0))),
174            0.0
175        );
176        assert_eq!(
177            proj.distance_to_segment((0.0, 0.5), ((0.0, 0.0), (0.0, 1.0))),
178            0.0
179        );
180        assert_eq!(
181            proj.distance_to_segment((0.0, 1.0), ((0.0, 0.0), (0.0, 1.0))),
182            0.0
183        );
184        assert_eq!(
185            proj.distance_to_segment((0.0, 0.0), ((0.0, 1.0), (0.0, 2.0))),
186            proj.distance((0.0, 0.0), (0.0, 1.0))
187        );
188        assert_eq!(
189            proj.distance_to_segment((0.0, 4.5), ((0.0, 1.0), (0.0, 3.0))),
190            proj.distance((0.0, 4.5), (0.0, 3.0))
191        );
192
193        // zero-length segment
194        assert_eq!(
195            proj.distance_to_segment((1.0, 2.0), ((3.0, 3.0), (3.0, 3.0))),
196            proj.distance((1.0, 2.0), (3.0, 3.0))
197        );
198
199        let proj = PlaneProjection::new(55.65);
200        assert_eq!(
201            proj.distance_to_segment(MALMO_C, (LUND_C, LUND_C)) as u32,
202            16373
203        );
204
205        // This is why it's tricky - point get projected not to a straight line in the lat/lon space,
206        // but to a segment in the projected space.
207        assert_eq!(
208            proj.distance_to_segment((55.67817981392954, 13.058789566271836), (MALMO_C, LUND_C))
209                as u32,
210            3615
211        );
212    }
213
214    #[test]
215    fn heading_test() {
216        let proj = PlaneProjection::new(55.65);
217        assert_eq!(proj.heading((55.70, 13.19), (55.80, 13.19)) as i32, 0);
218        assert_eq!(proj.heading((55.70, 13.19), (55.60, 13.19)) as i32, 180);
219        assert_eq!(proj.heading((55.70, 13.19), (55.70, 13.29)) as i32, 90);
220        assert_eq!(proj.heading((55.70, 13.19), (55.70, 13.09)) as i32, 270);
221
222        assert_eq!(proj.heading(MALMO_C, LUND_C,) as i32, 46);
223        assert_eq!(proj.heading(LUND_C, MALMO_C,) as i32, 180 + 46);
224    }
225}
plane_projection/lib.rs

plane_projection/
lib.rs
