Geodetic musings

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

  • Pyproj (with Proj4 v.4.8.0).
  • 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 [1]:
from pyproj import Geod
In [4]:
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 [2]:
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 [3]:
from geographiclib.geodesic import Geodesic
In [6]:
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 [7]:
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 [8]:
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 [11]:
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.4398682937 m

Pygc_____________________________9,406,468.4522715993 m

GeographicLib____________________9,406,468.4536480606 m

Difference Pyproj,Pygc___________0.012403305619955 m

Difference Pyproj,GeographicLib__0.013779766857624 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.0277870093 m

Pygc_____________________________13,189,756.0719679464 m

GeographicLib____________________13,189,756.0717643406 m

Difference Pyproj,Pygc___________0.044180937111378 m

Difference Pyproj,GeographicLib__0.043977331370115 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.7935455162 m

Pygc_____________________________7,510,240.7865020484 m

GeographicLib____________________7,510,240.7866536863 m

Difference Pyproj,Pygc___________0.007043467834592 m

Difference Pyproj,GeographicLib__0.006891829892993 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.2090491671 m

Pygc_____________________________9,871,047.1963992771 m

GeographicLib____________________9,871,047.1972853225 m

Difference Pyproj,Pygc___________0.012649890035391 m

Difference Pyproj,GeographicLib__0.011763844639063 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.4817370065 m

Pygc_____________________________5,853,723.4893133771 m

GeographicLib____________________5,853,723.4893777100 m

Difference Pyproj,Pygc___________0.007576370611787 m

Difference Pyproj,GeographicLib__0.007640703581274 m

Difference GeographicLib,Pygc____0.000064332969487 m

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

Out[11]:
[None, None, None, None, None]

Discussion

There are minor differences in outcomes though curiously the most puzzling is that the highest dissimilarity occur between Pyproj and GeographicLib which is using the same algorithms in newer versions (>= 4.9.0). There are a minimal difference between Pygc and Geographiclib.