Update geodesic routines

[cffk]

Update geodesic routines

Many moons ago (2012-12-21, to be precise), Paul Wessel had expressed some interest in improving the calculation of geodesics in GMT 6. Perhaps, it's time to revisit this.

GMT already depends (indirectly) on PROJ (via GDAL) and PROJ since version 5.0 incorporates my geodesic algorithms. So it would be relatively straightforward for GMT to use these. The benefits would be

• Distance calculation always converges
• Accurate to round-off
• Interface for computing way-points efficiently
• Accurate computation of ellipsoidal area of polygons
• Access to differential quantities (might be useful in computing buffers)

Are you willing to help implement and maintain this feature?

I'm not a stakeholder in GMT so it would probably be best if someone else make the code changes. However, I'm happy to pitch in if necessary.

If you can make PROJ 5.0 or later a requirement, that you would be able to remove code from GMT.

[PaulWessel]

Hi @cffk, no question this would simplify life in GMT. It has been on our perpetual to-do list for ages, but will require serious redesign. We will discuss this at an upcoming GMT community meeting to see how we might plan for this type of change.

[cffk]

OK. Let me know if I can help...

At the simplest level (shoehorning the PROJ4 version in as a direct substitute for Vincenty), it would involve

  CMakeLists.txt:
    find_package(PROJ4)
    set (HAVE_PROJ4 TRUE CACHE INTERNAL "System has PROJ4")
    include_directories (${PROJ4_INCLUDE_DIRS})
    list (APPEND GMT_OPTIONAL_LIBRARIES ${PROJ4_LIBRARIES})

  gmt_map.c:
    #ifdef HAVE_PROJ4
    #include <proj.h>
    #if PROJ_VERSION_MAJOR > 4
    #include <geodesic.h>
    #define HAVE_GEODESIC 1
    #endif
    #endif

    ...
    #if HAVE_GEODESIC
    /* Define
         gmtmap_vincenty_dist_meter
         gmtmap_az_backaz_vincenty
         gmtmap_translate_point_geodesic
       in terms of the geodesic routines in PROJ
     */
    #else
    /* Existing code */
    #endif

However, I realize that various issues will need to be worked out first:

• An supplementary name for the PROJ version, I suggest GMT_GEODESIC_ELLIPSOID (expecting the VINCENTY, RUDOE, and ANDOYER versions to be retired at some point).

• PROJ4 is transitioning to PROJ for the package name.

• Having CMake get the PROJ version number. This is easy with find_package(PROJ4 CONFIG), but might require examining proj.h with a MODULE lookup.

• Exposing the calculation of areas for geodesic polygons.