How reliable is Postgis's ST_Distance

Question

I am very confused about the measurements of ST_Distance in PostGIS. If I measure along the equator I get longer measurements using the SPHEROID instead of the sphere, but only up to 179 degrees (and this makes sense since distance on a spheroid is obviously longer). From 180 degrees + this turns around and SPHERE surpisingly yields longer results

SELECT  ST_Distance(ST_Point(0,0), ST_Point(179,0), True)/1000 AS spheroid_dist,
ST_Distance(ST_Point(0,0), ST_Point(179,0), False)/1000 AS sphere_dist;

outputs

spheroid_distance;sphere_distance
19926.1888494779;19903.9192724991

while

SELECT  ST_Distance(ST_Point(0,0), ST_Point(180,0), True)/1000 AS spheroid_dist,
    ST_Distance(ST_Point(0,0), ST_Point(180,0), False)/1000 AS sphere_dist;

outputs

spheroid_distance;sphere_distance
19903.5933909347;20015.1143522337

To make it more confusing I get different results for the spheroid measurement for different PostGIS versions. Above are results for Version 2.1

Version 2.3 yields the following results (for the second SELECT statement, i.e. 180 degrees):

spheroid_distance;sphere_distance
20003.9314586255;20015.1143522337

So what numbers are to be trusted? Can anyone shed some light?

Answer 1

To clear up the first bit of confusion, it is expected you would get different results on antipodal points with different PostGIS versions. The newer version of PostGIS has the correct result, reliable within less than 15 nanometers of precision.

Prior to PostGIS 2.2, spheroid distances were calculated using Vincenty's formulae, which works most of the time, except for nearly antipodal points where it fails to converge and you get the wrong answer (i.e. 19903.5933909347 is wrong and 20003.9314586255 is correct). PostGIS 2.2 and later use GeographicLib's inverse geodesic routine (see Karney, C.F.F. J Geod (2013)), which gets the correct answer with better precision in the sub-millimetre range.

---

The shape of the Earth is essentially an oblate spheroid, having a slightly smaller diameter from N–S pole (height), than the diameter along the equator (width).

With the first example from Null Island to 179°E on the equator, the shortest distance has a bearing eastward along the equator:

SELECT degrees(ST_Azimuth(A, B)),
    ST_Distance(A, B, true) > ST_Distance(A, B, false) spheroid_gt_sphere
FROM (SELECT ST_Point(0, 0)::geography A, ST_Point(179, 0)::geography B) d;
-[ RECORD 1 ]------+---
degrees            | 90
spheroid_gt_sphere | t

With the second example, the second point is on the equator at 180°E (or 180°W—same place) is an antipodal point, which is a rare but complicated problem. On a sphere, the two points have an infinite number of bearings between each other with the same shortest distance. However, on an oblate spheroid, there are exactly two solutions: traverse over N or S pole, because the distance is slightly shorter than along the equator. PostGIS 2.2 (or later) will only return one of these bearings:

SELECT degrees(ST_Azimuth(A, B)),
    ST_Distance(A, B, true) > ST_Distance(A, B, false) spheroid_gt_sphere
FROM (SELECT ST_Point(0, 0)::geography A, ST_Point(180, 0)::geography B) d;
-[ RECORD 1 ]------+--
degrees            | 0
spheroid_gt_sphere | f

Evaluating the azimuth between two points is important to understand why there are different distances between the sphere and spheroid methods. With the first, the distance is along the equator, which is wider for an oblate spheroid than a sphere. While the distance along the poles is shorter than on a sphere.