Add GEOSEARCH / GEOSEARCHSTORE commands (#8094)

Add commands to query geospatial data with bounding box.

Two new commands that replace the existing 4 GEORADIUS* commands.

GEOSEARCH key [FROMMEMBER member] [FROMLOC long lat] [BYRADIUS radius
unit] [BYBOX width height unit] [WITHCORD] [WITHDIST] [WITHASH] [COUNT
count] [ASC|DESC]

GEOSEARCHSTORE dest_key src_key [FROMMEMBER member] [FROMLOC long lat]
[BYRADIUS radius unit] [BYBOX width height unit] [WITHCORD] [WITHDIST]
[WITHASH] [COUNT count] [ASC|DESC] [STOREDIST]

- Add two types of CIRCULAR_TYPE and RECTANGLE_TYPE to achieve different searches
- Judge whether the point is within the rectangle, refer to:
geohashGetDistanceIfInRectangle

7 files changed, 317 insertions, 90 deletions

diff --git a/src/geo.c b/src/geo.c
index 7594bb2e9..13cbb1c50 100644
--- a/src/geo.c
+++ b/src/geo.c
@@ -144,31 +144,57 @@ double extractUnitOrReply(client *c, robj *unit) {
 
 /* Input Argument Helper.
  * Extract the distance from the specified two arguments starting at 'argv'
- * that should be in the form: <number> <unit>, and return the distance in the
- * specified unit on success. *conversions is populated with the coefficient
- * to use in order to convert meters to the unit.
- *
- * On error a value less than zero is returned. */
-double extractDistanceOrReply(client *c, robj **argv,
-                                     double *conversion) {
+ * that should be in the form: <number> <unit>, and return C_OK or C_ERR means success or failure
+ * *conversions is populated with the coefficient to use in order to convert meters to the unit.*/
+int extractDistanceOrReply(client *c, robj **argv,
+                              double *conversion, double *radius) {
     double distance;
     if (getDoubleFromObjectOrReply(c, argv[0], &distance,
                                    "need numeric radius") != C_OK) {
-        return -1;
+        return C_ERR;
     }
 
     if (distance < 0) {
         addReplyError(c,"radius cannot be negative");
-        return -1;
+        return C_ERR;
     }
+    if (radius) *radius = distance;
 
     double to_meters = extractUnitOrReply(c,argv[1]);
     if (to_meters < 0) {
-        return -1;
+        return C_ERR;
+    }
+
+    if (conversion) *conversion = to_meters;
+    return C_OK;
+}
+
+/* Input Argument Helper.
+ * Extract height and width from the specified three arguments starting at 'argv'
+ * that should be in the form: <number> <number> <unit>, and return C_OK or C_ERR means success or failure
+ * *conversions is populated with the coefficient to use in order to convert meters to the unit.*/
+int extractBoxOrReply(client *c, robj **argv, double *conversion,
+                         double *height, double *width) {
+    double h, w;
+    if ((getDoubleFromObjectOrReply(c, argv[0], &h, "need numeric height") != C_OK) ||
+        (getDoubleFromObjectOrReply(c, argv[1], &w, "need numeric width") != C_OK)) {
+        return C_ERR;
+    }
+
+    if (h < 0 || w < 0) {
+        addReplyError(c, "height or width cannot be negative");
+        return C_ERR;
+    }
+    if (height) *height = h;
+    if (width) *width = w;
+
+    double to_meters = extractUnitOrReply(c,argv[2]);
+    if (to_meters < 0) {
+        return C_ERR;
     }
 
     if (conversion) *conversion = to_meters;
-    return distance * to_meters;
+    return C_OK;
 }
 
 /* The default addReplyDouble has too much accuracy.  We use this
@@ -183,21 +209,22 @@ void addReplyDoubleDistance(client *c, double d) {
 }
 
 /* Helper function for geoGetPointsInRange(): given a sorted set score
- * representing a point, and another point (the center of our search) and
- * a radius, appends this entry as a geoPoint into the specified geoArray
- * only if the point is within the search area.
+ * representing a point, and a GeoShape, appends this entry as a geoPoint
+ * into the specified geoArray only if the point is within the search area.
  *
  * returns C_OK if the point is included, or REIDS_ERR if it is outside. */
-int geoAppendIfWithinRadius(geoArray *ga, double lon, double lat, double radius, double score, sds member) {
-    double distance, xy[2];
+int geoAppendIfWithinShape(geoArray *ga, GeoShape *shape, double score, sds member) {
+    double distance = 0, xy[2];
 
     if (!decodeGeohash(score,xy)) return C_ERR; /* Can't decode. */
     /* Note that geohashGetDistanceIfInRadiusWGS84() takes arguments in
      * reverse order: longitude first, latitude later. */
-    if (!geohashGetDistanceIfInRadiusWGS84(lon,lat, xy[0], xy[1],
-                                           radius, &distance))
-    {
-        return C_ERR;
+    if (shape->type == CIRCULAR_TYPE) {
+        if (!geohashGetDistanceIfInRadiusWGS84(shape->xy[0], shape->xy[1], xy[0], xy[1],
+                                               shape->t.radius*shape->conversion, &distance)) return C_ERR;
+    } else if (shape->type == RECTANGLE_TYPE) {
+        if (!geohashGetDistanceIfInRectangle(shape->bounds, shape->xy[0], shape->xy[1],
+                                             xy[0], xy[1], &distance)) return C_ERR;
     }
 
     /* Append the new element. */
@@ -222,7 +249,7 @@ int geoAppendIfWithinRadius(geoArray *ga, double lon, double lat, double radius,
  * using multiple queries to the sorted set, that we later need to sort
  * via qsort. Similarly we need to be able to reject points outside the search
  * radius area ASAP in order to allocate and process more points than needed. */
-int geoGetPointsInRange(robj *zobj, double min, double max, double lon, double lat, double radius, geoArray *ga) {
+int geoGetPointsInRange(robj *zobj, double min, double max, GeoShape *shape, geoArray *ga) {
     /* minex 0 = include min in range; maxex 1 = exclude max in range */
     /* That's: min <= val < max */
     zrangespec range = { .min = min, .max = max, .minex = 0, .maxex = 1 };
@@ -254,7 +281,7 @@ int geoGetPointsInRange(robj *zobj, double min, double max, double lon, double l
             ziplistGet(eptr, &vstr, &vlen, &vlong);
             member = (vstr == NULL) ? sdsfromlonglong(vlong) :
                                       sdsnewlen(vstr,vlen);
-            if (geoAppendIfWithinRadius(ga,lon,lat,radius,score,member)
+            if (geoAppendIfWithinShape(ga,shape,score,member)
                 == C_ERR) sdsfree(member);
             zzlNext(zl, &eptr, &sptr);
         }
@@ -275,7 +302,7 @@ int geoGetPointsInRange(robj *zobj, double min, double max, double lon, double l
                 break;
 
             ele = sdsdup(ele);
-            if (geoAppendIfWithinRadius(ga,lon,lat,radius,ln->score,ele)
+            if (geoAppendIfWithinShape(ga,shape,ln->score,ele)
                 == C_ERR) sdsfree(ele);
             ln = ln->level[0].forward;
         }
@@ -315,15 +342,15 @@ void scoresOfGeoHashBox(GeoHashBits hash, GeoHashFix52Bits *min, GeoHashFix52Bit
 /* Obtain all members between the min/max of this geohash bounding box.
  * Populate a geoArray of GeoPoints by calling geoGetPointsInRange().
  * Return the number of points added to the array. */
-int membersOfGeoHashBox(robj *zobj, GeoHashBits hash, geoArray *ga, double lon, double lat, double radius) {
+int membersOfGeoHashBox(robj *zobj, GeoHashBits hash, geoArray *ga, GeoShape *shape) {
     GeoHashFix52Bits min, max;
 
     scoresOfGeoHashBox(hash,&min,&max);
-    return geoGetPointsInRange(zobj, min, max, lon, lat, radius, ga);
+    return geoGetPointsInRange(zobj, min, max, shape, ga);
 }
 
 /* Search all eight neighbors + self geohash box */
-int membersOfAllNeighbors(robj *zobj, GeoHashRadius n, double lon, double lat, double radius, geoArray *ga) {
+int membersOfAllNeighbors(robj *zobj, GeoHashRadius n, GeoShape *shape, geoArray *ga) {
     GeoHashBits neighbors[9];
     unsigned int i, count = 0, last_processed = 0;
     int debugmsg = 0;
@@ -374,7 +401,7 @@ int membersOfAllNeighbors(robj *zobj, GeoHashRadius n, double lon, double lat, d
                 D("Skipping processing of %d, same as previous\n",i);
             continue;
         }
-        count += membersOfGeoHashBox(zobj, neighbors[i], ga, lon, lat, radius);
+        count += membersOfGeoHashBox(zobj, neighbors[i], ga, shape);
         last_processed = i;
     }
     return count;
@@ -454,51 +481,60 @@ void geoaddCommand(client *c) {
 
 #define RADIUS_COORDS (1<<0)    /* Search around coordinates. */
 #define RADIUS_MEMBER (1<<1)    /* Search around member. */
-#define RADIUS_NOSTORE (1<<2)   /* Do not acceot STORE/STOREDIST option. */
+#define RADIUS_NOSTORE (1<<2)   /* Do not accept STORE/STOREDIST option. */
+#define GEOSEARCH (1<<3)        /* GEOSEARCH command variant (different arguments supported) */
+#define GEOSEARCHSTORE (1<<4)   /* GEOSEARCHSTORE just accept STOREDIST option */
 
 /* GEORADIUS key x y radius unit [WITHDIST] [WITHHASH] [WITHCOORD] [ASC|DESC]
  *                               [COUNT count] [STORE key] [STOREDIST key]
- * GEORADIUSBYMEMBER key member radius unit ... options ... */
-void georadiusGeneric(client *c, int flags) {
-    robj *key = c->argv[1];
+ * GEORADIUSBYMEMBER key member radius unit ... options ...
+ * GEOSEARCH key [FROMMEMBER member] [FORMLOG long lat] [BYRADIUS radius unit]
+ *               [BYBOX width height unit] [WITHCORD] [WITHDIST] [WITHASH] [COUNT count] [ASC|DESC]
+ * GEOSEARCHSTORE dest_key src_key [FROMMEMBER member] [FORMLOG long lat] [BYRADIUS radius unit]
+ *               [BYBOX width height unit] [WITHCORD] [WITHDIST] [WITHASH] [COUNT count] [ASC|DESC] [STOREDIST]
+ *  */
+void georadiusGeneric(client *c, int srcKeyIndex, int flags) {
     robj *storekey = NULL;
     int storedist = 0; /* 0 for STORE, 1 for STOREDIST. */
 
     /* Look up the requested zset */
     robj *zobj = NULL;
-    if ((zobj = lookupKeyReadOrReply(c, key, shared.emptyarray)) == NULL ||
+    if ((zobj = lookupKeyReadOrReply(c, c->argv[srcKeyIndex], shared.emptyarray)) == NULL ||
         checkType(c, zobj, OBJ_ZSET)) {
         return;
     }
 
-    /* Find long/lat to use for radius search based on inquiry type */
+    /* Find long/lat to use for radius or box search based on inquiry type */
     int base_args;
-    double xy[2] = { 0 };
+    GeoShape shape = {0};
     if (flags & RADIUS_COORDS) {
         base_args = 6;
-        if (extractLongLatOrReply(c, c->argv + 2, xy) == C_ERR)
-            return;
+        shape.type = CIRCULAR_TYPE;
+        if (extractLongLatOrReply(c, c->argv + 2, shape.xy) == C_ERR) return;
+        if (extractDistanceOrReply(c, c->argv+base_args-2, &shape.conversion, &shape.t.radius) != C_OK) return;
     } else if (flags & RADIUS_MEMBER) {
         base_args = 5;
+        shape.type = CIRCULAR_TYPE;
         robj *member = c->argv[2];
-        if (longLatFromMember(zobj, member, xy) == C_ERR) {
+        if (longLatFromMember(zobj, member, shape.xy) == C_ERR) {
             addReplyError(c, "could not decode requested zset member");
             return;
         }
+        if (extractDistanceOrReply(c, c->argv+base_args-2, &shape.conversion, &shape.t.radius) != C_OK) return;
+    } else if (flags & GEOSEARCH) {
+        base_args = 2;
+        if (flags & GEOSEARCHSTORE) {
+            base_args = 3;
+            storekey = c->argv[1];
+        }
     } else {
         addReplyError(c, "Unknown georadius search type");
         return;
     }
 
-    /* Extract radius and units from arguments */
-    double radius_meters = 0, conversion = 1;
-    if ((radius_meters = extractDistanceOrReply(c, c->argv + base_args - 2,
-                                                &conversion)) < 0) {
-        return;
-    }
-
     /* Discover and populate all optional parameters. */
     int withdist = 0, withhash = 0, withcoords = 0;
+    int frommember = 0, fromloc = 0, byradius = 0, bybox = 0;
     int sort = SORT_NONE;
     long long count = 0;
     if (c->argc > base_args) {
@@ -525,18 +561,64 @@ void georadiusGeneric(client *c, int flags) {
                 i++;
             } else if (!strcasecmp(arg, "store") &&
                        (i+1) < remaining &&
-                       !(flags & RADIUS_NOSTORE))
+                       !(flags & RADIUS_NOSTORE) &&
+                       !(flags & GEOSEARCH))
             {
                 storekey = c->argv[base_args+i+1];
                 storedist = 0;
                 i++;
             } else if (!strcasecmp(arg, "storedist") &&
                        (i+1) < remaining &&
-                       !(flags & RADIUS_NOSTORE))
+                       !(flags & RADIUS_NOSTORE) &&
+                       !(flags & GEOSEARCH))
             {
                 storekey = c->argv[base_args+i+1];
                 storedist = 1;
                 i++;
+            } else if (!strcasecmp(arg, "storedist") &&
+                       (flags & GEOSEARCH) &&
+                       (flags & GEOSEARCHSTORE))
+            {
+                storedist = 1;
+            } else if (!strcasecmp(arg, "frommember") &&
+                      (i+1) < remaining &&
+                      flags & GEOSEARCH &&
+                      !fromloc)
+            {
+                if (longLatFromMember(zobj, c->argv[base_args+i+1], shape.xy) == C_ERR) {
+                    addReplyError(c, "could not decode requested zset member");
+                    return;
+                }
+                frommember = 1;
+                i++;
+            } else if (!strcasecmp(arg, "fromloc") &&
+                       (i+2) < remaining &&
+                       flags & GEOSEARCH &&
+                       !frommember)
+            {
+                if (extractLongLatOrReply(c, c->argv+base_args+i+1, shape.xy) == C_ERR) return;
+                fromloc = 1;
+                i += 2;
+            } else if (!strcasecmp(arg, "byradius") &&
+                       (i+2) < remaining &&
+                       flags & GEOSEARCH &&
+                       !bybox)
+            {
+                if (extractDistanceOrReply(c, c->argv+base_args+i+1, &shape.conversion, &shape.t.radius) != C_OK)
+                    return;
+                shape.type = CIRCULAR_TYPE;
+                byradius = 1;
+                i += 2;
+            } else if (!strcasecmp(arg, "bybox") &&
+                       (i+3) < remaining &&
+                       flags & GEOSEARCH &&
+                       !byradius)
+            {
+                if (extractBoxOrReply(c, c->argv+base_args+i+1, &shape.conversion, &shape.t.r.height,
+                        &shape.t.r.width) != C_OK) return;
+                shape.type = RECTANGLE_TYPE;
+                bybox = 1;
+                i += 3;
             } else {
                 addReply(c, shared.syntaxerr);
                 return;
@@ -552,17 +634,26 @@ void georadiusGeneric(client *c, int flags) {
         return;
     }
 
+    if ((flags & GEOSEARCH) && !(frommember || fromloc)) {
+        addReplyErrorFormat(c, "exactly one of FROMMEMBER or FROMLOC can be specified for %s", (char *)c->argv[0]->ptr);
+        return;
+    }
+
+    if ((flags & GEOSEARCH) && !(byradius || bybox)) {
+        addReplyErrorFormat(c, "exactly one of BYRADIUS and BYBOX can be specified for %s", (char *)c->argv[0]->ptr);
+        return;
+    }
+
     /* COUNT without ordering does not make much sense, force ASC
      * ordering if COUNT was specified but no sorting was requested. */
     if (count != 0 && sort == SORT_NONE) sort = SORT_ASC;
 
     /* Get all neighbor geohash boxes for our radius search */
-    GeoHashRadius georadius =
-        geohashGetAreasByRadiusWGS84(xy[0], xy[1], radius_meters);
+    GeoHashRadius georadius = geohashCalculateAreasByShapeWGS84(&shape);
 
     /* Search the zset for all matching points */
     geoArray *ga = geoArrayCreate();
-    membersOfAllNeighbors(zobj, georadius, xy[0], xy[1], radius_meters, ga);
+    membersOfAllNeighbors(zobj, georadius, &shape, ga);
 
     /* If no matching results, the user gets an empty reply. */
     if (ga->used == 0 && storekey == NULL) {
@@ -607,7 +698,7 @@ void georadiusGeneric(client *c, int flags) {
         int i;
         for (i = 0; i < returned_items; i++) {
             geoPoint *gp = ga->array+i;
-            gp->dist /= conversion; /* Fix according to unit. */
+            gp->dist /= shape.conversion; /* Fix according to unit. */
 
             /* If we have options in option_length, return each sub-result
              * as a nested multi-bulk.  Add 1 to account for result value
@@ -645,7 +736,7 @@ void georadiusGeneric(client *c, int flags) {
         for (i = 0; i < returned_items; i++) {
             zskiplistNode *znode;
             geoPoint *gp = ga->array+i;
-            gp->dist /= conversion; /* Fix according to unit. */
+            gp->dist /= shape.conversion; /* Fix according to unit. */
             double score = storedist ? gp->dist : gp->score;
             size_t elelen = sdslen(gp->member);
 
@@ -659,7 +750,7 @@ void georadiusGeneric(client *c, int flags) {
             zsetConvertToZiplistIfNeeded(zobj,maxelelen);
             setKey(c,c->db,storekey,zobj);
             decrRefCount(zobj);
-            notifyKeyspaceEvent(NOTIFY_ZSET,"georadiusstore",storekey,
+            notifyKeyspaceEvent(NOTIFY_ZSET,flags & GEOSEARCH ? "geosearchstore" : "georadiusstore",storekey,
                                 c->db->id);
             server.dirty += returned_items;
         } else if (dbDelete(c->db,storekey)) {
@@ -674,22 +765,30 @@ void georadiusGeneric(client *c, int flags) {
 
 /* GEORADIUS wrapper function. */
 void georadiusCommand(client *c) {
-    georadiusGeneric(c, RADIUS_COORDS);
+    georadiusGeneric(c, 1, RADIUS_COORDS);
 }
 
 /* GEORADIUSBYMEMBER wrapper function. */
 void georadiusbymemberCommand(client *c) {
-    georadiusGeneric(c, RADIUS_MEMBER);
+    georadiusGeneric(c, 1, RADIUS_MEMBER);
 }
 
 /* GEORADIUS_RO wrapper function. */
 void georadiusroCommand(client *c) {
-    georadiusGeneric(c, RADIUS_COORDS|RADIUS_NOSTORE);
+    georadiusGeneric(c, 1, RADIUS_COORDS|RADIUS_NOSTORE);
 }
 
 /* GEORADIUSBYMEMBER_RO wrapper function. */
 void georadiusbymemberroCommand(client *c) {
-    georadiusGeneric(c, RADIUS_MEMBER|RADIUS_NOSTORE);
+    georadiusGeneric(c, 1, RADIUS_MEMBER|RADIUS_NOSTORE);
+}
+
+void geosearchCommand(client *c) {
+    georadiusGeneric(c, 1, GEOSEARCH);
+}
+
+void geosearchstoreCommand(client *c) {
+    georadiusGeneric(c, 2, GEOSEARCH|GEOSEARCHSTORE);
 }
 
 /* GEOHASH key ele1 ele2 ... eleN
diff --git a/src/geohash.h b/src/geohash.h
index ed2ef9336..ebe52ff34 100644
--- a/src/geohash.h
+++ b/src/geohash.h
@@ -90,6 +90,25 @@ typedef struct {
     GeoHashBits south_west;
 } GeoHashNeighbors;
 
+#define CIRCULAR_TYPE 1
+#define RECTANGLE_TYPE 2
+typedef struct {
+    int type; /* search type */
+    double xy[2]; /* search center point, xy[0]: lon, xy[1]: lat */
+    double conversion; /* km: 1000 */
+    double bounds[4]; /* bounds[0]: min_lon, bounds[1]: min_lat
+                       * bounds[2]: max_lon, bounds[3]: max_lat */
+    union {
+        /* CIRCULAR_TYPE */
+        double radius;
+        /* RECTANGLE_TYPE */
+        struct {
+            double height;
+            double width;
+        } r;
+    } t;
+} GeoShape;
+
 /*
  * 0:success
  * -1:failed
diff --git a/src/geohash_helper.c b/src/geohash_helper.c
index 01fb2cb88..a1ddf5d31 100644
--- a/src/geohash_helper.c
+++ b/src/geohash_helper.c
@@ -82,10 +82,9 @@ uint8_t geohashEstimateStepsByRadius(double range_meters, double lat) {
     return step;
 }
 
-/* Return the bounding box of the search area centered at latitude,longitude
- * having a radius of radius_meter. bounds[0] - bounds[2] is the minimum
- * and maximum longitude, while bounds[1] - bounds[3] is the minimum and
- * maximum latitude.
+/* Return the bounding box of the search area by shape (see geohash.h GeoShape)
+ * bounds[0] - bounds[2] is the minimum and maximum longitude
+ * while bounds[1] - bounds[3] is the minimum and maximum latitude.
  *
  * This function does not behave correctly with very large radius values, for
  * instance for the coordinates 81.634948934258375 30.561509253718668 and a
@@ -100,34 +99,51 @@ uint8_t geohashEstimateStepsByRadius(double range_meters, double lat) {
  * Since this function is currently only used as an optimization, the
  * optimization is not used for very big radiuses, however the function
  * should be fixed. */
-int geohashBoundingBox(double longitude, double latitude, double radius_meters,
-                       double *bounds) {
+int geohashBoundingBox(GeoShape *shape, double *bounds) {
     if (!bounds) return 0;
+    double longitude = shape->xy[0];
+    double latitude = shape->xy[1];
+    double height = shape->conversion * (shape->type == CIRCULAR_TYPE ? shape->t.radius : shape->t.r.height/2);
+    double width = shape->conversion * (shape->type == CIRCULAR_TYPE ? shape->t.radius : shape->t.r.width/2);
 
-    bounds[0] = longitude - rad_deg(radius_meters/EARTH_RADIUS_IN_METERS/cos(deg_rad(latitude)));
-    bounds[2] = longitude + rad_deg(radius_meters/EARTH_RADIUS_IN_METERS/cos(deg_rad(latitude)));
-    bounds[1] = latitude - rad_deg(radius_meters/EARTH_RADIUS_IN_METERS);
-    bounds[3] = latitude + rad_deg(radius_meters/EARTH_RADIUS_IN_METERS);
+    const double long_delta = rad_deg(width/EARTH_RADIUS_IN_METERS/cos(deg_rad(latitude)));
+    const double lat_delta = rad_deg(height/EARTH_RADIUS_IN_METERS);
+    bounds[0] = longitude - long_delta;
+    bounds[2] = longitude + long_delta;
+    bounds[1] = latitude - lat_delta;
+    bounds[3] = latitude + lat_delta;
     return 1;
 }
 
-/* Return a set of areas (center + 8) that are able to cover a range query
- * for the specified position and radius. */
-GeoHashRadius geohashGetAreasByRadius(double longitude, double latitude, double radius_meters) {
+/* Calculate a set of areas (center + 8) that are able to cover a range query
+ * for the specified position and shape (see geohash.h GeoShape).
+ * the bounding box saved in shaple.bounds */
+GeoHashRadius geohashCalculateAreasByShapeWGS84(GeoShape *shape) {
     GeoHashRange long_range, lat_range;
     GeoHashRadius radius;
     GeoHashBits hash;
     GeoHashNeighbors neighbors;
     GeoHashArea area;
     double min_lon, max_lon, min_lat, max_lat;
-    double bounds[4];
     int steps;
 
-    geohashBoundingBox(longitude, latitude, radius_meters, bounds);
-    min_lon = bounds[0];
-    min_lat = bounds[1];
-    max_lon = bounds[2];
-    max_lat = bounds[3];
+    geohashBoundingBox(shape, shape->bounds);
+    min_lon = shape->bounds[0];
+    min_lat = shape->bounds[1];
+    max_lon = shape->bounds[2];
+    max_lat = shape->bounds[3];
+
+    double longitude = shape->xy[0];
+    double latitude = shape->xy[1];
+    /* radius_meters is calculated differently in different search types:
+     * 1) CIRCULAR_TYPE, just use radius.
+     * 2) RECTANGLE_TYPE, in order to calculate accurately, we should use
+     * sqrt((width/2)^2 + (height/2)^2), so that the box is bound by a circle,
+     * But the current code a simpler approach resulting in a smaller circle,
+     * which is safe because we search the 8 nearby boxes anyway. */
+    double radius_meters = shape->type == CIRCULAR_TYPE ? shape->t.radius :
+                            shape->t.r.width > shape->t.r.height ? shape->t.r.width/2 : shape->t.r.height/2;
+    radius_meters *= shape->conversion;
 
     steps = geohashEstimateStepsByRadius(radius_meters,latitude);
 
@@ -196,11 +212,6 @@ GeoHashRadius geohashGetAreasByRadius(double longitude, double latitude, double
     return radius;
 }
 
-GeoHashRadius geohashGetAreasByRadiusWGS84(double longitude, double latitude,
-                                           double radius_meters) {
-    return geohashGetAreasByRadius(longitude, latitude, radius_meters);
-}
-
 GeoHashFix52Bits geohashAlign52Bits(const GeoHashBits hash) {
     uint64_t bits = hash.bits;
     bits <<= (52 - hash.step * 2);
@@ -233,3 +244,15 @@ int geohashGetDistanceIfInRadiusWGS84(double x1, double y1, double x2,
                                       double *distance) {
     return geohashGetDistanceIfInRadius(x1, y1, x2, y2, radius, distance);
 }
+
+/* Judge whether a point is in the axis-aligned rectangle.
+ * bounds : see geohash.h GeoShape::bounds
+ * x1, y1 : the center of the box
+ * x2, y2 : the point to be searched
+ */
+int geohashGetDistanceIfInRectangle(double *bounds, double x1, double y1,
+                                    double x2, double y2, double *distance) {
+    if (x2 < bounds[0] || x2 > bounds[2] || y2 < bounds[1] || y2 > bounds[3]) return 0;
+    *distance = geohashGetDistance(x1, y1, x2, y2);
+    return 1;
+}
diff --git a/src/geohash_helper.h b/src/geohash_helper.h
index eb0dda38a..f28896a8d 100644
--- a/src/geohash_helper.h
+++ b/src/geohash_helper.h
@@ -49,14 +49,8 @@ typedef struct {
 
 int GeoHashBitsComparator(const GeoHashBits *a, const GeoHashBits *b);
 uint8_t geohashEstimateStepsByRadius(double range_meters, double lat);
-int geohashBoundingBox(double longitude, double latitude, double radius_meters,
-                        double *bounds);
-GeoHashRadius geohashGetAreasByRadius(double longitude,
-                                      double latitude, double radius_meters);
-GeoHashRadius geohashGetAreasByRadiusWGS84(double longitude, double latitude,
-                                           double radius_meters);
-GeoHashRadius geohashGetAreasByRadiusMercator(double longitude, double latitude,
-                                              double radius_meters);
+int geohashBoundingBox(GeoShape *shape, double *bounds);
+GeoHashRadius geohashCalculateAreasByShapeWGS84(GeoShape *shape);
 GeoHashFix52Bits geohashAlign52Bits(const GeoHashBits hash);
 double geohashGetDistance(double lon1d, double lat1d,
                           double lon2d, double lat2d);
@@ -66,5 +60,7 @@ int geohashGetDistanceIfInRadius(double x1, double y1,
 int geohashGetDistanceIfInRadiusWGS84(double x1, double y1, double x2,
                                       double y2, double radius,
                                       double *distance);
+int geohashGetDistanceIfInRectangle(double *bounds, double x1, double y1,
+                                    double x2, double y2, double *distance);
 
 #endif /* GEOHASH_HELPER_HPP_ */
diff --git a/src/server.c b/src/server.c
index 0c3bbff03..e4d236874 100644
--- a/src/server.c
+++ b/src/server.c
@@ -919,7 +919,7 @@ struct redisCommand redisCommandTable[] = {
 
     /* GEORADIUS has store options that may write. */
     {"georadius",georadiusCommand,-6,
-     "write @geo",
+     "write use-memory @geo",
      0,georadiusGetKeys,1,1,1,0,0,0},
 
     {"georadius_ro",georadiusroCommand,-6,
@@ -927,7 +927,7 @@ struct redisCommand redisCommandTable[] = {
      0,georadiusGetKeys,1,1,1,0,0,0},
 
     {"georadiusbymember",georadiusbymemberCommand,-5,
-     "write @geo",
+     "write use-memory @geo",
      0,georadiusGetKeys,1,1,1,0,0,0},
 
     {"georadiusbymember_ro",georadiusbymemberroCommand,-5,
@@ -946,6 +946,14 @@ struct redisCommand redisCommandTable[] = {
      "read-only @geo",
      0,NULL,1,1,1,0,0,0},
 
+    {"geosearch",geosearchCommand,-7,
+     "read-only @geo",
+      0,NULL,1,1,1,0,0,0},
+
+    {"geosearchstore",geosearchstoreCommand,-8,
+     "write use-memory @geo",
+      0,NULL,1,2,1,0,0,0},
+
     {"pfselftest",pfselftestCommand,1,
      "admin @hyperloglog",
       0,NULL,0,0,0,0,0,0},
diff --git a/src/server.h b/src/server.h
index 35b695a42..4cb5ab855 100644
--- a/src/server.h
+++ b/src/server.h
@@ -2546,6 +2546,8 @@ void geoaddCommand(client *c);
 void geohashCommand(client *c);
 void geoposCommand(client *c);
 void geodistCommand(client *c);
+void geosearchCommand(client *c);
+void geosearchstoreCommand(client *c);
 void pfselftestCommand(client *c);
 void pfaddCommand(client *c);
 void pfcountCommand(client *c);
diff --git a/tests/unit/geo.tcl b/tests/unit/geo.tcl
index 49e421ee9..b3e59d9b2 100644
--- a/tests/unit/geo.tcl
+++ b/tests/unit/geo.tcl
@@ -94,10 +94,43 @@ start_server {tags {"geo"}} {
         r georadius nyc -73.9798091 40.7598464 3 km asc
     } {{central park n/q/r} 4545 {union square}}
 
+    test {GEOSEARCH simple (sorted)} {
+        r geosearch nyc fromloc -73.9798091 40.7598464 bybox 6 6 km asc
+    } {{central park n/q/r} 4545 {union square} {lic market}}
+
+    test {GEOSEARCH FROMLOC and FROMMEMBER cannot exist at the same time} {
+        catch {r geosearch nyc fromloc -73.9798091 40.7598464 frommember xxx bybox 6 6 km asc} e
+        set e
+    } {ERR*syntax*}
+
+    test {GEOSEARCH FROMLOC and FROMMEMBER one must exist} {
+        catch {r geosearch nyc bybox 3 3 km asc desc withhash withdist withcoord} e
+        set e
+    } {ERR*exactly one of FROMMEMBER or FROMLOC*}
+
+    test {GEOSEARCH BYRADIUS and BYBOX cannot exist at the same time} {
+        catch {r geosearch nyc fromloc -73.9798091 40.7598464 byradius 3 km bybox 3 3 km asc} e
+        set e
+    } {ERR*syntax*}
+
+    test {GEOSEARCH BYRADIUS and BYBOX one must exist} {
+        catch {r geosearch nyc fromloc -73.9798091 40.7598464 asc desc withhash withdist withcoord} e
+        set e
+    } {ERR*exactly one of BYRADIUS and BYBOX*}
+
+    test {GEOSEARCH with STOREDIST option} {
+        catch {r geosearch nyc fromloc -73.9798091 40.7598464 bybox 6 6 km asc storedist} e
+        set e
+    } {ERR*syntax*}
+
     test {GEORADIUS withdist (sorted)} {
         r georadius nyc -73.9798091 40.7598464 3 km withdist asc
     } {{{central park n/q/r} 0.7750} {4545 2.3651} {{union square} 2.7697}}
 
+    test {GEOSEARCH withdist (sorted)} {
+        r geosearch nyc fromloc -73.9798091 40.7598464 bybox 6 6 km withdist asc
+    } {{{central park n/q/r} 0.7750} {4545 2.3651} {{union square} 2.7697} {{lic market} 3.1991}}
+
     test {GEORADIUS with COUNT} {
         r georadius nyc -73.9798091 40.7598464 10 km COUNT 3
     } {{central park n/q/r} 4545 {union square}}
@@ -120,6 +153,35 @@ start_server {tags {"geo"}} {
         r georadiusbymember nyc "wtc one" 7 km
     } {{wtc one} {union square} {central park n/q/r} 4545 {lic market}}
 
+    test {GEOSEARCH FROMMEMBER simple (sorted)} {
+        r geosearch nyc frommember "wtc one" bybox 14 14 km
+    } {{wtc one} {union square} {central park n/q/r} 4545 {lic market} q4}
+
+    test {GEOSEARCH vs GEORADIUS} {
+        r del Sicily
+        r geoadd Sicily 13.361389 38.115556 "Palermo" 15.087269 37.502669 "Catania"
+        r geoadd Sicily 12.758489 38.788135 "edge1"   17.241510 38.788135 "eage2"
+        set ret1 [r georadius Sicily 15 37 200 km asc]
+        assert_equal $ret1 {Catania Palermo}
+        set ret2 [r geosearch Sicily fromloc 15 37 bybox 400 400 km asc]
+        assert_equal $ret2 {Catania Palermo eage2 edge1}
+    }
+
+    test {GEOSEARCH non square, long and narrow} {
+        r del Sicily
+        r geoadd Sicily 12.75 37.00 "test1"
+        r geoadd Sicily 12.75 36.50 "test2"
+        r geoadd Sicily 13.00 36.50 "test3"
+        # box height=2km width=400km
+        set ret1 [r geosearch Sicily fromloc 15 37 bybox 2 400 km]
+        assert_equal $ret1 {test1}
+
+        # Add a western Hemisphere point
+        r geoadd Sicily -1 37.00 "test3"
+        set ret2 [r geosearch Sicily fromloc 15 37 bybox 2 3000 km asc]
+        assert_equal $ret2 {test1 test3}
+    }
+
     test {GEORADIUSBYMEMBER withdist (sorted)} {
         r georadiusbymember nyc "wtc one" 7 km withdist
     } {{{wtc one} 0.0000} {{union square} 3.2544} {{central park n/q/r} 6.7000} {4545 6.1975} {{lic market} 6.8969}}
@@ -178,6 +240,11 @@ start_server {tags {"geo"}} {
         set e
     } {*ERR*syntax*}
 
+    test {GEOSEARCHSTORE STORE option: syntax error} {
+        catch {r geosearchstore abc points fromloc 13.361389 38.115556 byradius 50 km store abc} e
+        set e
+    } {*ERR*syntax*}
+
     test {GEORANGE STORE option: incompatible options} {
         r del points
         r geoadd points 13.361389 38.115556 "Palermo" \
@@ -198,6 +265,11 @@ start_server {tags {"geo"}} {
         assert_equal [r zrange points 0 -1] [r zrange points2 0 -1]
     }
 
+    test {GEOSEARCHSTORE STORE option: plain usage} {
+        r geosearchstore points2 points fromloc 13.361389 38.115556 byradius 500 km
+        assert_equal [r zrange points 0 -1] [r zrange points2 0 -1]
+    }
+
     test {GEORANGE STOREDIST option: plain usage} {
         r del points
         r geoadd points 13.361389 38.115556 "Palermo" \
@@ -209,6 +281,14 @@ start_server {tags {"geo"}} {
         assert {[lindex $res 3] < 167}
     }
 
+    test {GEOSEARCHSTORE STOREDIST option: plain usage} {
+        r geosearchstore points2 points fromloc 13.361389 38.115556 byradius 500 km storedist
+        set res [r zrange points2 0 -1 withscores]
+        assert {[lindex $res 1] < 1}
+        assert {[lindex $res 3] > 166}
+        assert {[lindex $res 3] < 167}
+    }
+
     test {GEORANGE STOREDIST option: COUNT ASC and DESC} {
         r del points
         r geoadd points 13.361389 38.115556 "Palermo" \