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.

Answer 2

I'll elaborate on Lior Kogan's excellent answer by using a geometrical (instead of analytical) approach.

The great circle that contains the "line" lies on a plane passing through the centre of the spheroid. This plane is orthogonal to the vector obtained as cross-product of the vectors that pass through the origin and, respectively, p1 and p2.

Now, we are looking for the plane orthogonal to the one we have, passing for p0. This can be easily calculated, and the intersection of this plane with the spheroid should (WARNING: I'm in kinda of a hurry, and am not sure this step is mathematically sound) be the great circle orthogonal to the "line". The intersection of the arcs should be the point you are looking for, and is can be calculated as the interception of the line common to the two planes (cross product of the vectors orthogonal to each plane) and the spheroid.

Answer 3

thank you trying to teach me, but I have asked for a simple pragmatic solution. Not a math lecture.

Anyway, here is a simple answer:

Look here: http://www.movable-type.co.uk/scripts/latlong.html

Calculate "Bearing" (from start of the line to point distance from each I want to know) + "Cross-talk distance".

That is it. We are done in a couple of lines of code. No iterations. No additional libs. Lean and mean.

Answer 4

i havent worked with lon/lat before but.. i would create a coordinate system. center of earth being the origin. convert the lon/lat locations to the coordinate system do your calculations convert back to lon/lat

Answer 5

You have a line through earth-surface points A and B and a point C which you want to compute distance from.

You can compute the area of triangle ABC, divide it by distance between A and B and then multiply by 2.

function computeDistanceToLine(p, line) {
    var length = computeDistanceBetween(line.startPoint, line.endPoint);
    var triangleArea = computeTriangleArea(p, line.startPoint, line.endPoint);
    return 2 * triangleArea / length;
}

The algorithm for computing distance between two points is well-known. There are tons of implementations. One of them (as noticed in previous answers) can be found there http://www.movable-type.co.uk/scripts/latlong.html

For computing triangle area you can use some algorithm dependent on line lengths. For example

function computeTriangleArea(p0, p1, p2) {
    var r = 6378137;

    var d0 = computeDistanceBetween(p0, p1);
    var d1 = computeDistanceBetween(p1, p2);
    var d2 = computeDistanceBetween(p2, p0);

    var halfPerimeter = (d0 + d1 + d2) * 0.5;

    var t = Math.tan(halfPerimeter) *
        Math.tan((halfPerimeter - d0) * 0.5) *
        Math.tan((halfPerimeter - d1) * 0.5) *
        Math.tan((halfPerimeter - d2) * 0.5);

    return 4 * Math.atan(Math.sqrt(Math.abs(t))) * r * r;
}

Answer 6

Ilya Golota's answer is a correct approach but, turns out it always returns a really large number such as "3.0355243098445522e12" even for points that are 2 meters away from the line. I did this for computing the triangle area instead:

double d0 = computeDistanceBetween(p0, p1);
double d1 = computeDistanceBetween(p1, p2);
double d2 = computeDistanceBetween(p2, p0);

(Using Heron's formula)

double area = Math.sqrt(halfP*(halfP-d0)*(halfP-d1)*(halfP-d2));

And then plug in the area value here

double distanceToLine = 2*area/computeDistanceBetween(line.startPoint, line.endPoint)

This works well for me and returns the perpendicular distance from the point to the line in meters.