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
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.
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?
A: Short answer, curvature radii.
