I am required to calculate the distance between two points.

To reduce computational complexity for an embedded system, rather than use the more accurate Vincenty that uses the WSG84 ellipsoid, I have decided to use haversine to calculate the angular distance over a great circle arc then multiply that with the Earth radius.

I was comparing the accuracy between haversine vs Vincenty. When I use a geocentric radius, I get a worse accuracy compared with using the average radius.

Does anyone know why this is? I would expect smaller errors when using a more accurate radius.

``````import numpy as np
from pymap3d import vincenty

def simpleHav(lat1, long1, lat2, long2):
"""
Given 2 positions provide the distance (shortest distance) great circle arc.
Inputs in degrees lat long
Output is a length in metres
"""

AverageR = 6371000  # Earth Radius

r1 = 6378137
r2 = 6356752

R = np.sqrt(((r1**2*np.cos(rlat1))**2 + (r2**2*np.sin(rlat1))**2)/((r1*np.cos(rlat1))**2 + (r2*np.sin(rlat1))**2))

arclength = np.arccos(np.sin(rlat1)*np.sin(rlat2) + np.cos(rlat1)*np.cos(rlat2)*np.cos(rlong2-rlong1)  )
distance  = arclength * R
distance1 = arclength * AverageR

return distance, distance1

VincentyRange = vincenty.vdist(50,10, 51, 11)[0]

Haversine = simpleHav(50, 10,  51, 11)[0]
Haversine1 = simpleHav(50, 10, 51, 11)[1]

print("Error using GEOcentric Radius = " + str(VincentyRange - Haversine))
print("Error using Average Radius = " + str(VincentyRange - Haversine1))
``````

In this case where I have provided a start position of:

``````lat = 50
long = 10
``````

and an end position of

``````lat = 51
long = 11
``````

I get the following errors:

``````Error using GEOcentric Radius = 266.5363466117997

Error using Average Radius = 155.48858379435842
``````

Ok I figured it out after some research. Instead of using the Geocentric radius I instead calculate the radius Use the “Radius of Curvature formula at azimuth α” formula from the paper (page 5):

Here is the code updated....

``````import numpy as np
from pymap3d import vincenty

def simpleHav(lat1, long1, lat2, long2, Bearing):
"""
Given 2 positions provide the distance (shortest distance) great circle arc.
Inputs in degrees lat long
Output is a length in metres
"""

AverageR = 6371000  # Earth Radius

a = 6378137 #Semi Major Axis a
b = 6356752 #Semi Minor Axis b
e = np.sqrt(1-(b**2/a**2)) #eccentricity

GEOcentricRadius = np.sqrt(((a**2*np.cos(rlat1))**2 + (b**2*np.sin(rlat1))**2)/((a*np.cos(rlat1))**2 + (b*np.sin(rlat1))**2))

RN = a/np.sqrt(1-e**2*np.sin(rlat1)**2)         #Radius of Curvature in Prime Vertical, terminated by minor axis
RM = RN * ((1-e**2)/(1-e**2*np.sin(rlat1)**2))  #Radius of Curvature: in Meridian

arclength = np.arccos(np.sin(rlat1)*np.sin(rlat2) + np.cos(rlat1)*np.cos(rlat2)*np.cos(rlong2-rlong1)  )

distance  = arclength * AverageR

return distance, distance1, distance2

VincentyRange = vincenty.vdist(50, 10, 51, 11)[0]
Haversine          = simpleHav(50, 10, 51, 11,32.07)

print("Error using GEOcentric Radius = " + str(VincentyRange - Haversine[1]))
print("Error using Average Radius = " + str(VincentyRange - Haversine[0]))
print("Error using Radius of curvature = " + str(VincentyRange - Haversine[2]))
``````

The improvement is significant... If I provide the same inputs as before.

Start position:

``````lat = 50
long = 10
``````

End position:

``````lat = 51
long = 11

Error using GEOcentric Radius = 266.5363466117997
Error using Average Radius = 155.48858379435842
Error using Radius of curvature = 11.756586791569134
``````