0

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

1 Answer 1

2

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.