# 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.