15

MySQL says in the docs for ST_Distance_Sphere

Calculations use a spherical earth and a configurable radius. The optional radius argument should be given in meters. If omitted, the default radius is 6,370,986 meters. If the radius argument is present but not positive, an ER_WRONG_ARGUMENTS error occurs.

PostGIS says in the docs of ST_Distance_Sphere, (though the docs aren't accurate anymore)

Uses a spherical earth and radius of 6370986 meters.

Where did they get the default 6,370,986 meters from? WGS84 says major-axis radius is 6,378,137.0 m. PostGIS which now uses a Average Radius essentially uses 6371008.

Looking at the code

#define WGS84_MAJOR_AXIS 6378137.0
#define WGS84_INVERSE_FLATTENING 298.257223563
#define WGS84_MINOR_AXIS (WGS84_MAJOR_AXIS - WGS84_MAJOR_AXIS / WGS84_INVERSE_FLATTENING)
#define WGS84_RADIUS ((2.0 * WGS84_MAJOR_AXIS + WGS84_MINOR_AXIS ) / 3.0)

that means

-- SELECT 6378137.0 - 6378137.0 / 298.257223563;
WGS84_MINOR_AXIS = 6356752.314245179498
-- SELECT ( 2.0 * 6378137.0 + ( 6378137.0 - 6378137.0 / 298.257223563) ) / 3.0;
WGS84_RADIUS = 6371008.771415059833

Newer versions are much less efficient, more complex, and use Pro4j but they seem to do the same thing.

Still where does 6370986 come from?

  • 1
    It represents the mean earth radius, which should be (2*minorAxis+majorAxis)/3 ... though that value for WGS84 is still a few meters larger (6,371,008.771) – JGH Nov 20 '17 at 18:01
  • yep, that's the question why the discrepancy. – Evan Carroll Nov 20 '17 at 18:08
  • 2
    Some developer looked it up on the net? The postgis source may throw some light on it – Ian Turton Nov 20 '17 at 18:32
  • 2
    @IanTurton Most bugs can be reduced to "some developer did something and the source may throw light on it." I intended to do the work, figuring that would be what it takes if no one knew the story. See the answer below. – Evan Carroll Nov 20 '17 at 18:49
  • 1
    Maybe there was a typo and they meant 6370996...that's very close to Clarke 1866's authalic radius. – mkennedy Nov 21 '17 at 0:03
21

Ok, this is hilarriuusss. I tracked this down. In an old copy of lwgeom/lwgeom_spheroid.c in PostGIS 1.0.0rc4 you can see this,

/*
 * This algorithm was taken from the geo_distance function of the 
 * earthdistance package contributed by Bruno Wolff III.
 * It was altered to accept GEOMETRY objects and return results in
 * meters.
 */
PG_FUNCTION_INFO_V1(LWGEOM_distance_sphere);
Datum LWGEOM_distance_sphere(PG_FUNCTION_ARGS)
{
        const double EARTH_RADIUS = 6370986.884258304;

Moving on to the docs of earthdistance, you'll find this:

Note that unlike the cube-based part of the module, units are hardwired here: changing the earth() function will not affect the results of this operator.

And that hard-wired number: EARTH_RADIUS can be seen here

/* Earth's radius is in statute miles. */
static const double EARTH_RADIUS = 3958.747716;

So you can do a simple.

EARTH_RADIUS * MILES_TO_METERS = EARTH_RADIUS_IN_METERS
 3958.747716 * 1609.344        = 6370986.884258304

And you have your 6370986.884258304. Of course, just truncate that and store it in a long because why not.

So in essence, the radius in MySQL was lifted from a lazy-copy-job from PostGIS that converted a radius in miles to meters from an obscure constant from a random 20-year old PostgreSQL module.

earth_distance is a pre-PostGIS module by Bruce Momjian. I hereby proclaim 6370986 the Bmomjian Constant: a good nuff' approximation of Earth in meters to satisfy MySQL. Though maybe not for long.

  • 2
    But where did that very precise figure 3958.747716 come from, then? The closest I can find is 3958.74795, which is the number of U.S. survey miles in 6371 kilometers, but that still leaves about 37 cm unaccounted for ,,, – Henning Makholm Nov 20 '17 at 23:06
  • 1
    @HenningMakholm keep fighting the good fight no idea. ;) – Evan Carroll Nov 20 '17 at 23:07
  • 2
    Very nice find! – Paul Ramsey Nov 21 '17 at 1:29

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.