77

What level of error do I encounter when approximating the earth as a sphere? Specifically, when dealing with the location of points and, for example, the great circle distances between them.

Are there any studies on the average and worst case error compared to an ellipsoid? I'm wondering how much accuracy I'd be sacrificing if I go with a sphere for the sake of easier calculations.

My particular scenario involves directly mapping WGS84 coordinates as if they were coordinates on a perfect sphere (with the mean radius defined by the IUGG) without any transformation.

8
  • Are you specifically interested in a spherical model or are you interested in ellipsoid models? I imagine that the amount of error would vary greatly between a sphere and an ellipse.
    – Jay Laura
    May 15, 2012 at 18:30
  • 2
    A related analysis appears in this reply. To obtain an answer to your question, though, you need to specify how the earth is approximated as a sphere. Many approximations are in use. They all are tantamount to giving functions f' = u(f,l) and l' = v(f,l) where (f,l) are geographical coordinates of the sphere and (f',l') are geographical coordinates of the ellipsoid. See Section 1.7 ("Transformation...of the ellipsoid of revolution onto the surface of a sphere") in Bugayevskiy & Snyder, Map Projections, A Reference Manual. Taylor & Francis [1995].
    – whuber
    May 15, 2012 at 18:44
  • This is akin to the early debate over the Google/Bing EPSG 900913 projection (which uses WGS84 coordinates but projects tham as if they were on a sphere) and the errors probably account for EPSG initially rejecting the projection until giving in to pressure from developers. Without want to overly distract you, following up on some of this debate can add some additional breadth to the information in excellent link provided by whuber. May 15, 2012 at 19:40
  • @Jzl5325: Yup, I meant a strict sphere and not ellipsoid, edited the question to provide a bit more context too.
    – Jeff B
    May 15, 2012 at 19:49
  • 1
    I think you should read this: en.wikipedia.org/wiki/Haversine_formula
    – longtsing
    Apr 17, 2019 at 2:51

2 Answers 2

105

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

Residual errors

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

Contour plot

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

Globe plot

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]]
5
  • 31
    What a badass answer @whuber Dec 10, 2012 at 19:44
  • 1
    ok, but which spherical radius minimize the error ? for wgs84
    – yota
    Nov 23, 2021 at 8:34
  • @yota Minimize which error? The error in one particular calculation or some kind of average error?
    – whuber
    Nov 23, 2021 at 16:11
  • @whuber well, you are right... let says, which minimize the average error for the distance computed between two random points :) I found that the IUGG often uses R1 = (2*a+b)/3 (where a and b are the WGS84 semi axis)
    – yota
    Nov 24, 2021 at 16:14
  • @yota That's a useful rule of thumb. However, in many applications the points are not uniformly and independently distributed over the surface of the earth. If you assume a distribution appropriate for your application (for instance, the points might all be located within a given country) you can compute a spherical radius that provides a better approximation and less average error.
    – whuber
    Nov 24, 2021 at 16:29
25

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, R1 = (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 R1 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 R1. 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 R1 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((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 R1 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.)

Your Answer

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

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