9

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.

Glorfindel
  • 1,096
  • 2
  • 9
  • 14
bouke
  • 181
  • 1
  • 9
  • 1
    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
  • possible duplicate of Calculating the distance between a point and a virtual line of two lat/lngs – Martin F Apr 25 '15 at 17:27
  • @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 http://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
  • How do you propose to "Select the closest edge" without first "Projecting" your location onto each edge? – Martin F May 18 '15 at 00:59
  • @Martin That would be the most efficient way, actually. The closest point to a great circle (GC) is most easily found by computing the 3D distances between the points and the GC: that would be just one preliminary calculation for performing the actual projections. Once the closest point is found, only then does it actually need to be projected onto the GC. Because of these computational efficiencies, the present question is not an exact duplicate of the related ones you have found. – whuber May 18 '15 at 16:30

1 Answers1

7

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.

Figure

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.

Community
  • 1
whuber
  • 69,783
  • 15
  • 186
  • 281
  • 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
  • 1
    @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 04:07