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 $$\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).

Example usage:

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