Geo: Fix geohashEstimateStepsByRadius() step underestimation.

The returned step was in some case not enough towards normal
coordinates (for example when our search position was was already near the
margin of the central area, and we had to match, using the east or west
neighbor, a very far point). Example:

    geoadd points 67.575457940146066 -62.001317572780565 far
    geoadd points 66.685439060295664 -58.925040587282297 center
    georadius points 66.685439060295664 -58.925040587282297 200 km

In the above case the code failed to find a match (happens at smaller
latitudes too) even if far and center are at less than 200km.

Another fix introduced by this commit is a progressively larger area
towards the poles, since meridians are a lot less far away, so we need
to compensate for this.

The current implementation works comparably to the Tcl brute-force
stress tester implemented in the fuzzy test in the geo.tcl unit for
latitudes between -70 and 70, and is pretty accurate over +/-80 too,
with sporadic false negatives.

A more mathematically clean implementation is possible by computing the
meridian distance at the specified latitude and computing the step
according to it.

2 files changed, 18 insertions, 9 deletions

diff --git a/deps/geohash-int/geohash_helper.c b/deps/geohash-int/geohash_helper.c
index 7271c7b31..b6a00b7b5 100644
--- a/deps/geohash-int/geohash_helper.c
+++ b/deps/geohash-int/geohash_helper.c
@@ -55,16 +55,24 @@ static inline double rad_deg(double ang) { return ang / D_R; }
 
 /* You must *ONLY* estimate steps when you are encoding.
  * If you are decoding, always decode to GEO_STEP_MAX (26). */
-uint8_t geohashEstimateStepsByRadius(double range_meters) {
-    uint8_t step = 1;
-    while (range_meters > 0 && range_meters < MERCATOR_MAX) {
+uint8_t geohashEstimateStepsByRadius(double range_meters, double lat) {
+    if (range_meters == 0) return 26;
+    int step = 1;
+    while (range_meters < MERCATOR_MAX) {
         range_meters *= 2;
         step++;
     }
-    step--;
-    if (!step)
-        step = 26; /* if range = 0, give max resolution */
-    return step > 26 ? 26 : step;
+    step -= 2; /* Make sure range is included in the worst case. */
+    /* Wider range torwards the poles... Note: it is possible to do better
+     * than this approximation by computing the distance between meridians
+     * at this latitude, but this does the trick for now. */
+    if (lat > 67 || lat < -67) step--;
+    if (lat > 80 || lat < -80) step--;
+
+    /* Frame to valid range. */
+    if (step < 1) step = 1;
+    if (step > 26) step = 25;
+    return step;
 }
 
 int geohashBitsComparator(const GeoHashBits *a, const GeoHashBits *b) {
@@ -114,11 +122,12 @@ GeoHashRadius geohashGetAreasByRadius(double latitude, double longitude, double
     max_lat = bounds[2];
     max_lon = bounds[3];
 
-    steps = geohashEstimateStepsByRadius(radius_meters);
+    steps = geohashEstimateStepsByRadius(radius_meters,latitude);
 
     geohashGetCoordRange(&lat_range, &long_range);
     geohashEncode(&lat_range, &long_range, latitude, longitude, steps, &hash);
     geohashNeighbors(&hash, &neighbors);
+    geohashGetCoordRange(&lat_range, &long_range);
     geohashDecode(lat_range, long_range, hash, &area);
 
     if (area.latitude.min < min_lat) {
diff --git a/deps/geohash-int/geohash_helper.h b/deps/geohash-int/geohash_helper.h
index c10a02c6a..9354a3238 100644
--- a/deps/geohash-int/geohash_helper.h
+++ b/deps/geohash-int/geohash_helper.h
@@ -48,7 +48,7 @@ typedef struct {
 } GeoHashRadius;
 
 int GeoHashBitsComparator(const GeoHashBits *a, const GeoHashBits *b);
-uint8_t geohashEstimateStepsByRadius(double range_meters);
+uint8_t geohashEstimateStepsByRadius(double range_meters, double lat);
 int geohashBoundingBox(double latitude, double longitude, double radius_meters,
                         double *bounds);
 GeoHashRadius geohashGetAreasByRadius(double latitude,