# When calcuating distance between points on Earth why are my Haversine vs. Geodesic calculations wildy diverging?

I am getting wildly diverging distances using two approximations to calculate distance between points on Earth's surface. I am using the Haversine (vectorized) approximation (spherical earth) and the more precise (presumably) geopy.distance.geodesic(ellipsoidal earth) . See here for discussion of Haversine n Geodesic

As you can see I am off by five percent as the distances between points becomes large. Is this divergence due to rounding error in Haversine? Do I indeed trust the Geodesic? Here is code:

``````import numpy as np
lat = np.linspace(35,45,100)
lon = np.linspace(-120,-110,100)

data = pd.DataFrame({'Latitude':lat,'Longitude':lon})

def Haversine(v):
"""
distance between two lat,lon coordinates
using the Haversine formula. Assumes one
radius. r = 3,950 to 3,963 mi
"""
from timeit import default_timer as timer
start = timer()
R = 3958 # radius at 40 deg 750 m elev

dlat = v[:, 0, np.newaxis] - v[:, 0]
dlon = v[:, 1, np.newaxis] - v[:, 1]
c = np.cos(v[:,0,None]
a = np.sin(dlat / 2.0) ** 2 + c @ c.T * np.sin(dlon / 2.0) ** 2

c = 2 * np.arcsin(np.sqrt(a))
result = R * c
print(round((timer() - start),3))
return result

def slowdistancematrix(data):

from geopy.distance import geodesic
distance = np.zeros((data.shape,data.shape))
for i in range(data.shape):

lat_lon_i = data.Latitude.iloc[i],data.Longitude.iloc[i]

for j in range(i):

lat_lon_j = data.Latitude.iloc[j],data.Longitude.iloc[j]

distance[i,j] = geodesic(lat_lon_i, lat_lon_j).miles
distance[j,i] = distance[i,j] # make use of symmetry

return distance

distanceG = slowdistancematrix(data)
distanceH = Haversine(data.values)

plt.scatter(distanceH.ravel(),distanceG.ravel()/distanceH.ravel(),s=.5)
plt.ylabel('Geodesic/Haversine')
plt.xlabel('Haversine distance (miles)')
plt.title('all points in distance matrix')

``````

I would rather use the vectorized version becuase it is fast. However,the 5% is too big for me to be comfortable with it. Supposedly Haversine is only suppose to be off by .5%.

Update Solved. Matrix algebra error.

• I suspect that part of the problem is that you're using a sphere that's in-between the ellipsoid being used. If you used the major auxiliary sphere (R = semimajor axis), routes below say mid-latitudes would match better, in exchange for more divergence closer to the poles. Oct 15, 2019 at 20:30
• I am using the radius at the mid center of the Latitudes I want to calcuate. If I change this to either the max or the min it doesn't change the divergence but it does move up or down the mean. Oct 15, 2019 at 21:33
• Which radius? The radius of curvature over a point on an ellipsoidal revolution surface varies with the azimuth. Oct 15, 2019 at 22:39
• The radius for the spherical approximation. Oct 15, 2019 at 22:40
• If you approximate the sphere according to the radius of curvature of the normal meridian section at that midpoint, or according to the radius of curvature of the normal section to the normal meridian, at the same point, or according to the radius of curvature at that point of the normal section in a certain azimuth, all different radius, you will be getting different distances, because you are approaching different spheres supported on the same point. Oct 15, 2019 at 22:44

I appreciate your code because it contributes to geodesy, by allowing to make some interesting analyzes in a practical way, regarding spherical approximations for the calculation of length of geodesic lines.

The definition of the spherical radius that approximates the ellipsoid at one point, is a definition that is, at least, ambiguous. What I propose is to precisely define a sphere and it radius R, that can be correctly calculated deterministically.

When we perform geodetic calculations, we have no way of modeling the topographic form of the Earth, then we make a first approximation: the ellipsoid of revolution, on which we can perform geodetic calculations (and from where we can project plane representations).

Differential calculus on the surface of the ellipsoid is not easy, so we can look for one more approximation: the spherical one.
But which sphere are we talking about?
Usually, when we calculate a single geodetic line, we choose one or more special spheres, based on radii of curvature, for that approximation.
When we need to approximate a sphere that allows us to calculate the geodesic from any point to any other, we also choose a particular sphere, whose definition does not involve any point at any latitude.

But in this case we want to define another one, so let's define it like this:

The sphere that has the same center as the ellipsoid, and contains a point P.
P is any point of the ellipsoid, which has latitude lat.
R is the radius of that sphere.

Will that approximation be good?
It depends on our definition of good approximation. But let's experiment.

Note 1:
I assure you that the sphere of radius R does not match, strictly speaking, meridian arcs with arcs of circumference at equal latitudes of P.
Also, I assure you that it is not necessarily tangential to the ellipsoid at any point.
I can only assure you the definition: that sphere contains P.

Note 2:
Without going into details, considering the elevation of P can only bring us problems. If you want to use a larger radius, decrease the latitude of P. The diagram represents a fragment of the meridian section that contains the point P, which when rotating generates an ellipsoid of revolution; and the associated sphere, both centered on O.
a and b are the ellipsoid parameters.
n is normal to the ellipsoid that passes through P.
P has zero longitude and phi latitude.
x and z are geocentric coordinates of P.
R is the radius of the O-centered sphere that contains P.

The square of the eccentricity of the ellipsoid is:

The x coordinate is:

The z coordinate is:

Finally, the radius of the sphere is:

Thus, for the WGS-84 ellipsoid and a point at 40º latitude, R = 6369344.86 meters = 3957.73725 miles.
The sphere of radius R = 3958 miles passes through a point P at latitude 38.84323º. We do not know where you took that value. No problem because it can be considered a new approximation.
The latitude = 35.44539 of P corresponds to the volumetric radius R = 3958.7564 miles.

So, let me to modify your code to do some tests:

``````import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from timeit import default_timer as timer

fromlat = -90
tolat = 90
fromlon = 0
tolon = 180
n = 200 # number of points
latP = 35.44539 # latitude of P (35.44539 for the volumetric radius)

lat = np.linspace(fromlat,tolat,n)
lon = np.linspace(fromlon,tolon,n)
data = pd.DataFrame({'Latitude':lat,'Longitude':lon})

aellips = 6378137 # a parameter of the ellipsoid (6378137 for WGS-84)
bellips = 6356752.3 # b parameter of the ellipsoid (6356752.3 for WGS-84)
e2 = 1 - bellips ** 2 / aellips ** 2
xP = (aellips * np.cos(latPradians)) / (np.sqrt(1 - e2 * np.sin(latPradians) ** 2))
zP = (aellips * (1 - e2) * np.sin(latPradians)) / (np.sqrt(1 - e2 * np.sin(latPradians) ** 2))
Rmeters = np.sqrt(xP ** 2 + zP ** 2)
R = Rmeters / 1609.34 # Radius of the sphere (miles)

print("\nFrom lat =", fromlat, "to lat =", tolat, "(decimal degrees),")
print("from lon =", fromlon, "to lon =", tolon, "(decimal degrees),")
print("with", n, "points.")
print("\nLatitude of P =", latP, "(decimal degrees),")
print("Radius of the sphere that cointains P =", round(R,5), "(miles).")

def Haversine(v,R):
"""
distance between two lat,lon coordinates
using the Haversine formula.
"""
start = timer()

dlat = v[:, 0, np.newaxis] - v[:, 0]
dlon = v[:, 1, np.newaxis] - v[:, 1]
c = np.cos(v[:,0,None])
a = np.sin(dlat / 2.0) ** 2 + c @ c.T * np.sin(dlon / 2.0) ** 2
c = 2 * np.arcsin(np.sqrt(a))
distance = R * c

print("\nHaversine =", np.count_nonzero(distance), "distances calculated in",
round((timer() - start),3), "seconds.")
return distance

def slowdistancematrix(data):
from geopy.distance import geodesic
start = timer()
distance = np.zeros((data.shape,data.shape))
for i in range(data.shape):

lat_lon_i = data.Latitude.iloc[i],data.Longitude.iloc[i]

for j in range(i):

lat_lon_j = data.Latitude.iloc[j],data.Longitude.iloc[j]

distance[i,j] = geodesic(lat_lon_i, lat_lon_j).miles
distance[j,i] = distance[i,j] # make use of symmetry

print("Geodesic =", np.count_nonzero(distance), "distances calculated in",
round((timer() - start),3), "seconds.")
return distance

def skip_diag_strided(A):
m = A.shape
strided = np.lib.stride_tricks.as_strided
s0,s1 = A.strides
return strided(A.ravel()[1:], shape=(m-1,m), strides=(s0+s1,s1)).reshape(m,-1)

distanceH = Haversine(data.values,R)
distanceG = slowdistancematrix(data)

distH = skip_diag_strided(distanceH) # avoid zero diagonal
distG = skip_diag_strided(distanceG) # avoid zero diagonal
d = distG - distH # absolute difference
pd = 100 * (distG-distH) / distG # porcentual difference

plt.figure(1)
plt.scatter(distG.ravel(),d.ravel(),s=.5)
plt.ylabel('Geodesic - Haversine distance [miles]')
plt.xlabel('Geodesic distance [miles]')
plt.title('Absolute Difference')
plt.figure(2)
plt.scatter(distG.ravel(),pd.ravel(),s=.5)
plt.ylabel('100*(Geodesic-Haversine)/Geodesic [%]')
plt.xlabel('Geodesic distance [miles]')
plt.title('Porcentual Difference')
``````

Ok, first things first:

Q: Is it true that Haversine's formula returns a maximum porcentual difference of 0.5% between distances from any to any point on Earth using the volumetric radius?
A: Yes, it seems to be true.

Q: Is the approximation of the radius of 3958 miles good for calculating the distances between the question points?
A: It seems not to be bad.

Q: It could be better?
A: I think so.

Q: Excuse me, I liked the program. Can we verify that the geodesic paths on the equator are equivalent to those for a sphere of radius a, and in turn observe the error produced by the Haversine formula for antipodal points?
A: Yes, sure.

Q: Finally, why doesn't it work the same way for trips to and from points at a parallel with latitude other than zero?  