LatLon.nearestPointOnSegment() not working properly?

[eilamg]

Hi mrJean1,

I think I've found an issue with the nearestPointOnSegment method in sphericalNvector.LatLon.

>>> from PyGeodesy.geodesy.sphericalNvector import LatLon
>>> a = LatLon(0, 0)
>>> b = LatLon(0, 1)
>>> x = LatLon(1, 0)
>>> x.nearestPointOnSegment(a, b)
LatLon(00°00′00.0″N, 180°00′00.0″E)

Cheers, Eilam

[mrJean1]

Correct, there was an issue and that has been fixed. Also, there were no tests for that method and there are now, including your example. Thank you for reporting the issue.

[eilamg]

:+1:

[eilamg]

Hi Jean,

Further testing with new version shows the issue is only partially resolved. I've written this simple function to calculate the distance of a point from a segment.

def dist_p_ab(p, a, b):
    """Calculate distance (in radians) from point p to segment (a, b)"""

    n = p.nearestPointOnSegment(a, b)
    d = p.distanceTo(n, radius=1)

    return d

I've then tested this with a segment (a = LatLon(20, -40); b = LatLon(40, 60)) and a grid of latitudes and longitudes spanning the entire globe (lat = np.arange(-90, 90, 1); lon = np.arange(-180, 180, 1)), and rendered the results as an image using matplotlib.

The following (unrelated) code also illustrates the issue:

>>> from PyGeodesy.geodesy.sphericalNvector import LatLon
>>> LatLon(10, -140).nearestPointOnSegment(LatLon(0, 20), LatLon(0, 40))
LatLon(00°00′00.0″N, 140°00′00.0″W)

It seems like current implementation is failing for a subset of points.

Cheers, Eilam

[mrJean1]

This is suspicious and I'll need to look into this further. It may take a few weeks, though.

[mrJean1]

Courtesy of the original author, antipodal points (on the same hemisphere) are now handled in the sphericalNvector.LatLon method isWithinExtent.

[eilamg]

That was significantly shorter than a few weeks :smiley:

[chrisveness]

Nice visualisation!