Distance from Lat/Lng point to Minor Arc segment

Question

I need to calculate the shortest distance from a lat/lng GPS point P to a line segment described by 2 other lat/lng GPS points A and B.

'Cross-track distance' helps me to calculate the shortest distance between P and the great circle described by A and B.

However, this is not what I want. I need need the distance between P and the line segment of A-B, not the entire great circle.

I have used the following implementation from http://www.movable-type.co.uk/scripts/latlong.html

Formula:    dxt = asin( sin(δ13) ⋅ sin(θ13−θ12) ) ⋅ R
where:
δ13 is (angular) distance from start point to third point
θ13 is (initial) bearing from start point to third point
θ12 is (initial) bearing from start point to end point
R is the earth’s radius

The following images hopefully demonstrate the problem I am trying to solve:

In the first image the Cross-Track distance, indicated by the green line is correct and indeed the shortest distance to the line segment AB.

In the second image the problem with cross-track distance is shown, In this case I would want the shortest distance to be the simple distance AP, but Cross-Track distance gives me the distance indicated by the red line.

How do I change my algoritm to take this into account, or check whether or not point X is within AB. Is it possible to do this computationally? Or is iterative the only possible (expensive) solution? (take N points along AB and calculate the min distance from P to all these points)

For simplicity purposes all lines in the images are straight. In reality, these are minor arcs on a great circle

What are the maximum distances in that application? some hundred meters, or up to some thousand kilomters? — AlexWien, Sep 25 '15 at 16:09
The line segments are waypoints in a (driven) route. I would say the distances between them would be in the 100-1000m range, roughly. — ChrisDekker, Sep 25 '15 at 16:30
Did you already translate the cross track distance to matlab? If so could you put the code into your question? — Daniel, Sep 26 '15 at 16:23
In the plane, you'd just check if all angles in the ABP triangle were <= 90 degrees. Is there something similar that could be applied from spherical trig? — Chad Kennedy, Sep 28 '15 at 14:52
Sure Chad, was thinking the same thing at first, but I could come up with several scenario's where that would not produce accurate results, including points spanning different hemispheres and instances where the minor arc segment was very small — ChrisDekker, Sep 30 '15 at 02:57
You ever find a solution? I'm in the same boat as you. I used the algorithms from that site but I can't figure out how to pick the cross-track or closest end-point. If I could get the actual coordinates of the cross-track point I could figure it out. If the cross-track point is the right one (like in your first figure), the X coordinates would be between the two segment coordinates. So X-lat would be be between A-lat and B-lat, same with the longitudes. However I don't think that site has an example of that and I don't understand the lat-long math well enough to figure it out myself. — Zip184, Feb 18 '16 at 13:55

score 30 · Accepted Answer · edited Apr 30 '20 at 08:08

First, some nomenclature:
Our arc is drawn from p1 to p2.
Our third point is p3.
The imaginary point that intersects the great circle is p4.
p1 is defined by lat1,lon1; p2 by lat2,lon2; etc.
dis12 is the distance from p1 to p2; etc.
bear12 is the bearing from p1 to p2; etc.
dxt is cross-track distance.
dxa is cross-arc distance, our goal!

Notice that the cross-track formula relies on the relative bearing, bear13-bear12

We have 3 cases to deal with.

Case 1: The relative bearing is obtuse. So, dxa=dis13.

Case 2.1: The relative bearing is acute, AND p4 falls on our arc. So, dxa=dxt.

Case 2.2: The relative bearing is acute,AND p4 falls beyond our arc. So, dxa=dis23

The algorithm:

Step 1: If relative bearing is obtuse, dxa=dis13
Done!
Step 2: If relative bearing is acute:
2.1: Find dxt.
2.3: Find dis12.
2.4: Find dis14.
2.4: If dis14>dis12, dxa=dis23.
Done!
2.5: If we reach here, dxa=abs(dxt)

MATLAB code:

function [ dxa ] = crossarc( lat1,lon1,lat2,lon2,lat3,lon3 )
%// CROSSARC Calculates the shortest distance in meters 
%// between an arc (defined by p1 and p2) and a third point, p3.
%// Input lat1,lon1,lat2,lon2,lat3,lon3 in degrees.
    lat1=deg2rad(lat1); lat2=deg2rad(lat2); lat3=deg2rad(lat3);
    lon1=deg2rad(lon1); lon2=deg2rad(lon2); lon3=deg2rad(lon3);

    R=6371000; %// Earth's radius in meters
    %// Prerequisites for the formulas
    bear12 = bear(lat1,lon1,lat2,lon2);
    bear13 = bear(lat1,lon1,lat3,lon3);
    dis13 = dis(lat1,lon1,lat3,lon3);

    diff = abs(bear13-bear12);
    if diff > pi
        diff = 2 * pi - diff;
    end
    %// Is relative bearing obtuse?
    if diff>(pi/2)
        dxa=dis13;
    else
        %// Find the cross-track distance.
        dxt = asin( sin(dis13/R)* sin(bear13 - bear12) ) * R;

        %// Is p4 beyond the arc?
        dis12 = dis(lat1,lon1,lat2,lon2);
        dis14 = acos( cos(dis13/R) / cos(dxt/R) ) * R;
        if dis14>dis12
            dxa=dis(lat2,lon2,lat3,lon3);
        else
            dxa=abs(dxt);
        end   
    end
end

function [ d ] = dis( latA, lonA, latB, lonB )
%DIS Finds the distance between two lat/lon points.
R=6371000;
d = acos( sin(latA)*sin(latB) + cos(latA)*cos(latB)*cos(lonB-lonA) ) * R;
end

function [ b ] = bear( latA,lonA,latB,lonB )
%BEAR Finds the bearing from one lat/lon point to another.
b=atan2( sin(lonB-lonA)*cos(latB) , ...
    cos(latA)*sin(latB) - sin(latA)*cos(latB)*cos(lonB-lonA) );
end

Sample outputs: Demonstrate all cases. See maps below.

>> crossarc(-10.1,-55.5,-15.2,-45.1,-10.5,-62.5)
ans =
   7.6709e+05
>> crossarc(40.5,60.5,50.5,80.5,51,69)
ans =
   4.7961e+05
>> crossarc(21.72,35.61,23.65,40.7,25,42)
ans =
   1.9971e+05

Those same outputs on the map!:

Demonstrates case 1:

Demonstrates case 2.1:

Demonstrates case 2.2:

Credit to: http://www.movable-type.co.uk/scripts/latlong.html
for the formulas
and: http://www.darrinward.com/lat-long/?id=1788764
for generating the map images.

Great explanation, thanks! Could you add how I can find the lat/lon coords of the point along the segment? I could use the complex formula "Destination point given distance and bearing from start point" from Movable-type.co.uk, but I think I will be doing some expensive calculations double then. — ChrisDekker, Apr 29 '16 at 04:54
Thank you! If I am not mistaken current code has slight bug, this condition is too weak `if abs(bear13-bear12)>(pi/2)`, because one can have mirror-like cases. Thus first it is best to compute difference -- `diff = abs(bear13-bear12)`. And then check if the difference is greater than 180 degrees -- `if diff>PI then diff = 2*PI-diff` (this is not normalization, but mirror-flip). The rest follows `if diff>pi/2...` — greenoldman, Jan 13 '19 at 14:19
Great answer, thank you for such a detailed explanation! However, @greenoldman's comment is correct and this answer has to be edited. Here's an example to highlight the problem: `crossarc(-29.7762, 30.98883,-29.77483, 30.98877, -29.775664, 30.988693)` and `crossarc(-29.77483, 30.98877,-29.7762, 30.98883, -29.775664, 30.988693)` compute the distance from the same point, to the same segment (but with p1 and p2 switched) and they return different values. — Cristian Lupascu, Apr 29 '20 at 13:20
Will the "If dis14>dis12" be enough for the third case though? What if we have a right angle at the second point? — VuVirt, Mar 15 '23 at 08:26

score 1 · Answer 2 · edited May 23 '17 at 12:02

1

For 100 - 1000m spherical problems, it is easy to just convert to cartesian space, using a equirectangular projection.
Then it continues with school mathematics:
Use the function "distance from line segment" which is easy to find ready implemented. This fucntion uses (and sometimes returns) a relative forward/backward position for the projected point X on the line A,B. The value is

in the interval [0,1] if the projected point is inside the line segment.
it is negative if X is outside before A,
it is >1 if outside after B.
If the relative position is between 0,1 the normal distance is taken, if outside the shorter distance of the both start and line-end points, A,B.

An example of such / or very similar an cartesian implementaion is Shortest distance between a point and a line segment

edited May 23 '17 at 12:02

Community

1
1

answered Sep 28 '15 at 14:30

AlexWien

28,470
6
53
83

Any idea what the accuracy drop-off is when distances increase? Considering we are treating distances on a sphere as if they were in a 2D plane. – ChrisDekker Sep 30 '15 at 02:55
Yes it is based by the simple equirectangular projection, which treats the earth as cyclinder, instead of a sphere. For small distancres a plane is the same as the surface of a sphere (as long as the coordinate system is transformed, to have equals units on both axes, which is not the case for the longitude. There are better projections, like UTM,which is used inhihking maps, etc. but they are more complex. – AlexWien Sep 30 '15 at 14:40

Distance from Lat/Lng point to Minor Arc segment

2 Answers2

Linked