Python Haversine Accuracy – Geocentric Radius vs Average Earth Radius

haversinepythonvincenty-formulae

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
    

    rlat1  = np.radians(lat1)
    rlong1 = np.radians(long1)
    rlat2  = np.radians(lat2)
    rlong2 = np.radians(long2)

    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

Best Answer

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):

http://clynchg3c.com/Technote/geodesy/radiigeo.pdf

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
    
    rlat1  = np.radians(lat1)
    rlong1 = np.radians(long1)
    rlat2  = np.radians(lat2)
    rlong2 = np.radians(long2)
    rBearing = np.radians(Bearing)

    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 
    RadiusofCurvature = 1/(np.cos(rBearing)**2/RM + np.sin(rBearing)**2/RN) #Radius of Curvature at azimuth


    arclength = np.arccos(np.sin(rlat1)*np.sin(rlat2) + np.cos(rlat1)*np.cos(rlat2)*np.cos(rlong2-rlong1)  )
        
    distance  = arclength * AverageR
    distance1 = arclength * GEOcentricRadius
    distance2 = arclength * RadiusofCurvature
    
    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
Related Question