maps_engine_rust/
geo.rs

1//! Geographic math — distances, bearings, destinations on the sphere.
2
3use std::f64::consts::PI;
4
5/// Mean Earth radius in meters.
6pub const EARTH_R: f64 = 6_371_008.8;
7
8#[inline]
9fn rad(d: f64) -> f64 {
10    d * PI / 180.0
11}
12#[inline]
13fn deg(r: f64) -> f64 {
14    r * 180.0 / PI
15}
16
17/// A geographic coordinate.
18#[derive(Debug, Clone, Copy, PartialEq)]
19pub struct LonLat {
20    pub lon: f64,
21    pub lat: f64,
22}
23
24impl LonLat {
25    pub fn new(lon: f64, lat: f64) -> Self {
26        LonLat { lon, lat }
27    }
28}
29
30/// Great-circle distance in meters (haversine).
31pub fn haversine(a: LonLat, b: LonLat) -> f64 {
32    let (la1, la2, dlo, dla) = (
33        rad(a.lat),
34        rad(b.lat),
35        rad(b.lon - a.lon),
36        rad(b.lat - a.lat),
37    );
38    let h = (dla / 2.0).sin().powi(2) + la1.cos() * la2.cos() * (dlo / 2.0).sin().powi(2);
39    2.0 * EARTH_R * h.sqrt().asin()
40}
41
42/// Initial bearing from `a` to `b` in degrees clockwise from north.
43pub fn bearing(a: LonLat, b: LonLat) -> f64 {
44    let (la1, la2, dlo) = (rad(a.lat), rad(b.lat), rad(b.lon - a.lon));
45    let y = dlo.sin() * la2.cos();
46    let x = la1.cos() * la2.sin() - la1.sin() * la2.cos() * dlo.cos();
47    (deg(y.atan2(x)) + 360.0) % 360.0
48}
49
50/// Destination point from `a`, travelling `dist_m` meters on `bearing_deg`.
51pub fn destination(a: LonLat, dist_m: f64, bearing_deg: f64) -> LonLat {
52    let (la1, lo1, br, d) = (rad(a.lat), rad(a.lon), rad(bearing_deg), dist_m / EARTH_R);
53    let la2 = (la1.sin() * d.cos() + la1.cos() * d.sin() * br.cos()).asin();
54    let lo2 = lo1 + (br.sin() * d.sin() * la1.cos()).atan2(d.cos() - la1.sin() * la2.sin());
55    LonLat::new(deg(lo2), deg(la2))
56}
57
58/// Total length of a polyline in meters.
59pub fn polyline_length(points: &[LonLat]) -> f64 {
60    points.windows(2).map(|w| haversine(w[0], w[1])).sum()
61}
62
63#[cfg(test)]
64mod tests {
65    use super::*;
66
67    #[test]
68    fn hanoi_to_haiphong_distance() {
69        let hanoi = LonLat::new(105.85, 21.02);
70        let haiphong = LonLat::new(106.68, 20.86);
71        let d = haversine(hanoi, haiphong);
72        assert!((d - 88_000.0).abs() < 5_000.0, "got {d}");
73    }
74
75    #[test]
76    fn bearing_due_north() {
77        let a = LonLat::new(105.0, 21.0);
78        let b = LonLat::new(105.0, 22.0);
79        assert!((bearing(a, b) - 0.0).abs() < 0.5);
80    }
81
82    #[test]
83    fn destination_roundtrip() {
84        let a = LonLat::new(105.85, 21.02);
85        let b = destination(a, 10_000.0, 90.0);
86        let d = haversine(a, b);
87        assert!((d - 10_000.0).abs() < 1.0, "got {d}");
88        assert!((bearing(a, b) - 90.0).abs() < 0.5);
89    }
90
91    #[test]
92    fn polyline_sums_segments() {
93        let pts = [
94            LonLat::new(0.0, 0.0),
95            LonLat::new(0.0, 1.0),
96            LonLat::new(0.0, 2.0),
97        ];
98        let total = polyline_length(&pts);
99        let one = haversine(pts[0], pts[1]);
100        assert!((total - 2.0 * one).abs() < 1.0);
101    }
102}
maps_engine_rust/geo.rs

maps_engine_rust/
geo.rs
