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

Tue Feb 5 11:10:36 EST 2013

Mikael:

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
locus".)

(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
complicated.

   --Charles

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