Calculate the distance between a point and a line segment on an ellipsoid (or WGS84 coordinates)?

Question

Hi,

I have a line segment AB between two WGS84 coordinates on the earth's surface (or, alternatively, on an ellipsoid), and need to calculate the distance between a point P and the point nearest to P that is on line segment AB. How can I do this?

Any help is greatly appreciated.

Regards, Jochen

Answer 1

I've never dealt with WGS84 coordinates and am rusty in this type of math, but I'll give you what came to mind. I'm not really answering this question, but this was way too much to put in a comment, and it has psudo code. I just thought that the problem was interesting and wanted to play around with it a little bit. Anyway, here goes...

Imagine the set of all spheres whose center is P.

Now, only one of those sphere's should share exactly 1 point with AB. If you have an equation describing AB then you should be able to find the sphere centered at that shares exactly one point with AB.

There may be a swifter way to do this than trial and error, but a binary search is what comes to mind. This should find the straight line distance to AB. If you want the disatance between P and AB I'll sorta cover that after the code. psudocode:

   Radius_lo = 0
   Radius_hi = minimum( distance(P, A), distance(P, B) )
   Radius_test = Raduis_hi // you may want to add a miniscule amount to this to keep from floating point inprecision from effectively rounding down which could lead to no actual intersection
   while true
      intersections = intersection_points( sphere( P, Radius_test), AB )
      if ( count(intersections) == 1 ) // or close enough.  With floating pt math you'll likely need to recognize close enough or it will never exit.
            return  Radius_test
      if ( count(intersections) == 2 )
            // Here you may do one of two things, both of which are binary searches.
            // 1. The first and simplest would be to divide average _test with _lo and try again
            Radius_hi = Radius_test
            Radius_test = (Radius_test + Radius_lo) / 2
            Continue
            // 2. The second would be to attempt to divide the segment fromed by the two intersection points, which is more complicated, but would almost always result in fewer times through the loop
            X = midpoint(intersection[0], intersection[1]) // midpoint would have to know how to find the midpoint of them along the elipse that they live on
            Radius_test = distance( P, X)
            Continue
     if ( count(intersections) == 0)
           // Gone too far.  If you had done (1) above, then do
           Radius_lo = Radius_test
           Radius_test = (Radius_test + Radius_hi) / 2
           Continue
           // If on the otherhand you had done (2) above then you shouldn't ever have this case
           error("0 intersections!")
     // Now, this should be for any intersection count other than 0, 1, or 2.  This shouldn't happen so
     error("bad intersection count")
  endloop // while true

This will find the straight line distance between P and AB as long as the eliptical segment AB is closer to P than any other part of the elipse that AB lies on. If this isn't true it should be easy enough to test for that and just use the closer of A or B as the closest point.

Ok, so you actually wanted the surface distance between P and AB. This is more complicated, but you can likely alter my algorithm to work with that.