1

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

enter image description hereenter image description here

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
    v = np.radians(v)

    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[0],data.shape[0]))
    for i in range(data.shape[0]):

        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. – mkennedy Oct 15 '19 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. – Bstampe Oct 15 '19 at 21:33
  • Which radius? The radius of curvature over a point on an ellipsoidal revolution surface varies with the azimuth. – Gabriel De Luca Oct 15 '19 at 22:39
  • The radius for the spherical approximation. – Bstampe Oct 15 '19 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. – Gabriel De Luca Oct 15 '19 at 22:44
2

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.


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

2

The x coordinate is:

3

The z coordinate is:

4

Finally, the radius of the sphere is:

5

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
latPradians = np.radians(latP)
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()

    v = np.radians(v)

    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[0],data.shape[0]))
    for i in range(data.shape[0]):

        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[0]
    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.

6


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.

7


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

8


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.

9


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

10

| improve this answer | |
2

There was a matrix algebra error in the Haversine formula. I updated the code in the question. I am getting much better agreement between Haversine and geodesic now:

enter image description here

On my actual dataset:

enter image description here enter image description here

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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