[Proj] Intersection test for two geodesics goes wrong for very long distances

Charles Karney charles.karney at sri.com
Tue Feb 5 11:10:36 EST 2013


Here are a few quick observations...

(1) For an oblate ellipsoid, a geodesic starting at lat1, lon1, azi1
stops being a shortest path when the spherical arc length is 180d.  The
Geod utility (with the -a flag) can compute this point with

   echo lat1 lon1 azi1 180 | Geod -a

(2) The geodesic with the supplementary azimuth 180-azi1 stops being a
short geodesic at the same point.

(3) For a given lat1 lon1 the locus of "end points" as azi1 is varied
[-180,180] is a line segment at latitude -lat1 centered on lon1+180d.
The end-points of this segment are given by azi1 = +/- 90d.  (The jargon
for this locus in the differential geometry community is "the cut

(4) The intersection point for two geodesics (defined by their starting
points and azimuths) might be at a point where neither of the resulting
geodesic segments is a shortest path.  Here's an example with an
ellipsoid f = 1/2 (presumably similar examples occur for f = 1/300)

              lat1 lon1  azi1   lat2 lon2 azi2
geodesic a    0  -180  89.9999  0    0  89.9999
geodesic b    0   180 -89.9999  0    0 -89.9999

You can compute these geodesics with, e.g. (the -E is needed because f
is not small),

   echo 0 -180 89.9999 360 | Geod -a -e 6.4e6 1/2 -E
   0.00000000 0.00000000 89.99990000

Neither of these geodesics is a shortest path, because an equatorial
geodesics ceases being the shortest path when the longitude difference
exceeds (b/a)*180d = 90d (in this example).

(5) I've only thought about the intersection problem in the simple case
where the end points and the intersection point all lie in a single
[shd]emi-ellipsoid (for f small).  In that case, the ellipsoidal
gnomonic projection allows you to find the intersection point in a few
iterations (I think the convergence is quadratic, a la Newton).

(6) I suspect that solving the general problem will end up being fairly


On 2013-02-05 08:03, Mikael Rittri wrote:
> Hello all,
> I have had a strange problem with geodesics. While I think I have
> solved it in a practical sense, by saying: "don't go too near the
> antipode", I wonder exactly how near the antipode it is safe to go,
> and why.
>      I post here in the hope that someone would enjoy working on
> this kind of problem (hi, Charles).
> The problem was to find the intersection of two radials, where
> a radial is a geodetic line that starts at a given point with
> a given azimuth, but doesn't have a definite end point.
> While Carmenta Engine doesn't support this computation directly,
> there is support for intersecting two geodetic lines defined by
> their start and end points. So, my idea was to compute an arbitrary
> end point for each radial, by going forward long enough on the radial.
> (This is just the forward or direct geodetic problem.)
> The question then was, how far to go forward? Obviously, one cannot
> go longer than 180 degrees or about 20 000 km on the Earth, since
> there would then be a shorter geodetic route from start to end on
> the backside of the Earth. In fact, since we want to compute on
> an ellipsoid, one cannot go quite 180 degrees either. I thought the
> maximal safe distance would be  pi * a * (1 - f), where a is the
> equatorial radius and f is the flattening. On the WGS84 ellipsoid,
> this would be pi * 6378137 * (1 - 1/298.257223563) = 19 970 326 meters.
> (Because, if one starts somewhere on the equator, and goes east
>   along the equator to an end point longer than this distance, but
>   shorter than 180 degrees, then the shortest route from start to end
>   doesn't follow the equator but takes a shortcut nearer a pole,
>   exploiting the flattening of the ellipsoid.)
> So, we tried to let "forward distance" = 19 950 000 meters, and it seemed
> to work for a while. But then we found an example where the intersecting
> method wrongly claimed that there was no intersection. It turned out
> that this wasn't caused by numerical instability, but by something more
> fundamental.
> You see, we have some initial tests in the intersecting method, which
> should say when there is definitely no intersection. One condition
> is like this, in ASCII graphics:
> B1
> | .
> |   .
> |     .
> |       . A2             A1
> |         -----------------
> |       .
> |     .
> |   .
> | .
> |
> B2
> One geodesic goes from A1 to A2, the other goes from B1 to B2. From
> A2, take the azimuths toward A1, B1 and B2. Then if the azimuth
> toward A1 doesn't go between the two azimuths toward B1 and B2,
> an intersection would be impossible, I thought. Because A1 and A2
> would be on the same side of the geodesic B1B2.
>     Of course, one has to define "between-ness" of azimuths modulo
> 360 degrees. But the idea is that the "inside" of the angle B1 A2 B2
> would be the side that is less than 180 degrees, and if the line
> from A2 toward A1 doesn't go on the inside, there can be no intersection.
> I still feel that this _ought_ to be true. However, here is a counterexample:
> The first radial is defined by
> A1 = 17.58139°S 40.17167°W
> A1_faz = -90° // A1 forward azimuth
> The second radial is defined by
> B1 = 22.23556°S 49.91139°W
> B1_faz = 30° // B1 forward azimuth
> They do intersect quite near the start points.
> However, with "forward distance" = 19 950 000 m, the forward computation
> gives the end points:
> A2 = 17.58126738°N 140.62400233°E
> B2 = 22.60088604°N 129.58209305°E
> Now,
> Azimuth at A2 toward B1 = 117.78924816°
> Azimuth at A2 toward B2 = -62.47659403°
> And 117.78924816° - (-62.47659403°) = 180.26584219° which is larger than 180°.
> This means that the inside of the angle B1 A2 B2 is south of the angle edges.
> But
> Azimuth at A2 toward A1 = 90.06652107°
> which is north of the A2-toward-B1 azimuth. So, the initial test will say
> (wrongly) that there is no intersection.
> In this particular example, I have found that the code finds the intersection
> if I lower the forward distance to 19 940 632 m. But there must be some worst
> case that requires an even shorter distance, and I don't know how to find
> the worst case (except by a lot of trial-and-error).
> I have visualized the counterexample in several map projections and waited
> for my brain to say "Aha!", but in vain. I think the intuition behind the
> initial test is that the line A1A2 is on one (oblique) hemisphere, and the
> triangle B1B2A2 is on the other hemisphere. It feels right, but clearly it
> isn't true.
> The pragmatic solution seems to be to use a much shorter forward distance,
> something like 15 000 000 m perhaps, but I would really like to understand
> what's going on here.
> Credit: the values above have been computed on the GeographicLib site,
> http://geographiclib.sourceforge.net/cgi-bin/Geod
> Regards,
> Mikael Rittri
> Carmenta
> Sweden
> http://www.carmenta.com

Charles Karney <charles.karney at sri.com>
SRI International, Princeton, NJ 08543-5300

Tel: +1 609 734 2312
Fax: +1 609 734 2662

More information about the Proj mailing list