Project location on a (great circle) path

Question

I've been searching this SE site for quite a few hours now, and I'm still struggling to find a solution to my question. My goal is that given a way in OSM and my location (lat/lon coordinates), I want to find the closest location (lat/lon coordinates) on that way. The point can be anywhere on the way, not limited to the points used to define the way.

So I'm thinking of the following algorithm:

1. Separate Path into separate edges, each edge connecting only two points.
2. Select the closest edge.
3. Project my location onto that edge.

Now there are many questions about calculating the distance between a location and a path:

• WGS point to WGS line segment (great circle) distance
• Calculating the distance between a point and a virtual line of two lat/lngs
• How to approximate point to segment distance on sphere?

Also a very similar question of which I cannot get the calculations right or verified:

• Find the nearest point on a line

There is also some info from Dr. Math about that subject. However I cannot seem to find an algorithm to calculate the location in step 3. As I haven't touched (vector) algebra in quite a while, I don't quite understand the logic in those answers.

Can someone show an algorithm to do this? A solution in any reasonable programming language is fine with me.

Answer 1

Using a spherical model of the earth may give adequate accuracy and leads to simple fast calculations.

Convert all coordinates into earth-centered (3D) cartesian coordinates. For example, the formula

(cos(lon)*cos(lat), sin(lon)*cos(lat), sin(lat))

will do. (It uses a distance measure in which the earth's radius is one unit, which is convenient.)

Writing X0 = (x0,y0,z0) for the start point and X1 = (x1,y1,z1) for the destination point, which define the great circle (provided X0 is distinct from X1 and the two are not diametrically opposite), let U be the normalized cross product of X0 and X1. This is computed in two steps:

V = (xv, yv, zv) = (y0*z1 - z0*y1, z0*x1 - x0*z1, x0*y1 - y0*x1)

The length of V is

|V| = sqrt(xv^2 + yv^2 + zv^2)

Normalization stretches V to unit length:

U = (xu, yu, zu) = V / |V| = (xv/|V|, yv/|V|, zv/|V|).

The oriented 3D distance between any point X = (x,y,z) and the plane of this great circle is just the dot product of X with Z, given by

d = X * U = x*xu + y*yu + z*zu

The closest point in terms of the distance on the earth's surface is the one that is closest to the plane: thus, it has the smallest absolute value of d.

This figure shows a great circle (in black) determined by the two white points and 2000 random points on the sphere colored and shaded according to their absolute 3D distance to the plane of that great circle; that is, |d|.

Having found a closest point, project it to the great circle by first projecting it to the great circle's plane (in 3D) and then extending that radially outward to the earth's surface. The projection merely subtracts d * U:

X' = (x', y', z') = X - d*U = (x - d*xu, y - d*yu, z - d*zu).

The radial projection simply renormalizes X' in the same way V was renormalized to U:

X'' = X' / |X'|.

(This will be problematic if |X'| = 0, which happens when the closest point is one of the poles of the great circle. Include a test in the code for this condition, if it could happen, and deal with it separately, using the sign of d to identify which pole.)

If desired, convert the coordinates of X'' back to (lat, lon) using the usual formulas.