Geodetic musings

This is a litle benchmark of three Python libraries to compute geodesic distances:

  • Pyproj (with Proj4 v.4.9.1).
  • Pygc.
  • GeographicLib.

Pyproj

Geodesic distance computation with Pyproj (python interface to PROJ4 C library).

https://github.com/jswhit/pyproj

PROJ4 C library (>= v.4.9.0) routines used to compute geodesic distances are a simple transcription from C. Karney Geographiclib C++ Library.

Geodesic distance calculations using Charles Karney geodesic algorithms: C. F. F. Karney, Algorithms for Geodesics, J. Geodesy 87(1), 43–55 (Jan. 2013)

https://trac.osgeo.org/proj/browser/trunk/proj/src/geodesic.h

PROJ4 C library (< v.4.9.0) routines used to compute geodesic distances are: Paul D. Thomas, Spheroidal Geodesics, Reference Systems, and Local Geometry. U.S. Naval Oceanographic Office, p. 162 Engineering Library 526.3 T36s (1970)

Default values used (WGS84): https://github.com/jswhit/pyproj/blob/master/lib/pyproj/__init__.py

  • Semi-major axis: 6378137.0
  • Inverse flattening: 298.257223563

Therefore:

  • Semi-minor axis: 6356752.314245179
  • Flattening: 0.0033528106647474805
In [2]:
from pyproj import Geod
In [3]:
def grcrcl1(lon_1, lat_1, lon_2, lat_2):
    
    g = Geod(ellps='WGS84')
    
    az, az_inv, dist = g.inv(lon_1, lat_1, lon_2, lat_2)
    
    return dist

Pygc

Geodesic distance computation with Pygc.

https://github.com/axiom-data-science/pygc

Geodesic distance calculations using Vincenty's formulae: http://en.wikipedia.org/wiki/Vincenty%27s_formulae

Default values used (WGS84): https://github.com/axiom-data-science/pygc/blob/master/pygc/gc.py

  • Semi-major axis: 6378137.0
  • Semi-minor axis: 6356752.3142

Therefore:

  • Flattening: 0.0033528106718309896
  • Inverse flattening: 298.2572229328697
In [4]:
from pygc import great_distance
In [5]:
def grcrcl2(startlong, startlat, endlong, endlat):
    
    res = great_distance(start_latitude=startlat, start_longitude=startlong, 
                        end_latitude=endlat, end_longitude=endlong)
    
    return float(res['distance'])

GeographicLib

Geodesic distance computation with GeographicLib (python interface to GeographicLib C++ library).

http://geographiclib.sourceforge.net/html/other.html#python

Geodesic distance calculations using Charles Karney geodesic algorithms: C. F. F. Karney, Algorithms for Geodesics, Journal of Geodesy 87(1), 43–55 (Jan. 2013)

Default values used (WGS84): http://geographiclib.sourceforge.net/html/Constants_8hpp_source.html

  • Semi-major axis: 6378137.0
  • Inverse flattening: 298.257223563

Therefore:

  • Semi-minor axis: 6356752.314245179
  • Flattening: 0.0033528106647474805
In [6]:
from geographiclib.geodesic import Geodesic
In [7]:
def grcrcl3(lon_1, lat_1, lon_2, lat_2):
    
    res = Geodesic.WGS84.Inverse(lat_1,lon_1,lat_2,lon_2)
    
    return res['s12']

Results

Function to print results from geodesic distance functions:

In [8]:
def printResults(input_coords):
    
    lon_1, lat_1, lon_2, lat_2 = input_coords
    
    dist1 = grcrcl1(lon_1, lat_1, lon_2, lat_2)
    dist2 = grcrcl2(lon_1, lat_1, lon_2, lat_2)
    dist3 = grcrcl3(lon_1, lat_1, lon_2, lat_2)
    
    print "Distances for: lon_1 {0}, lat_1 {1}, lon_2 {2}, lat_2 {3}\n".format(lon_1, lat_1, lon_2, lat_2)
    print "Pyproj___________________________{:,.10f} m\n".format(dist1)
    print "Pygc_____________________________{:,.10f} m\n".format(dist2)
    print "GeographicLib____________________{:,.10f} m\n".format(dist3)
    print "Difference Pyproj,Pygc___________{:.15f} m\n".format(abs(dist1 - dist2))
    print "Difference Pyproj,GeographicLib__{:.15f} m\n".format(abs(dist1 - dist3))
    print "Difference GeographicLib,Pygc____{:.15f} m\n".format(abs(dist3 - dist2))
    print "------------------------------------------\n"

Creating test data. Test data is a list of lists. Each inner list contains:

  • Start longitude.
  • Start latitude.
  • End longitude.
  • End latitude.
In [9]:
data = [
        [-3.6,40.5,-118.4,33.9],
        [-6.,37.,-145.,11.],
        [-150.,37.,140.,11.],
        [-50.,7.,40.,11.],
        [-100.,80.,-140.,30.]
    ]

Printing three functions results with created test data:

In [10]:
map(printResults, data)
Distances for: lon_1 -3.6, lat_1 40.5, lon_2 -118.4, lat_2 33.9

Pyproj___________________________9,406,468.4536480606 m

Pygc_____________________________9,406,468.4522715993 m

GeographicLib____________________9,406,468.4536480606 m

Difference Pyproj,Pygc___________0.001376461237669 m

Difference Pyproj,GeographicLib__0.000000000000000 m

Difference GeographicLib,Pygc____0.001376461237669 m

------------------------------------------

Distances for: lon_1 -6.0, lat_1 37.0, lon_2 -145.0, lat_2 11.0

Pyproj___________________________13,189,756.0717643406 m

Pygc_____________________________13,189,756.0719679464 m

GeographicLib____________________13,189,756.0717643406 m

Difference Pyproj,Pygc___________0.000203605741262 m

Difference Pyproj,GeographicLib__0.000000000000000 m

Difference GeographicLib,Pygc____0.000203605741262 m

------------------------------------------

Distances for: lon_1 -150.0, lat_1 37.0, lon_2 140.0, lat_2 11.0

Pyproj___________________________7,510,240.7866536863 m

Pygc_____________________________7,510,240.7865020484 m

GeographicLib____________________7,510,240.7866536863 m

Difference Pyproj,Pygc___________0.000151637941599 m

Difference Pyproj,GeographicLib__0.000000000000000 m

Difference GeographicLib,Pygc____0.000151637941599 m

------------------------------------------

Distances for: lon_1 -50.0, lat_1 7.0, lon_2 40.0, lat_2 11.0

Pyproj___________________________9,871,047.1972853243 m

Pygc_____________________________9,871,047.1963992771 m

GeographicLib____________________9,871,047.1972853225 m

Difference Pyproj,Pygc___________0.000886047258973 m

Difference Pyproj,GeographicLib__0.000000001862645 m

Difference GeographicLib,Pygc____0.000886045396328 m

------------------------------------------

Distances for: lon_1 -100.0, lat_1 80.0, lon_2 -140.0, lat_2 30.0

Pyproj___________________________5,853,723.4893777100 m

Pygc_____________________________5,853,723.4893133771 m

GeographicLib____________________5,853,723.4893777100 m

Difference Pyproj,Pygc___________0.000064332969487 m

Difference Pyproj,GeographicLib__0.000000000000000 m

Difference GeographicLib,Pygc____0.000064332969487 m

------------------------------------------

Out[10]:
[None, None, None, None, None]

Discussion

As it was expected there are no difference betwen Pyproj and GeographicLib outcomes. However there are a minimal difference between Pygc and Pyproj/Geographiclib.