17

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?

Evan Carroll
  • 7,071
  • 2
  • 32
  • 58
  • 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 00:03

1 Answers1

23

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.

Evan Carroll
  • 7,071
  • 2
  • 32
  • 58