[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:  PostGIS 2.2.0
Component:  postgis  |     Version:  trunk        
 Keywords:           |  
---------------------+------------------------------------------------------

Comment(by mwtoews):

 Here are some more performance and accuracy tests on Debian x86_64. A
 suite of 500k points is
 [http://geographiclib.sourceforge.net/html/geodesic.html#testgeod
 available here]. Descriptions of the fields are provided in the link, and
 most fields are regarded to be either exact, or accurate to 10^-18 .

 The `GeodTest.dat.gz` file can be `gunzip`'ed and loaded into a PostGIS
 database:
 {{{
 create unlogged table testgeod (
   id serial primary key,
   lat1 numeric, -- degrees
   lon1 numeric, -- always 0
   azi1 numeric, -- azimuth, in degrees
   lat2 numeric,
   lon2 numeric,
   azi2 numeric,
   s12 numeric,  -- distance, in metres
   a12 numeric,
   m12 numeric,
   "S12" numeric  -- area under the geodesic, in square metres
 );

 copy testgeod (lat1, lon1, azi1, lat2, lon2, azi2, s12, a12, m12, "S12")
 from '/path/to/GeodTest.dat'
 delimiter ' ';

 alter table testgeod add column geog1 geography(Point,4326);
 alter table testgeod add column geog2 geography(Point,4326);
 alter table testgeod add column distance double precision;
 alter table testgeod add column azimuth double precision;

 update testgeod set
   geog1 = ST_MakePoint(lon1, lat1)::geography,
   geog2 = ST_MakePoint(lon2, lat2)::geography,
   azimuth = radians(azi1),
   distance = s12;
 }}}

 == Direct geodesic calculations ==

 {{{
 -- test performance of st_project
 drop table if exists foo;
 select id, st_project(geog1, s12, azimuth) into unlogged foo from
 testgeod;

 -- test accuracy of st_project
 select f.id, lat1, lon1, lat2, lon2,
   st_y(st_project::geometry) as st_project_lat, st_x(st_project::geometry)
 as st_project_lon,
   s12, azi1, st_distance(st_project, geog2) AS st_distance_error
 from foo f
 join testgeod g on g.id=f.id
 order by st_distance(st_project, geog2) desc
 limit 10;

 -- percent below 10 nm accuracy
 select (count(*)/5e5*100)::numeric(8,4) || '%'
 from foo f
 join testgeod g on g.id=f.id
 where st_distance(st_project, geog2) < 1e-8
 }}}

 For the performance, compare 2490 ms for current, and 3423 ms for patch,
 or about 1.37 slower.

 For precision of distance in comparison with an exact number, compare a
 maximum of 0.002373693 m to 1.3e-08 m. Looking at % that is less than 10
 nm accuracy, compare 24.8922% to 99.9522%. At the micrometer and
 millimetre scale, the current routines have accuracy for 28.0476% and
 96.3272% of the 500k points, respectively, while !GeographicLib has 100%
 accuracy at the 13 nm level. Considering only a small performance hit,
 there is a gain of precision with !GeographicLib for those that care for
 sub millimetre accuracy.

 I'll continue on tomorrow for a similar comparison with the inverse
 geodesic functions, ST_Azimuth and ST_Distance. I'm still figuring out how
 to test S12 for ST_Area.

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