geo_coding/distance/
earth.rs

1const WGS_84_A: f64 = 6_378_137.0;
2const WGS_84_B: f64 = 6_356_752.314_2;
3
4#[inline]
5fn to_normal_vector(location: &[f64; 2]) -> [f64; 3] {
6    // https://en.wikipedia.org/wiki/N-vector
7    let longitude = location[0];
8    let latitude = location[1];
9    let (sin_longitude, cos_longitude) = longitude.to_radians().sin_cos();
10    let (sin_latitude, cos_latitude) = latitude.to_radians().sin_cos();
11    [
12        cos_latitude * cos_longitude,
13        cos_latitude * sin_longitude,
14        sin_latitude,
15    ]
16}
17
18#[inline]
19fn dot(a: [f64; 3], b: [f64; 3]) -> f64 {
20    a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
21}
22
23#[inline]
24fn cross(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
25    [
26        a[1] * b[2] - a[2] * b[1], //
27        a[2] * b[0] - a[0] * b[2], //
28        a[0] * b[1] - a[1] * b[0], //
29    ]
30}
31
32#[inline]
33fn length(a: [f64; 3]) -> f64 {
34    dot(a, a).sqrt()
35}
36
37fn to_f64(location: &[i64; 2]) -> [f64; 2] {
38    let longitude = location[0] as f64 * 1e-9;
39    let latitude = location[1] as f64 * 1e-9;
40    [longitude, latitude]
41}
42
43/// Returns geographical distance between two points.
44///
45/// The first coordinate is longitude, the second coordinate is latitude.
46/// Both are in nanodegrees.
47///
48/// The distance is computed by converting each point to surface normals.
49///
50/// # References
51///
52/// - <https://en.wikipedia.org/wiki/N-vector>
53#[cfg_attr(docsrs, doc(cfg(feature = "std")))]
54pub fn earth_distance(a: &[i64; 2], b: &[i64; 2]) -> u64 {
55    earth_distance_f64(&to_f64(a), &to_f64(b)) as u64
56}
57
58pub(crate) fn earth_distance_f64(a: &[f64; 2], b: &[f64; 2]) -> f64 {
59    const R_AVG: f64 = (WGS_84_A + WGS_84_B) * 0.5;
60    let n1 = to_normal_vector(a);
61    let n2 = to_normal_vector(b);
62    R_AVG * length(cross(n1, n2)).atan2(dot(n1, n2))
63}
geo_coding/distance/earth.rs

geo_coding/distance/
earth.rs
