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

 Well, I figured out the NaN issue, and is detailed in ticket #2932. The
 issue was due to a shortcut of logic that determines the sphere area if it
 is deemed too difficult to solve on an ellipsoid of revolution. I've moved
 the `USE_GEODESIC` def to `liblwgeom/liblwgeom_internal.h`, which is
 referenced in `postgis/geography_measurement.c`. The updated PR updates
 this. Should the `USE_GEODESIC` def appear elsewhere?

 The previous attempts to find the area "S12" are mixed, as it is not easy
 to define the quadrilateral AFHB using a planimeter approach, since the FH
 geodesics don't necessarily follow the equator (since the shortest path is
 sometimes across the poles for edge length (1-f)*180 deg or 179.396494
 deg). It would require `geod_gendirect` to properly determine "S12", which
 isn't useful for ST_Area. So I'm scrapping this test approach for area,
 and I don't have any more rigid way to test and compare methods except for
 the previous reports comparing country sizes. I've individually tested
 polygons with GeographicLib's Planimeter with exact option (`-E`) show
 essentially the same area as normal calculations without `-E`.

 However, with respects to GeographicLib's area calculation, it looks much
 more improved in several ways:
  1. it actually calculates on the ellipsoid of a revolution if the user
 requests it (i.e. `use_spheroid=true`, whereas the current implementation
 silently calculates on a sphere for some brittle cases)
  2. areas can be calculated around poles, and across the equator, etc.,
 and is generally much more robust for different spheroid definitions.

 E.g., area of Antarctica:
 {{{
 SELECT ST_Area('POLYGON ((-74 -72.9, -102 -71.9, -102 -74.9, -131 -74.3,
 -163 -77.5, 163 -77.4, 172 -71.7, 140 -65.9, 113 -65.7, 88 -66.6, 59
 -66.9, 25 -69.8, -4 -70, -14 -71, -33 -77.3, -46 -77.9, -61 -74.7, -74
 -72.9))'::geography);
 }}}

 previously, this would have silently short-cut to a sphere calculation,
 regardless of `use_spheroid` option. The exact calculation is
 13376856682207.375 m²

 Now what? Possible ideas: redefine ST_Length and ST_Perimeter for
 geography. Possibly look at ST_Segmentize. And round a few asserts in
 `liblwgeom/cunit/cu_geodetic.c` so they work with current and
 GeographicLib's results, since they are most likely very close, but
 certainly not exactly the same.

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