How accurate is approximating the Earth as a sphere?

In short, the distance can be in error up to roughly 22km or 0.3%, depending on the points in question. That is:

• The error can be expressed in several natural, useful ways, such as (i) (residual) error, equal to the difference between the two calculated distances (in kilometers), and (ii) relative error, equal to the difference divided by the "correct" (ellipsoidal) value. To produce numbers convenient to work with, I multiply these ratios by 1000 to express the relative error in parts per thousand.

• The errors depend on the endpoints. Due to the rotational symmetry of the ellipsoid and sphere and their bilateral (north-south and east-west) symmetries, we may place one of the endpoints somewhere along the prime meridian (longitude 0) in the northern hemisphere (latitude between 0 and 90) and the other endpoint in the eastern hemisphere (longitude between 0 and 180).

To explore these dependencies, I have plotted the errors between endpoints at (lat,lon) = (mu,0) and (x,lambda) as a function of latitude x between -90 and 90 degrees. (All points are nominally at an ellipsoid height of zero.) In the figures, rows correspond to values of mu at {0, 22.5, 45, 67.5} degrees and columns to values of lambda at {0, 45, 90, 180} degrees. This gives us a good view of the spectrum of possibilities. As expected, their maximum sizes are approximately the flattening (around 1/300) times the major axis (around 6700 km), or about 22 km.

Errors

Relative errors

Contour plot

Another way to visualize the errors is to fix one endpoint and let the other vary, contouring the errors that arise. Here, for example, is a contour plot where the first endpoint is at 45 degrees north latitude, 0 degrees longitude. As before, error values are in kilometers and positive errors mean the spherical calculation is too large:

It might be easier to read when wrapped around the globe:

The red dot in the south of France shows the location of the first endpoint.

For the record, here is the Mathematica 8 code used for the calculations:

WGS84[x_, y_] := GeoDistance @@ (GeoPosition[Append[#, 0], "WGS84"] & /@ {x, y});
sphere[x_, y_] := GeoDistance @@
   (GeoPosition[{GeodesyData["WGS84", {"ReducedLatitude", #[[1]]}], #[[2]], 0}, "WGS84"] & /@ {x, y});

And one of the plotting commands:

With[{mu = 45}, ContourPlot[(sphere[{mu, 0}, {x, y}] - WGS84[{mu, 0}, {x, y}]) / 1000, 
                   {y, 0, 180}, {x, -90, 90}, ContourLabels -> True]]

I've explored this question recently. I think people want to know

1. what spherical radius should I use?
2. what is the resulting error?

A reasonable metric for the quality of the approximation is the maximum absolute relative error in the great-circle distance

err = |s_sphere - s_ellipsoid| / s_ellipsoid

with the maximum evaluated over all possible pairs of points.

If the flattening f is small, the spherical radius which minimizes err is very close to (a + b)/2 and the resulting error is about

err = 3*f/2 = 0.5% (for WGS84)

(evaluated with 10^6 randomly chosen pairs of points). It is sometimes suggested to use (2*a + b)/3 as the spherical radius. This results in a slightly larger error, err = 5*f/3 = 0.56% (for WGS84).

Geodesics whose length is most underestimated by the spherical approximation lie near a pole, e.g., (89.1,0) to (89.1,180). Geodesics whose length is most overestimated by the spherical approximation are meridional near the equator, e.g., (-0.1,0) to (0.1,0).

ADDENDUM: Here's another way of approaching this problem.

Select pairs of uniformly distributed points on the ellipsoid. Measure the ellipsoidal distance s and the distance on a unit sphere t. For any pair of points, s / t gives an equivalent spherical radius. Average this quantity over all the pairs of points and this gives a mean equivalent spherical radius. There's a question of exactly how the average should be done. However all the choices I tried

1. <s>/<t>
2. <s/t>
3. sqrt(<s^2>/<t^2>)
4. <s^2>/<s*t>
5. <s^2/t>/<s>

all came out within a few meters of the IUGG recommended mean radius, R₁ = (2 a + b) / 3. Thus, this value minimizes the RMS error in spherical distance calculations. (However it results in a slightly larger maximum relative error compared to (a + b) / 2; see above.) Given that R₁ is likely to be used for other purposes (area calculations and the like), there a good reason to stick with this choice for distance calculations.

The bottom line:

• For any kind of systematic work, where you can tolerate a 1% error in distance calculations, use a sphere of radius R₁. The maximum relative error is 0.56%. Use this value consistently when you approximate the earth with a sphere.
• It you need additional accuracy, solve the ellipsoidal geodesic problem.
• For back of the envelope calculations, use R₁ or 6400 km or 20000/pi km or a. These result in a maximum relative error of about 1%.

ANOTHER ADDENDUM: You can squeeze a little more accuracy out of the great circle distance by using μ = tan⁻¹((1 − f)^3/2 tanφ) (a poor man's rectifying latitude) as the latitude in the great circle calculation. This reduces the maximum relative error from 0.56% to 0.11% (using R₁ as the radius of the sphere). (It's not clear whether it's really worth taking this approach as opposed to computing the ellipsoidal geodesic distance directly.)