Implement `HaversineDistance` for all combinations of geometry types

[frewsxcv]

Right now it's only implemented between two Points

https://docs.rs/geo/0.20.1/geo/algorithm/haversine_distance/trait.HaversineDistance.html

[michaelkirk]

Does anyone have any thoughts on how to go about this? Or have any examples of other implementations?

[urschrei]

I had a quick look earlier but couldn't find anything obvious. The first thing to note is that the closest point between a Point and a Line isn't necessarily a vertex – we'd need to have a look at the EuclideanDistance logic to see whether we currently determine the closest point in those cases, or whether we just calculate the distance directly (not at a computer rn, but will check).

[michaelkirk]

from http://www.faqs.org/faqs/graphics/algorithms-faq/

Subject 1.02: How do I find the distance from a point to a line?

Has a good euclidean answer (JTS uses this under the hood for geom_a.distance(geom_b)).

We have implemented this algorithm as geo_types::private_utils::point_line_euclidean_distance

I don't see how it could work in non-euclidean space though.

[urschrei]

Oh, I wonder whether we could use ClosestPoint to get the destination point…

[urschrei]

I think this is a correct way to work it out:

impl<T> HaversineDistance<T, Point<T>> for Polygon<T>
where
    T: GeoFloat + FromPrimitive,
{
    fn haversine_distance(&self, rhs: &Point<T>) -> T {
        match self.closest_point(rhs) {
            crate::Closest::Intersection(_) => T::zero(),
            crate::Closest::SinglePoint(sp) => rhs.haversine_distance(&sp),
            crate::Closest::Indeterminate => todo!(),
        }
    }
}

But we can't do this in practice because the HaversineDistance trait doesn't return a Result but we have to account for the possibility of the closest_point result being an intersection or indeterminate. Ugh.

[michaelkirk]

I think this only kicks the can down the road a bit, since closest_point is itself computed with assumptions on euclidean space.

We'd be getting the haversine meter distance to the point that is closest in terms of degrees, but may not be closest in terms of meters - right?

[michaelkirk]

I'm thinking we'd need something specially formulated for the spherish model of the earth. Maybe I finally need to read the chamberlain duquette paper?

https://trs.jpl.nasa.gov/bitstream/handle/2014/41271/07-0286.pdf

[michaelkirk]

I wasn't familiar with the term before, but "cross track distance" is the distance between a point and a great circle.

http://www.movable-type.co.uk/scripts/latlong.html#cross-track https://gis.stackexchange.com/questions/209540/projecting-cross-track-distance-on-great-circle

geographiclib (karney geodesics): https://sourceforge.net/p/geographiclib/discussion/1026621/thread/299518a3e4/

Maybe that's useful.

[urschrei]

Yeah, that's a good point. It may be that there is no straightforward way to perform this calculation with our existing toolbox, because in order for it to be correct (and even then, only sort-of correct because the earth isn't a sphere) you have to implement all the necessary functions using spherical trigonometry, and I for one do not have the chops.

Taking a step back, there is a reason that geodetic calculations generally punt on this and try to project geometries into planar coordinates whenever possible: it avoids a lot of complexity. It's easy to say "we should be able to calculate the closest point between a point and a geometry on a sphere", but the more I think about it, the less confident I am that that's true; we now have to contend with primitives such as curves and ensuring that planar and non-planar geometries aren't mixed.

I am not saying "we shouldn't try to do this", just…perhaps we need to think about it a bit more carefully.

[urschrei]

That GIS StackExchange answer about projecting points onto the great circle is a lot. I think that's probably Jason Davies.

[urschrei]

*Looks at https://link.springer.com/content/pdf/10.1007/s00190-012-0578-z.pdf and immediately passes out

[michaelkirk]

One idea for point_line distance that relies on proj, but otherwise might not introduce too much new machinery:

fn haversine_distance(point: Point, line: Line) -> meters {
    // https://en.wikipedia.org/wiki/Azimuthal_equidistant_projection
    // https://proj.org/operations/projections/aeqd.html
    // ... has the useful properties that all points on the map are at proportionally 
    // correct distances from the center point,
    let proj = Proj::new(some_azimuthal_equidistant_projection, center = point);

    let point_projected  = point.transformed(proj); // == (0, 0) since we defined the point as the center of the projection
    let line_projected = line.transformed(proj);

    let closest_point_on_line = point_projected.closest_point(line_projected);
    let unprojected_closest_point_on_line = closest_point_on_line.inverse_transform(proj);

    point.haversine_distance(unprojected_closest_point_on_line)
}

Edit: We'd probably want some kind of intersection check to short circuit the above and return 0. IIRC chamberlain/duquette have a formula for "is point in polygon" in https://trs.jpl.nasa.gov/bitstream/handle/2014/41271/07-0286.pdf

That doesn't solve line/line distance though.

[urschrei]

Ooooh. That's still a lot of progress, though.

[michaelkirk]

*Looks at https://link.springer.com/content/pdf/10.1007/s00190-012-0578-z.pdf and immediately passes out

You know it's going to be a good paper when a bunch of the citations are from the 19th century (and one from the 18th!).

[michaelkirk]

Some maybe related discussion on point-to-polygon distances https://github.com/locationtech/jts/issues/929

[mbcltd]

I have had a crack at implementing cross track distance based on this http://www.movable-type.co.uk/scripts/latlong.html?from=48.86,-122.0992&to=48.8599,-122.1449.

Hopefully it's in keeping with the style of HaversineDistance.

use num_traits::FromPrimitive;
use geo_types::{CoordFloat, Line, Point};
use crate::{Bearing, HaversineDistance, MEAN_EARTH_RADIUS};

pub trait CrossTrackDistance<T, Rhs> {
    fn cross_track_distance(&self, rhs: &Rhs) -> T;
}

impl<T> CrossTrackDistance<T, Line<T>> for Point<T>
    where
        T: CoordFloat + FromPrimitive,
{
    fn cross_track_distance(&self, rhs: &Line<T>) -> T {
        cross_track_distance(self, rhs)
    }
}

impl<T> CrossTrackDistance<T, Point<T>> for Line<T>
    where
        T: CoordFloat + FromPrimitive,
{
    fn cross_track_distance(&self, rhs: &Point<T>) -> T {
        cross_track_distance(rhs, self)
    }
}

fn cross_track_distance<T: CoordFloat + FromPrimitive>(point: &Point<T>, line: &Line<T>) -> T {
    let mean_earth_radius = T::from(MEAN_EARTH_RADIUS).unwrap();
    let path_start = Point(line.start);
    let path_end = Point(line.end);
    let l_delta_13: T = path_start.haversine_distance(point) / mean_earth_radius;
    let theta_13: T = path_start.bearing(*point).to_radians();
    let theta_12: T = path_start.bearing(path_end).to_radians();
    let l_delta_xt: T = (l_delta_13.sin() * (theta_12 - theta_13).sin()).asin();
    mean_earth_radius * l_delta_xt
}

I have two reasonable tests which both pass:

#[cfg(test)]
mod test {
    use geo_types::Line;
    use crate::CrossTrackDistance;
    use crate::HaversineDistance;
    use crate::Point;

    #[test]
    fn should_equal_haversine_distance_for_short_lines() {
        let a = Point::new(0., 0.);
        let b = Line::new([1., 0.], [1., 0.1]);
        let c = Point::new(1., 0.);
        assert_relative_eq!(
            a.cross_track_distance(&b),
            a.haversine_distance(&c),
            epsilon = 1.0e-6
        );
    }

    #[test]
    fn distance1_test() {
        let a = Point::new(-0.7972, 53.2611);
        let b = Line::new([-1.7297, 53.3206], [0.1334, 53.1887]);
        assert_relative_eq!(
            a.cross_track_distance(&b),
            307.549995,
            epsilon = 1.0e-6
        );
    }
}

Happy to do a PR and fix up the docs if this looks like a useful stepping stone to the ultimate goal of a calculating a minimum distance between any set of Points and any set of Points or Lines.

[mbcltd]

I've just done another PR for the distance between a point and a bounded line based on the cross track distance. Next step I think is to implement a something for every geometry to every other depending on whether the two geometries should be decomposed into points or lines. Brute force should perform in O(n^2), more efficient algorithms could then come later.

[frewsxcv]

An attempt was made in https://github.com/georust/geo/pull/1015 if anyone wants to pick up where this one left off