ST_Azimuth and ST_Distance on geography antipodes or near-antipodes

[mwtoews]

I've been looking at ​this question, and have been thinking that there is something wrong with PostGIS' calculation of azimuth and/or distance on an oblate spheroid.

First consider the azimuth from the south pole to the north pole. There are actually an infinite number of directions with equal distances:

SELECT degrees(ST_Azimuth(A, B)), ST_Distance(A, B)/1000.0 AS distance_km
FROM (
  SELECT 'POINT (-90 -90)'::geography AS A, 'POINT (90 90)'::geography AS B
) AS f;
-[ RECORD 1 ]-----------------
degrees     | 90
distance_km | 20003.9314586236

Sure 90°, is just the same as any other (e.g. what direction do you go from the south pole to the north pole?). POINT (10 -90) and POINT (20 90) give 5°, suggesting that the formula is something like half of the differences in longitude, which appears a bit strange. Slightly more consistent would be to just find the "average" longitude, or always use 0° (from N). Lastly, distance is spot on (half meridional circumference), so nothing out of the ordinary.

Second, consider antipodes on the equator. There should be exactly two azimuth solutions to this: north or south (0° or 180°), each spaced approx 20003.93 km apart (half meridional circumference), however I don't see this:

SELECT degrees(ST_Azimuth(A, B)), ST_Distance(A, B)/1000.0 AS distance_km
FROM (
  SELECT 'POINT (-90 0)'::geography AS A, 'POINT (90 0)'::geography AS B
) AS f;
-[ RECORD 1 ]-----------------
degrees     | 270
distance_km | 19903.5933909347

The answer is the same for any pairs of longitude spaced 180° apart on the equator (y=0). I was expecting either 0° or 180°, and a distance of approx 20003.9314586236 km. I don't know where a distance of 19903 would have come from, as it is inconsistent with an azimuth of a western-bearing, which would have used the half equatorial circumference of approx 20037.5085 km. It seems that there could be an axis order issue in the algorithm (but I haven't looked).

Thirdly, as the antipodes on the equator have exactly two possible azimuth solutions, take near-antipodal points that are slightly above or below the equator. These should yield exactly one azimuth solution right? Moving point A by a few fractions of a degree;

SELECT degrees(ST_Azimuth(A, B)), ST_Distance(A, B)/1000.0 AS distance_km
FROM (
  SELECT 'POINT (-90 0.0000001)'::geography AS A, 'POINT (90 0)'::geography AS B
) AS f;
-[ RECORD 1 ]-----------------
degrees     | 336.928683929735
distance_km | 20003.9314475662

The expected bearing should be 0°. The distance jumped about 100 m, even though I've only moved the point about 0.011 m. This expected distance should be the same half meridional circumference.

Using POSTGIS="2.1.1 r12113" GEOS="3.4.2-CAPI-1.8.2 r3924" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.10.0, released 2013/04/24" LIBXML="2.7.8" LIBJSON="UNKNOWN" TOPOLOGY RASTER

[mwtoews]

Oh, and a fourth case to consider is any non-pole antipode pair, which should have the same results as the equator example, that is two azimuth solutions 0 or 180. I'm not sure what the expected distance should be. It could be half the meridional circumference, but it might be some other number close to it.

[mwtoews]

This function seems to get the correct answers: ​http://geographiclib.sourceforge.net/html/C/geodesic_8h.html#a19bc3d000428010ad9d8509174e672c9

Perhaps I'll investigate further to see what's ticking underneath in GeographicLib, and see what happens if I rewrite spheroid_direction and spheroid_distance functions in lwspheroid.c

[mwtoews]

The geographic functions in GeographicLib are published here:

3. F. F. Karney, Algorithms for geodesics, J. Geodesy 87(1), 43–55 (Jan. 2013); DOI:​10.1007/s00190-012-0578-z

With addenda ​here.

The license is MIT/X11, so some of the relevant algorithms can be brought over to lwspheroid.c with attribution to Charles Karney.

Further note that the current implementation of algorithms is based on ​Geocentric Datum of Australia Technical Manual, with ​a web app calculator. However, the web app also has the correct answers the this bug report, so I'm not sure which one is best, except to say that the PostGIS implementation is definitely incorrect.

[mwtoews]

Ok, I think I've unearthed the substance of this ticket. From (Vincenty 1975a):

The inverse formulae may give no solution over a line between two nearly antipodal points. This will occur when λ is greater than π in absolute value.

And similar, from Karney (2013):

Vincenty’s method fails to converge for nearly antipodal points.

Furthermore:

Vincenty (1975a), who uses the iterative method of Helmert (1880, §5.13) to solve the inverse problem, was aware of its failure to converge for nearly antipodal points. In an unpublished report (​Vincenty 1975b), he gives a modification of his method which deals with this case. Unfortunately, this sometimes requires many thousands of iterations to converge, whereas Newton’s method as described here only requires a few iterations.

[pramsey]

Since we moved to Karney's algorithm's in #2918, this should resolve over time as people move to Proj 4.9