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.