TL;DR - Vincenty's inverse formula provides an accurate method for calcualting the distance between two latitude/longitude pairs. Learn how to implement your own with Python.

Recently the work I have been doing requires a higher degree of accuracy than which the haversine method I was using could provide. This prompted me to implement a Python version of the Vincenty's inverse formula. Unlike the Haversine method (which I posted about previously) of directly calculating the great-circle distance between two points on a perfectly spherical Earth, Vincenty's formulae is an iterative method which more realistically assumes Earth as an oblate spheroid.

Before proceeding further, several values need to be defined:

VariableDefinition
$$a = 6378137.0$$length of semi-major axis of the ellipsoid aka radius of the equator in meters (WGS-84)
$$f = 1 / 298.257223563$$flattening of the ellipsoid (WGS-84)
$$b = \left ( 1 - f \right ) a$$length of semi-minor axis of the ellipsoid (radius at the poles, 6356752.314245 meters in WGS-84)
$$\phi_1, \phi_2$$latitude of the points
$$U_1 = \arctan \left [ \left ( 1 - f \right ) \tan \phi_1 \right ]$$reduced latitude (latitude on the auxiliary sphere)
$$U_2 = \arctan \left [ \left ( 1 - f \right ) \tan \phi_2 \right ]$$reduced latitude (latitude on the auxiliary sphere)
$$L = L_2 - L_1$$difference in longitude of two points
$$\lambda_1, \lambda_2$$longitude of the points on the auxiliary sphere
$$\alpha_1, \alpha_2$$forward azimuths at the points
$$\alpha$$azimuth at the equator
$$s$$ellipsoidal distance between the two points
$$\sigma$$arc length between points on the auxiliary sphere

Using the values defined above, and the eight equations below, the process of iterating towards a solution can begin. Starting with the first equation on top, sequentially progress towards the last equation at the bottom where the output is $$\lambda$$. Once through the first iteration, simply use this newly calculated value of $$\lambda$$ and dump it back into the first equation and repeat the process. The end game here is to minimize the value of $$\lambda$$.

$\sin \sigma = \sqrt{(\cos U_2 \sin\lambda)^2 + (\cos U_1 \sin U_2 - \sin U_1 \cos U_2 \cos\lambda)^2}$

$\cos \sigma = \sin U_1 \sin U_2 + \cos U_1 \cos U_2 \cos \lambda$

$\sigma = \arctan \frac{\sin \sigma}{\cos \sigma}$

$\sin \alpha = \frac{ \cos U_1 \cos U_2 \sin \lambda }{ \sin \sigma }$

$\cos^2 \alpha = 1 - \sin^2 \alpha$

$\cos \left( 2 \sigma_m \right) = \cos \sigma - \frac{2 \sin U_1 sin U_2 }{ \cos^2 \alpha }$

$C = \frac{f}{16} \cos^2 \alpha \left[ 4 + f \left( 4 - 3 \cos^2 \alpha \right) \right]$

$\lambda = L + (1 - C) f \sin \alpha \left \{ \sigma + C \sin \sigma \left [ \cos \left ( 2 \sigma_m \right) + C \cos \sigma \left ( -1 + 2\cos^2 \left ( 2 \sigma_m \right ) \right ) \right ] \right \}$

When the difference between the current value of $$\lambda$$ and the value of $$\lambda$$ from the previous iteration is less than the convergence tolerance, we can move on to the final stage (a convergence tolerance of 10e-12 translates to an accuracy of approximately 0.06mm).

$u^2 = \cos^2 \alpha \frac{a^2 - b^2}{b^2}$

$A = 1 + \frac{u^2}{16384} \left \{ 4096 + u^2 \left [ -768 + u^2 \left ( 320 - 175 u^2 \right ) \right ] \right \}$

$B = \frac{u^2}{1024} \left \{ 256 + u^2 \left [ -128 + u^2 \left ( 74 - 47 u^2 \right ) \right ] \right \}$

$\Delta \sigma = B \sin \sigma \left \{ \cos \left ( 2 \sigma_m \right ) + \frac{1}{4} B \left [ \cos \sigma \left ( -1 + 2 \cos^2 \left ( 2 \sigma_m \right ) \right ) - \frac{1}{6} B \cos \left ( 2 \sigma_m \right ) \left ( -3 + 4 \sin^2 \sigma \right ) \left ( -3 + 4 \cos^2 \left ( 2 \sigma_m \right ) \right ) \right ] \right \}$

$s = b A \left ( \sigma - \Delta \sigma \right )$

$\alpha_1 = \arctan \left ( \frac{ \cos U_2 \sin \lambda }{ \cos U_1 \sin U_2 - \sin U_1 \cos U_2 \cos \lambda } \right )$

$\alpha_2 = \arctan \left ( \frac{ \cos U_1 \sin \lambda }{ - \sin U_1 \cos U_2 + \cos U_1 \sin U_2 \cos \lambda } \right )$

Below is the Python implementation of the above equations. Note, I included a maximum number of iterations (maxIter=200).

#------------------------------------------------------------------------------+
#
#   Nathan A. Rooy
#   2016-SEP-30
#   Solve the inverse Vincenty's formulae
#   https://en.wikipedia.org/wiki/Vincenty%27s_formulae
#
#------------------------------------------------------------------------------+

#--- IMPORT DEPENDENCIES ------------------------------------------------------+

from __future__ import division
from math import atan
from math import atan2
from math import cos
from math import sin
from math import sqrt
from math import tan

#--- MAIN ---------------------------------------------------------------------+

class vincenty_inverse:
def __init__(self,coord1,coord2,maxIter=200,tol=10**-12):

#--- CONSTANTS ------------------------------------+

a=6378137.0                             # radius at equator in meters (WGS-84)
f=1/298.257223563                       # flattening of the ellipsoid (WGS-84)
b=(1-f)*a

phi_1,L_1,=coord1                       # (lat=L_?,lon=phi_?)
phi_2,L_2,=coord2

Lambda=L                                # set initial value of lambda to L

sin_u1=sin(u_1)
cos_u1=cos(u_1)
sin_u2=sin(u_2)
cos_u2=cos(u_2)

#--- BEGIN ITERATIONS -----------------------------+
self.iters=0
for i in range(0,maxIter):
self.iters+=1

cos_lambda=cos(Lambda)
sin_lambda=sin(Lambda)
sin_sigma=sqrt((cos_u2*sin(Lambda))**2+(cos_u1*sin_u2-sin_u1*cos_u2*cos_lambda)**2)
cos_sigma=sin_u1*sin_u2+cos_u1*cos_u2*cos_lambda
sigma=atan2(sin_sigma,cos_sigma)
sin_alpha=(cos_u1*cos_u2*sin_lambda)/sin_sigma
cos_sq_alpha=1-sin_alpha**2
cos2_sigma_m=cos_sigma-((2*sin_u1*sin_u2)/cos_sq_alpha)
C=(f/16)*cos_sq_alpha*(4+f*(4-3*cos_sq_alpha))
Lambda_prev=Lambda
Lambda=L+(1-C)*f*sin_alpha*(sigma+C*sin_sigma*(cos2_sigma_m+C*cos_sigma*(-1+2*cos2_sigma_m**2)))

# successful convergence
diff=abs(Lambda_prev-Lambda)
if diff<=tol:
break

u_sq=cos_sq_alpha*((a**2-b**2)/b**2)
A=1+(u_sq/16384)*(4096+u_sq*(-768+u_sq*(320-175*u_sq)))
B=(u_sq/1024)*(256+u_sq*(-128+u_sq*(74-47*u_sq)))
delta_sig=B*sin_sigma*(cos2_sigma_m+0.25*B*(cos_sigma*(-1+2*cos2_sigma_m**2)-(1/6)*B*cos2_sigma_m*(-3+4*sin_sigma**2)*(-3+4*cos2_sigma_m**2)))

self.m=b*A*(sigma-delta_sig)                 # output distance in meters
#self.km=self.meters/1000                    # output distance in kilometers
#self.mm=self.meters*1000                    # output distance in millimeters
#self.miles=self.meters*0.000621371          # output distance in miles
#self.n_miles=self.miles*(6080.20/5280)      # output distance in nautical miles
#self.ft=self.miles*5280                     # output distance in feet
#self.inches=self.feet*12                    # output distance in inches
#self.yards=self.feet/3                      # output distance in yards

if __name__ == "__vincenty_inverse__":
main()

Example usage:

>>> from vincenty import vincenty_inverse
>>> vincenty_inverse([39.152501,-84.412977],[39.152505,-84.412946]).m
2.7161912585815897

And that's it! Not too difficult right? The final code can be found on my GitHub (here).