[Proj] Improvements to geod and invgeod

[Charles Karney]

geod and invgeod, programs for performing geodesic calculations, have
been part of proj.4 since at least 1995.  There are several problems
with them:

* the overall accuracy is limited;

* the results of inverse calculations are sometimes wildly wrong;

* there's no documented API for the geodesic capabilities;

* there's no provision from computing various quantities associated with
   geodesics (in particular the area of a geodesic polygon).

These problems will be addressed in the next release of proj.4 with the
inclusion of the algorithms for geodesics from GeographicLib.  These are
documented in

   C. F. F. Karney,
   Algorithms for geodesics,
   J. Geodesy 87(1), 43-55 (Jan. 2013)
   http://dx.doi.org/10.1007/s00190-012-0578-z
   Addenda: http://geographiclib.sf.net/geod-addenda.html

Points to note:

* The changes have been checked into subversion but are not part of a
   release yet.  To try them out, check out proj.4 from subversion with

     svn checkout http://svn.osgeo.org/metacrs/proj/trunk/proj

   As of revision 2334, the man page geodesic.3 is missing.  In order to
   have everything build without complaint, create a dummy one with
   something like

     touch man/man3/geodesic.3

* The interfaces to geod and invgeod are unchanged (and hence there are
   minimal changes to the man page geod.1).

* geod and invgeod now yield accurate results for all legal inputs
   provided that |f| < 1/50 where f is the flattening of the ellipsoid.

* An API for geodesic calculations is now included in the proj library.
   The header file is geodesic.h and the routine names begin with
   "geod_".  This API is (will be) briefly documented in geodesic.3;
   however this just directs you to the doxygen-generated documentation
   at

     http://geographiclib.sf.net/html/C/index.html

* There are a couple of noteworthy differences between the APIs for
   geod_ and pj_:

   + the ellipsoid is specified by the equatorial radius and flattening
     (instead of the equatorial radius and eccentricity);

   + latitudes, longitudes, azimuths are specified in degrees instead of
     radians (this allows an accurate treatment of special angles).

* The geodesic API includes most of the functionality from
   GeographicLib; this includes:

   + efficient computation of multiple points along a specific geodesic,

   + the computation of the reduced length and geodesic scale,

   + the computation of the area of geodesic polygons.

   Not included is the solution of the geodesics in terms of elliptic
   integrals (needed for the accurate solution of geodesics when |f| >
   1/50, e.g., for Jupiter).  (This is provided by the GeodesicExact
   class in GeographicLib.)

* The solutions of geodesic problems are now sufficiently fast,
   accurate, and reliable, that it is possible to implement some
   projections directly in terms of geodesics.  Obvious candidates are
   the azimuthal equidistant and Cassini-Soldner.  The scales for these
   projections can be expressed in terms of the reduced length and
   geodesic scale.

* Several libraries that depend on proj.4, e.g., GRASS, use approximate
   methods for computing polygonal areas.  They can now use proj.4 to
   compute areas accurately.

* The geodesic utility in GeographicLib used to be named Geod which
   caused a conflict with proj.4's geod utility on systems which do not
   pay attention to the case of file names.  Starting with GeographicLib
   1.30, I have started to transition the name to GeodSolve.

-- 
Charles Karney <charles.karney at sri.com>
SRI International, Princeton, NJ 08543-5300

Tel: +1 609 734 2312
Fax: +1 609 734 2662