Potential bug in calculations in metersToLongitudeDegrees

[sruditsky]

Version info

Firebase: 3.1.0

GeoFire: 4.1.1

I think there is a small bug in calculations in metersToLongitudeDegrees (or to be more precise in the value of g_E2 this function is using):

These are the relevant snippets of code from geofire.js:

// g_E2 == (g_EARTH_EQ_RADIUS^2-g_EARTH_POL_RADIUS^2)/(g_EARTH_EQ_RADIUS^2)

var g_E2 = 0.00669447819799;

var metersToLongitudeDegrees = function(distance, latitude) {
  var radians = degreesToRadians(latitude);
  var num = Math.cos(radians)*g_EARTH_EQ_RADIUS*Math.PI/180;
  var denom = 1/Math.sqrt(1-g_E2*Math.sin(radians)*Math.sin(radians));
  var deltaDeg = num*denom;
  ...

As part of its calculation metersToLongitudeDegrees finds deltaDeg, which is the length of a one-degree arc of a parallel at the given latitude.

deltaDeg is calculated as:

cos(latitude) * g_EARTH_EQ_RADIUS*PI/180
/
sqrt(1 - g_E2 * sin(latitude) * sin(latitude))

R, the radius of the same parallel, should be equal to deltaDeg * 360 / (2 * PI), i.e.

R == 
(cos(latitude) * g_EARTH_EQ_RADIUS)
/
(sqrt(1 - g_E2 * sin(latitude) * sin(latitude)))

Would the Earth be a perfect sphere with radius of g_EARTH_EQ_RADIUS, then the radius at a particular parallel would be:

R_sphere == 
cos(latitude) * g_EARTH_EQ_RADIUS

However, because the polar radius of Earth is smaller than the equatorial it should be R < R_sphere.

The issue is that the calculation of R above (i.e. the one used by metersToLongitudeDegrees) renders R > R_sphere (for 1 > g_E2 > 0).

The following is my attempt to understand what is going on:

If each meridian is an ellipse with semi-axes equal to R_eq and R_pol then:

R² / R_eq² + R² * tan²(latitude) / R_pol² = 1

From here we get:

R² * (1 / R_eq² + tan²(latitude) / R_pol²) = 1
R = 1 / sqrt (1 / R_eq² + tan²(latitude) / R_pol²)
R = R_eq / sqrt (1 + R_eq² * tan²(latitude) / R_pol²)
R = R_eq / sqrt (1 + R_eq² * sin²(latitude) / (R_pol² * cos²(latitude)))
R = R_eq * cos(latitude) / sqrt (cos²(latitude) + R_eq² * sin²(latitude) / R_pol²)
R = R_eq * cos(latitude) / sqrt (1 - sin²(latitude) + R_eq² * sin²(latitude) / R_pol²)
R = R_eq * cos(latitude) / sqrt (1 - sin²(latitude) * (1 - R_eq² / R_pol²))
R = R_eq * cos(latitude) / sqrt (1 - sin²(latitude) * ( (R_pol² - R_eq²) / R_pol²))

g_E2_2 = (R_pol² - R_eq²) / R_pol²

R = cos(latitude) * R_eq / sqrt(1 - g_E2_2 * sin²(latitude))

This is very close to what metersToLongitudeDegrees is using, but in g_E2_2 the equatorial and polar radii are flipped compared to g_E2.

The calculation of g_E2_2 gives about -0.00673950125, which is a negative value and as such solves the issue described above.

[asciimike]

Any chance you've calculated the error ratio between using said values to see what % different values computed with one vs the other are? Are we talking off by millimeters, meters, kilometers...?

[sruditsky]

It depends where on Earth you make the measurement -- the closer to the poles the bigger the mistake.

The following shows dependency between the latitude and the difference in result between the two methods:

Latitude	Difference
10°	~0.02%
45°	~0.3%
60°	~0.5%
80°	~0.6%

---

[asciimike]

Ok, so it sounds like it's a pretty small difference. We use Geohashes under the hood, so behavior at the poles is pretty broken anyways (we also don't anticipate too many penguins or polar bears using Firebase ;)

At normal latitudes (~15-30º), 5 bits of precision is probably good enough to place you within 2 meters, and that's pretty much within the error here (assume ~0.1% = 1m/km, which would be within the geohash area of a 5 character hash). Geohash precision listed below:

# characters	km error
1	±2500
2	±630
3	±78
4	±20
5	±2.4
6	±0.61
7	±0.076
8	±0.019
9	±0.0024
10	±0.00060
11	±0.000074

---

[sruditsky]

I do agree that the common geofire query -- "give me everything in N-km radius" -- is not expected to be extremely precise.

However, then, the natural (while probably rhetorical) question is what on earth (no pun intended) the correction for the planet being a spheroid doing in the code in the first place?

I may be wrong, but the following seems to suggest that somebody really wanted to fix something:

// g_E2 == (g_EARTH_EQ_RADIUS^2-g_EARTH_POL_RADIUS^2)/(g_EARTH_EQ_RADIUS^2)
// The exact value is used here to avoid rounding errors
var g_E2 = 0.00669447819799;