{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#Geodetic musings\n",
    "\n",
    "This is a litle benchmark of three Python libraries to compute geodesic distances:\n",
    "- Pyproj (with Proj4 v.4.8.0).\n",
    "- Pygc.\n",
    "- GeographicLib."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##Pyproj\n",
    "Geodesic distance computation with Pyproj\n",
    "(python interface to PROJ4 C library).\n",
    "\n",
    "https://github.com/jswhit/pyproj\n",
    "\n",
    "PROJ4 C library (>= v.4.9.0) routines used to compute \n",
    "geodesic distances are a simple transcription from \n",
    "C. Karney Geographiclib C++ Library.\n",
    "\n",
    "Geodesic distance calculations using Charles Karney \n",
    "geodesic algorithms:\n",
    "    C. F. F. Karney, Algorithms for Geodesics, \n",
    "        J. Geodesy 87(1), 43–55 (Jan. 2013)\n",
    "\n",
    "https://trac.osgeo.org/proj/browser/trunk/proj/src/geodesic.h\n",
    "\n",
    "PROJ4 C library (< v.4.9.0) routines used to compute \n",
    "geodesic distances are:\n",
    "    Paul D. Thomas, Spheroidal Geodesics, \n",
    "        Reference Systems, and Local Geometry.\n",
    "        U.S. Naval Oceanographic Office, p. 162\n",
    "        Engineering Library 526.3 T36s (1970)\n",
    "\n",
    "Default values used (WGS84):\n",
    "https://github.com/jswhit/pyproj/blob/master/lib/pyproj/__init__.py\n",
    "   \n",
    "   - Semi-major axis: 6378137.0\n",
    "   - Inverse flattening: 298.257223563\n",
    "   \n",
    "Therefore:\n",
    "   - Semi-minor axis: 6356752.314245179\n",
    "   - Flattening: 0.0033528106647474805"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from pyproj import Geod"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def grcrcl1(lon_1, lat_1, lon_2, lat_2):\n",
    "    \n",
    "    g = Geod(ellps='WGS84')\n",
    "    \n",
    "    az, az_inv, dist = g.inv(lon_1, lat_1, lon_2, lat_2)\n",
    "    \n",
    "    return dist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##Pygc\n",
    "Geodesic distance computation with Pygc.\n",
    "\n",
    "https://github.com/axiom-data-science/pygc\n",
    "\n",
    "Geodesic distance calculations using Vincenty's formulae:\n",
    "    http://en.wikipedia.org/wiki/Vincenty%27s_formulae\n",
    "\n",
    "Default values used (WGS84):\n",
    "https://github.com/axiom-data-science/pygc/blob/master/pygc/gc.py\n",
    "    \n",
    "   - Semi-major axis: 6378137.0\n",
    "   - Semi-minor axis: 6356752.3142\n",
    "\n",
    "Therefore:\n",
    "   - Flattening: 0.0033528106718309896\n",
    "   - Inverse flattening: 298.2572229328697"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from pygc import great_distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def grcrcl2(startlong, startlat, endlong, endlat):\n",
    "    \n",
    "    res = great_distance(start_latitude=startlat, start_longitude=startlong, \n",
    "                        end_latitude=endlat, end_longitude=endlong)\n",
    "    \n",
    "    return float(res['distance'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##GeographicLib\n",
    "Geodesic distance computation with GeographicLib \n",
    "(python interface to GeographicLib C++ library).\n",
    "\n",
    "http://geographiclib.sourceforge.net/html/other.html#python\n",
    "\n",
    "Geodesic distance calculations using Charles Karney \n",
    "geodesic algorithms:\n",
    "    C. F. F. Karney, Algorithms for Geodesics, \n",
    "    Journal of Geodesy 87(1), 43–55 (Jan. 2013)\n",
    "\n",
    "Default values used (WGS84):\n",
    "http://geographiclib.sourceforge.net/html/Constants_8hpp_source.html\n",
    "   - Semi-major axis: 6378137.0\n",
    "   - Inverse flattening: 298.257223563\n",
    "\n",
    "Therefore:\n",
    "   - Semi-minor axis: 6356752.314245179\n",
    "   - Flattening: 0.0033528106647474805"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from geographiclib.geodesic import Geodesic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def grcrcl3(lon_1, lat_1, lon_2, lat_2):\n",
    "    \n",
    "    res = Geodesic.WGS84.Inverse(lat_1,lon_1,lat_2,lon_2)\n",
    "    \n",
    "    return res['s12']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##Results\n",
    "Function to print results from geodesic distance functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def printResults(input_coords):\n",
    "    \n",
    "    lon_1, lat_1, lon_2, lat_2 = input_coords\n",
    "    \n",
    "    dist1 = grcrcl1(lon_1, lat_1, lon_2, lat_2)\n",
    "    dist2 = grcrcl2(lon_1, lat_1, lon_2, lat_2)\n",
    "    dist3 = grcrcl3(lon_1, lat_1, lon_2, lat_2)\n",
    "    \n",
    "    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)\n",
    "    print \"Pyproj___________________________{:,.10f} m\\n\".format(dist1)\n",
    "    print \"Pygc_____________________________{:,.10f} m\\n\".format(dist2)\n",
    "    print \"GeographicLib____________________{:,.10f} m\\n\".format(dist3)\n",
    "    print \"Difference Pyproj,Pygc___________{:.15f} m\\n\".format(abs(dist1 - dist2))\n",
    "    print \"Difference Pyproj,GeographicLib__{:.15f} m\\n\".format(abs(dist1 - dist3))\n",
    "    print \"Difference GeographicLib,Pygc____{:.15f} m\\n\".format(abs(dist3 - dist2))\n",
    "    print \"------------------------------------------\\n\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Creating test data. Test data is a list of lists. Each inner list contains:\n",
    " - Start longitude.\n",
    " - Start latitude.\n",
    " - End longitude.\n",
    " - End latitude."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data = [\n",
    "        [-3.6,40.5,-118.4,33.9],\n",
    "        [-6.,37.,-145.,11.],\n",
    "        [-150.,37.,140.,11.],\n",
    "        [-50.,7.,40.,11.],\n",
    "        [-100.,80.,-140.,30.]\n",
    "    ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Printing three functions results with created test data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Distances for: lon_1 -3.6, lat_1 40.5, lon_2 -118.4, lat_2 33.9\n",
      "\n",
      "Pyproj___________________________9,406,468.4398682937 m\n",
      "\n",
      "Pygc_____________________________9,406,468.4522715993 m\n",
      "\n",
      "GeographicLib____________________9,406,468.4536480606 m\n",
      "\n",
      "Difference Pyproj,Pygc___________0.012403305619955 m\n",
      "\n",
      "Difference Pyproj,GeographicLib__0.013779766857624 m\n",
      "\n",
      "Difference GeographicLib,Pygc____0.001376461237669 m\n",
      "\n",
      "------------------------------------------\n",
      "\n",
      "Distances for: lon_1 -6.0, lat_1 37.0, lon_2 -145.0, lat_2 11.0\n",
      "\n",
      "Pyproj___________________________13,189,756.0277870093 m\n",
      "\n",
      "Pygc_____________________________13,189,756.0719679464 m\n",
      "\n",
      "GeographicLib____________________13,189,756.0717643406 m\n",
      "\n",
      "Difference Pyproj,Pygc___________0.044180937111378 m\n",
      "\n",
      "Difference Pyproj,GeographicLib__0.043977331370115 m\n",
      "\n",
      "Difference GeographicLib,Pygc____0.000203605741262 m\n",
      "\n",
      "------------------------------------------\n",
      "\n",
      "Distances for: lon_1 -150.0, lat_1 37.0, lon_2 140.0, lat_2 11.0\n",
      "\n",
      "Pyproj___________________________7,510,240.7935455162 m\n",
      "\n",
      "Pygc_____________________________7,510,240.7865020484 m\n",
      "\n",
      "GeographicLib____________________7,510,240.7866536863 m\n",
      "\n",
      "Difference Pyproj,Pygc___________0.007043467834592 m\n",
      "\n",
      "Difference Pyproj,GeographicLib__0.006891829892993 m\n",
      "\n",
      "Difference GeographicLib,Pygc____0.000151637941599 m\n",
      "\n",
      "------------------------------------------\n",
      "\n",
      "Distances for: lon_1 -50.0, lat_1 7.0, lon_2 40.0, lat_2 11.0\n",
      "\n",
      "Pyproj___________________________9,871,047.2090491671 m\n",
      "\n",
      "Pygc_____________________________9,871,047.1963992771 m\n",
      "\n",
      "GeographicLib____________________9,871,047.1972853225 m\n",
      "\n",
      "Difference Pyproj,Pygc___________0.012649890035391 m\n",
      "\n",
      "Difference Pyproj,GeographicLib__0.011763844639063 m\n",
      "\n",
      "Difference GeographicLib,Pygc____0.000886045396328 m\n",
      "\n",
      "------------------------------------------\n",
      "\n",
      "Distances for: lon_1 -100.0, lat_1 80.0, lon_2 -140.0, lat_2 30.0\n",
      "\n",
      "Pyproj___________________________5,853,723.4817370065 m\n",
      "\n",
      "Pygc_____________________________5,853,723.4893133771 m\n",
      "\n",
      "GeographicLib____________________5,853,723.4893777100 m\n",
      "\n",
      "Difference Pyproj,Pygc___________0.007576370611787 m\n",
      "\n",
      "Difference Pyproj,GeographicLib__0.007640703581274 m\n",
      "\n",
      "Difference GeographicLib,Pygc____0.000064332969487 m\n",
      "\n",
      "------------------------------------------\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[None, None, None, None, None]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "map(printResults, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##Discussion\n",
    "\n",
    "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."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}