algorithm - Distance from point to line on Earth

Question

I need something as simple as

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

But that works with Lon/Lat.

Answer 1

We'll start with some assumptions:

• We'll model the figure of the earth as an oblate spheroid.
• By "line" we mean the minor arc of a great circle.
• We are looking for a highly accurate solution (for the suggested model)

So, rephrasing your question:

Given three earth-surface points - p0, p1 and p2, find an earth-surface point on the minor arc of the great circle defined by p1 and p2, which is the closest to p0.

As an infrastructure for the solution, We'll need:

• An accurate function that finds a target point based on an initial point, initial azimuth and distance.
• An accurate function that measures the distance between two points.

I suggest using GeographicLib's Direct and Inverse functions repectively, which are the most accurate implementations I know of.

Since that mathematics involved in oblate spheroid calculations is highly nonlinear, we'll build an iterative solution.

As a first step, we'll try to understand how a graph where the X axis is a point on the minor arc of the great circle defined by p1 and p2, and the Y axis is the distance from p0 to that point - may look like:

There are several options to how such graph may look like: The function may be monotonically increasing or monotonically decreasing. It may also contain a single point where its 1st derivative may be 0. It may be a minima (quite trivial), but it may as well be a maxima (for example - if Lat(p0)=0, Lat(p1)=100 and Lat(p2)=-100). However, in all cases there are 0 or 1 points where the derivative changes sign.

Understanding this, we can now build an iterative algorithm. In each iteration:

We'll calculate the dist(p0,p1), dist(p0,p2) and also dist(p0,pM) where M is the mid-point between p1 and p2 on the minor arc of a great circle definded by p1 and p2. Now. we'll check:

• if (dist(p0,p1) <= dist(p0,pM)) && (dist(p0,pM)<=dist(p0,p2)) - p0 is closer to p1 than it is to p2
• if (dist(p0,p2) <= dist(p0,pM)) && (dist(p0,pM)<=dist(p0,p1)) - p0 is closer to p2 than it is to p1
• if (dist(p0,p1) <= dist(p0,p2)) && (dist(p0,p2)<=dist(p0,pM)) - p0 is p1
• if (dist(p0,p2) <= dist(p0,p1)) && (dist(p0,p1)<=dist(p0,pM)) - p0 is p2

Otherwise, we can't determine if the minimum is closer to p1 or to p2, so we'll use two more points to check: We'll define pL as the mid-point between p1 and pM, and pN as the mid-point between pM and p2. Now,

• if (dist(p0,pL) <= dist(p0,pM)) - p0 is closer to p1
• if (dist(p0,pN) <= dist(p0,pM)) - p0 is closer to p2

Otherwise - p0 is between pL and pN.

So in each iteration we are halving the arc-length in which we are looking for a solution.

Using this method, we can get a 1 cm accuracy in less than 30 iterations.