[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):

 Performance tests are looking just over twice as slow for the geodesic C
 functions from GeographicLib compared to the current algorithms in
 PostGIS, equally for `ST_Azimuth` and `ST_Distance` which call the same
 function `geod_inverse`. I've also tested modifications to `geodesic.c`'s
 `#define GEOGRAPHICLIB_GEODESIC_ORDER` from 6 to 3, which speeds things up
 with less precision, but they are still less than twice as slow as the
 existing functions.

 The geodesic length and area functions are not tested to compare, but I'll
 investigate only if PSC think it is worth doing. There could be
 improvements to speed by changing internal coordinate units for
 `GEOGRAPHIC_POINT` in degrees, rather than constantly flipping them to
 radians (I'm yet to find any geodesic tool that uses radians), but this
 would require broader internal changes.

 Some of the commands used to do some testing are provided here:
 {{{
 -- Sanity check to see if PostGIS's or GeographicLib's functions are being
 used
 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;

 CREATE OR REPLACE FUNCTION ST_GriddedPoints(
         nlon integer, nlat integer,
         lonsize float8, latsize float8,
         lon0 float8, lat0 float8)
     RETURNS SETOF geography AS
 'SELECT ST_MakePoint(i / ($1::float8 - 1) * $3 + $5, j / ($2::float8 - 1)
 * $4 + $6)::geography
 FROM
   generate_series(0, $1 - 1) AS i,
   generate_series(0, $2 - 1) AS j' LANGUAGE sql IMMUTABLE STRICT;

 CREATE TEMP TABLE pts (gid serial primary key, geog
 geography(Point,4326));

 TRUNCATE pts;
 INSERT INTO pts(geog)
 SELECT g
 FROM ST_GriddedPoints(31, 21, 10, 10, 0, 0) AS g;

 SELECT ST_Azimuth(p1.geog, p2.geog), ST_Distance(p1.geog, p2.geog)
 INTO UNLOGGED foo
 FROM pts p1, pts p2;
 drop table foo;
 }}}
 Also, I've only tested this on a Debian x64-based computer. I haven't
 tested on any other platform, but would be keen to hear the comparisons.

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