[postgis-tickets] [PostGIS] #2913: ST_Azimuth and ST_Distance on geography antipodes or near-antipodes

[PostGIS]

#2913: ST_Azimuth and ST_Distance on geography antipodes or near-antipodes
---------------------+------------------------------------------------------
 Reporter:  mwtoews  |       Owner:  pramsey
     Type:  defect   |      Status:  new    
 Priority:  medium   |   Milestone:         
Component:  postgis  |     Version:  2.1.x  
 Keywords:           |  
---------------------+------------------------------------------------------
 I've been looking at [http://stackoverflow.com/q/25553230/327026 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

-- 
Ticket URL: <http://trac.osgeo.org/postgis/ticket/2913>
PostGIS <http://trac.osgeo.org/postgis/>
The PostGIS Trac is used for bug, enhancement & task tracking, a user and developer wiki, and a view into the subversion code repository of PostGIS project.