1 /*
2  * Copyright (c) 2013-2014, yinqiwen <yinqiwen@gmail.com>
3  * Copyright (c) 2014, Matt Stancliff <matt@genges.com>.
4  * Copyright (c) 2015-2016, Salvatore Sanfilippo <antirez@gmail.com>.
5  * All rights reserved.
6  *
7  * Redistribution and use in source and binary forms, with or without
8  * modification, are permitted provided that the following conditions are met:
9  *
10  *  * Redistributions of source code must retain the above copyright notice,
11  *    this list of conditions and the following disclaimer.
12  *  * Redistributions in binary form must reproduce the above copyright
13  *    notice, this list of conditions and the following disclaimer in the
14  *    documentation and/or other materials provided with the distribution.
15  *  * Neither the name of Redis nor the names of its contributors may be used
16  *    to endorse or promote products derived from this software without
17  *    specific prior written permission.
18  *
19  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS
23  * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
29  * THE POSSIBILITY OF SUCH DAMAGE.
30  */
31 
32 /* This is a C++ to C conversion from the ardb project.
33  * This file started out as:
34  * https://github.com/yinqiwen/ardb/blob/d42503/src/geo/geohash_helper.cpp
35  */
36 
37 #include "fmacros.h"
38 #include "geohash_helper.h"
39 #include "debugmacro.h"
40 #include <math.h>
41 
42 #define D_R (M_PI / 180.0)
43 #define R_MAJOR 6378137.0
44 #define R_MINOR 6356752.3142
45 #define RATIO (R_MINOR / R_MAJOR)
46 #define ECCENT (sqrt(1.0 - (RATIO *RATIO)))
47 #define COM (0.5 * ECCENT)
48 
49 /// @brief The usual PI/180 constant
50 const double DEG_TO_RAD = 0.017453292519943295769236907684886;
51 /// @brief Earth's quatratic mean radius for WGS-84
52 const double EARTH_RADIUS_IN_METERS = 6372797.560856;
53 
54 const double MERCATOR_MAX = 20037726.37;
55 const double MERCATOR_MIN = -20037726.37;
56 
deg_rad(double ang)57 static inline double deg_rad(double ang) { return ang * D_R; }
rad_deg(double ang)58 static inline double rad_deg(double ang) { return ang / D_R; }
59 
60 /* This function is used in order to estimate the step (bits precision)
61  * of the 9 search area boxes during radius queries. */
geohashEstimateStepsByRadius(double range_meters,double lat)62 uint8_t geohashEstimateStepsByRadius(double range_meters, double lat) {
63     if (range_meters == 0) return 26;
64     int step = 1;
65     while (range_meters < MERCATOR_MAX) {
66         range_meters *= 2;
67         step++;
68     }
69     step -= 2; /* Make sure range is included in most of the base cases. */
70 
71     /* Wider range towards the poles... Note: it is possible to do better
72      * than this approximation by computing the distance between meridians
73      * at this latitude, but this does the trick for now. */
74     if (lat > 66 || lat < -66) {
75         step--;
76         if (lat > 80 || lat < -80) step--;
77     }
78 
79     /* Frame to valid range. */
80     if (step < 1) step = 1;
81     if (step > 26) step = 26;
82     return step;
83 }
84 
85 /* Return the bounding box of the search area by shape (see geohash.h GeoShape)
86  * bounds[0] - bounds[2] is the minimum and maximum longitude
87  * while bounds[1] - bounds[3] is the minimum and maximum latitude.
88  * since the higher the latitude, the shorter the arc length, the box shape is as follows
89  * (left and right edges are actually bent), as shown in the following diagram:
90  *
91  *    \-----------------/          --------               \-----------------/
92  *     \               /         /          \              \               /
93  *      \  (long,lat) /         / (long,lat) \              \  (long,lat) /
94  *       \           /         /              \              /            \
95  *         ---------          /----------------\            /--------------\
96  *  Northern Hemisphere       Southern Hemisphere         Around the equator
97  */
geohashBoundingBox(GeoShape * shape,double * bounds)98 int geohashBoundingBox(GeoShape *shape, double *bounds) {
99     if (!bounds) return 0;
100     double longitude = shape->xy[0];
101     double latitude = shape->xy[1];
102     double height = shape->conversion * (shape->type == CIRCULAR_TYPE ? shape->t.radius : shape->t.r.height/2);
103     double width = shape->conversion * (shape->type == CIRCULAR_TYPE ? shape->t.radius : shape->t.r.width/2);
104 
105     const double lat_delta = rad_deg(height/EARTH_RADIUS_IN_METERS);
106     const double long_delta_top = rad_deg(width/EARTH_RADIUS_IN_METERS/cos(deg_rad(latitude+lat_delta)));
107     const double long_delta_bottom = rad_deg(width/EARTH_RADIUS_IN_METERS/cos(deg_rad(latitude-lat_delta)));
108     /* The directions of the northern and southern hemispheres
109      * are opposite, so we choice different points as min/max long/lat */
110     int southern_hemisphere = latitude < 0 ? 1 : 0;
111     bounds[0] = southern_hemisphere ? longitude-long_delta_bottom : longitude-long_delta_top;
112     bounds[2] = southern_hemisphere ? longitude+long_delta_bottom : longitude+long_delta_top;
113     bounds[1] = latitude - lat_delta;
114     bounds[3] = latitude + lat_delta;
115     return 1;
116 }
117 
118 /* Calculate a set of areas (center + 8) that are able to cover a range query
119  * for the specified position and shape (see geohash.h GeoShape).
120  * the bounding box saved in shaple.bounds */
geohashCalculateAreasByShapeWGS84(GeoShape * shape)121 GeoHashRadius geohashCalculateAreasByShapeWGS84(GeoShape *shape) {
122     GeoHashRange long_range, lat_range;
123     GeoHashRadius radius;
124     GeoHashBits hash;
125     GeoHashNeighbors neighbors;
126     GeoHashArea area;
127     double min_lon, max_lon, min_lat, max_lat;
128     int steps;
129 
130     geohashBoundingBox(shape, shape->bounds);
131     min_lon = shape->bounds[0];
132     min_lat = shape->bounds[1];
133     max_lon = shape->bounds[2];
134     max_lat = shape->bounds[3];
135 
136     double longitude = shape->xy[0];
137     double latitude = shape->xy[1];
138     /* radius_meters is calculated differently in different search types:
139      * 1) CIRCULAR_TYPE, just use radius.
140      * 2) RECTANGLE_TYPE, we use sqrt((width/2)^2 + (height/2)^2) to
141      * calculate the distance from the center point to the corner */
142     double radius_meters = shape->type == CIRCULAR_TYPE ? shape->t.radius :
143             sqrt((shape->t.r.width/2)*(shape->t.r.width/2) + (shape->t.r.height/2)*(shape->t.r.height/2));
144     radius_meters *= shape->conversion;
145 
146     steps = geohashEstimateStepsByRadius(radius_meters,latitude);
147 
148     geohashGetCoordRange(&long_range,&lat_range);
149     geohashEncode(&long_range,&lat_range,longitude,latitude,steps,&hash);
150     geohashNeighbors(&hash,&neighbors);
151     geohashDecode(long_range,lat_range,hash,&area);
152 
153     /* Check if the step is enough at the limits of the covered area.
154      * Sometimes when the search area is near an edge of the
155      * area, the estimated step is not small enough, since one of the
156      * north / south / west / east square is too near to the search area
157      * to cover everything. */
158     int decrease_step = 0;
159     {
160         GeoHashArea north, south, east, west;
161 
162         geohashDecode(long_range, lat_range, neighbors.north, &north);
163         geohashDecode(long_range, lat_range, neighbors.south, &south);
164         geohashDecode(long_range, lat_range, neighbors.east, &east);
165         geohashDecode(long_range, lat_range, neighbors.west, &west);
166 
167         if (geohashGetDistance(longitude,latitude,longitude,north.latitude.max)
168             < radius_meters) decrease_step = 1;
169         if (geohashGetDistance(longitude,latitude,longitude,south.latitude.min)
170             < radius_meters) decrease_step = 1;
171         if (geohashGetDistance(longitude,latitude,east.longitude.max,latitude)
172             < radius_meters) decrease_step = 1;
173         if (geohashGetDistance(longitude,latitude,west.longitude.min,latitude)
174             < radius_meters) decrease_step = 1;
175     }
176 
177     if (steps > 1 && decrease_step) {
178         steps--;
179         geohashEncode(&long_range,&lat_range,longitude,latitude,steps,&hash);
180         geohashNeighbors(&hash,&neighbors);
181         geohashDecode(long_range,lat_range,hash,&area);
182     }
183 
184     /* Exclude the search areas that are useless. */
185     if (steps >= 2) {
186         if (area.latitude.min < min_lat) {
187             GZERO(neighbors.south);
188             GZERO(neighbors.south_west);
189             GZERO(neighbors.south_east);
190         }
191         if (area.latitude.max > max_lat) {
192             GZERO(neighbors.north);
193             GZERO(neighbors.north_east);
194             GZERO(neighbors.north_west);
195         }
196         if (area.longitude.min < min_lon) {
197             GZERO(neighbors.west);
198             GZERO(neighbors.south_west);
199             GZERO(neighbors.north_west);
200         }
201         if (area.longitude.max > max_lon) {
202             GZERO(neighbors.east);
203             GZERO(neighbors.south_east);
204             GZERO(neighbors.north_east);
205         }
206     }
207     radius.hash = hash;
208     radius.neighbors = neighbors;
209     radius.area = area;
210     return radius;
211 }
212 
geohashAlign52Bits(const GeoHashBits hash)213 GeoHashFix52Bits geohashAlign52Bits(const GeoHashBits hash) {
214     uint64_t bits = hash.bits;
215     bits <<= (52 - hash.step * 2);
216     return bits;
217 }
218 
219 /* Calculate distance using haversine great circle distance formula. */
geohashGetDistance(double lon1d,double lat1d,double lon2d,double lat2d)220 double geohashGetDistance(double lon1d, double lat1d, double lon2d, double lat2d) {
221     double lat1r, lon1r, lat2r, lon2r, u, v;
222     lat1r = deg_rad(lat1d);
223     lon1r = deg_rad(lon1d);
224     lat2r = deg_rad(lat2d);
225     lon2r = deg_rad(lon2d);
226     u = sin((lat2r - lat1r) / 2);
227     v = sin((lon2r - lon1r) / 2);
228     return 2.0 * EARTH_RADIUS_IN_METERS *
229            asin(sqrt(u * u + cos(lat1r) * cos(lat2r) * v * v));
230 }
231 
geohashGetDistanceIfInRadius(double x1,double y1,double x2,double y2,double radius,double * distance)232 int geohashGetDistanceIfInRadius(double x1, double y1,
233                                  double x2, double y2, double radius,
234                                  double *distance) {
235     *distance = geohashGetDistance(x1, y1, x2, y2);
236     if (*distance > radius) return 0;
237     return 1;
238 }
239 
geohashGetDistanceIfInRadiusWGS84(double x1,double y1,double x2,double y2,double radius,double * distance)240 int geohashGetDistanceIfInRadiusWGS84(double x1, double y1, double x2,
241                                       double y2, double radius,
242                                       double *distance) {
243     return geohashGetDistanceIfInRadius(x1, y1, x2, y2, radius, distance);
244 }
245 
246 /* Judge whether a point is in the axis-aligned rectangle, when the distance
247  * between a searched point and the center point is less than or equal to
248  * height/2 or width/2 in height and width, the point is in the rectangle.
249  *
250  * width_m, height_m: the rectangle
251  * x1, y1 : the center of the box
252  * x2, y2 : the point to be searched
253  */
geohashGetDistanceIfInRectangle(double width_m,double height_m,double x1,double y1,double x2,double y2,double * distance)254 int geohashGetDistanceIfInRectangle(double width_m, double height_m, double x1, double y1,
255                                     double x2, double y2, double *distance) {
256     double lon_distance = geohashGetDistance(x2, y2, x1, y2);
257     double lat_distance = geohashGetDistance(x2, y2, x2, y1);
258     if (lon_distance > width_m/2 || lat_distance > height_m/2) {
259         return 0;
260     }
261     *distance = geohashGetDistance(x1, y1, x2, y2);
262     return 1;
263 }
264