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:
| Variable | Definition |
|---|---|
| \( 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 . Once through the first iteration, simply use this newly calculated value of and dump it back into the first equation and repeat the process. The end game here is to minimize the value of .
$$ \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 and the value of 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 radians
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
u_1=atan((1-f)*tan(radians(phi_1)))
u_2=atan((1-f)*tan(radians(phi_2)))
L=radians(L_2-L_1)
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).