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 = remainder(a.1 - b.1, 360.0) * 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
60/// Returns the IEEE 754-2008 floating-point remainder of the division x / y,
61/// thus normalizing the x in range [-y/2.0, y/2.0]
62#[inline(always)]
63fn remainder(x: f64, y: f64) -> f64 {
64    x - y * (x / y).round()
65}
66
67#[cfg(test)]
68mod tests {
69    use super::*;
70
71    #[test]
72    fn reminder_test() {
73        assert_eq!(remainder(0.0, 360.0), 0.0);
74        assert_eq!(remainder(358.0, 360.0), -2.0);
75        assert_eq!(remainder(-358.0, 360.0), 2.0);
76        assert_eq!(remainder(-180.0, 360.0), 180.0);
77        assert_eq!(remainder(180.0, 360.0), -180.0);
78    }
79
80    #[test]
81    fn distance_test() {
82        // From Lund C to Malmo C
83        let proj = PlaneProjection::new(55.65);
84        assert_eq!(
85            proj.distance(
86                (55.704141722528554, 13.191304107330561),
87                (55.60330902847681, 13.001973666557435)
88            ) as u32,
89            16373
90        );
91
92        let proj = PlaneProjection::new(51.05);
93        assert_eq!(
94            proj.distance((50.823194, 6.186389), (51.301389, 6.953333)) as u32,
95            75646
96        );
97    }
98}
plane_projection/lib.rs

plane_projection/
lib.rs
