# Project location on a (great circle) path

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:

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

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.

• Since it seems critical to your "rejection" of the other questions, please elaborate on "project my location onto that edge". The projection may not be on the edge. I believe that issue is addressed in the other questions. (Well done, for the research, BTW.) – Martin F Apr 25 '15 at 17:16
• @MartinF that question calculates the distance from a point to a line, but not the closest point on the line itself. – bouke Apr 25 '15 at 18:36
• There is a solution at gis.stackexchange.com/a/23500/3195 though it is maybe hard to comprehend. – Martin F Apr 25 '15 at 20:40
• Ah yes thanks, I've updated ref no. 3. The 'solution' in that particular question links to a general explanation of the problem field. While this might be sufficient for well-grounded mathematicians, I don't quite understand the mathematics in that paper. – bouke Apr 25 '15 at 22:40

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.

• One question. Consider the non too unusual case where we can choose any X1 and X0 (on the great circle), from an accuracy point of view, is it better to pick X1 and X0 close or far apart (again provided that X0 is distinct from X1 and the two are not diametrically opposite)? – user189035 Mar 12 '16 at 13:05
• @user189035 Pick them 90 degrees apart. When they are very close, their cross product is numerically uncertain: there is a lot of cancellation in the subtractions, leading to loss of significant figures. – whuber Mar 13 '16 at 4:07