[postgis-tickets] [PostGIS] #2918: Use GeographicLib functions for ST_Azimuth, ST_Distance and ST_Project

[PostGIS]

#2918: Use GeographicLib functions for ST_Azimuth, ST_Distance and ST_Project
---------------------+------------------------------------------------------
 Reporter:  mwtoews  |       Owner:  pramsey
     Type:  patch    |      Status:  new    
 Priority:  medium   |   Milestone:         
Component:  postgis  |     Version:  trunk  
 Keywords:           |  
---------------------+------------------------------------------------------

Comment(by mwtoews):

 I haven't done any tests yet, but I'll design a few. Basically they would
 find the azimuth and distances from a cross product of point grids over a
 few extents: global, continental, and regional. I'll try a few latitudes
 to see if one algorithm converges quicker than another, because both old
 and new algorithms are iterative.

 The implementation works on different spheroid definitions, such as Mars,
 e.g. comparing equatorial antipodes:

 {{{
 SELECT ST_Distance_Spheroid(A, B, earth)/1000.0 AS earth_km,
 ST_Distance_Spheroid(A, B, mars)/1000.0 AS mars_km
 FROM (
   SELECT 'POINT(-90 0)'::geometry AS A, 'POINT(90 0)'::geometry AS B,
     'SPHEROID["WGS
 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]]'::spheroid AS earth,
     'SPHEROID["Mars_2000_IAU_IAG",3396190.0,169.89444722361179]'::spheroid
 AS mars
 ) AS f;

      earth_km     |     mars_km
 ------------------+------------------
  20003.9314586254 | 10638.0685065677
 (1 row)
 }}}

 One limitation, however, is that the geography type has a hard-coded
 spheroid, so this doesn't work:
 {{{
 INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text,
 srtext) values ( 949900, 'iau2000', 49900, '+proj=longlat +a=3396190
 +b=3376200 +no_defs ', 'GEOGCS["Mars
 2000",DATUM["D_Mars_2000",SPHEROID["Mars_2000_IAU_IAG",3396190.0,169.89444722361179]],PRIMEM["Greenwich",0],UNIT["Decimal_Degree",0.0174532925199433]]');
 SELECT degrees(ST_Azimuth(A, B)), ST_Distance(A, B)/1000.0 AS distance_km
 FROM (
   SELECT
     'SRID=949900;POINT (-90 0)'::geography AS A,
     'SRID=949900;POINT (90 0)'::geography AS B
 ) AS f;
 }}}
 Gives Earth numbers, and not Martian numbers.

-- 
Ticket URL: <http://trac.osgeo.org/postgis/ticket/2918#comment:3>
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.