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 torwards 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 centered at latitude,longitude
86 * having a radius of radius_meter. bounds[0] - bounds[2] is the minimum
87 * and maxium longitude, while bounds[1] - bounds[3] is the minimum and
88 * maximum latitude.
89 *
90 * This function does not behave correctly with very large radius values, for
91 * instance for the coordinates 81.634948934258375 30.561509253718668 and a
92 * radius of 7083 kilometers, it reports as bounding boxes:
93 *
94 * min_lon 7.680495, min_lat -33.119473, max_lon 155.589402, max_lat 94.242491
95 *
96 * However, for instance, a min_lon of 7.680495 is not correct, because the
97 * point -1.27579540014266968 61.33421815228281559 is at less than 7000
98 * kilometers away.
99 *
100 * Since this function is currently only used as an optimization, the
101 * optimization is not used for very big radiuses, however the function
102 * should be fixed. */
geohashBoundingBox(double longitude,double latitude,double radius_meters,double * bounds)103 int geohashBoundingBox(double longitude, double latitude, double radius_meters,
104 double *bounds) {
105 if (!bounds) return 0;
106
107 bounds[0] = longitude - rad_deg(radius_meters/EARTH_RADIUS_IN_METERS/cos(deg_rad(latitude)));
108 bounds[2] = longitude + rad_deg(radius_meters/EARTH_RADIUS_IN_METERS/cos(deg_rad(latitude)));
109 bounds[1] = latitude - rad_deg(radius_meters/EARTH_RADIUS_IN_METERS);
110 bounds[3] = latitude + rad_deg(radius_meters/EARTH_RADIUS_IN_METERS);
111 return 1;
112 }
113
114 /* Return a set of areas (center + 8) that are able to cover a range query
115 * for the specified position and radius. */
geohashGetAreasByRadius(double longitude,double latitude,double radius_meters)116 GeoHashRadius geohashGetAreasByRadius(double longitude, double latitude, double radius_meters) {
117 GeoHashRange long_range, lat_range;
118 GeoHashRadius radius;
119 GeoHashBits hash;
120 GeoHashNeighbors neighbors;
121 GeoHashArea area;
122 double min_lon, max_lon, min_lat, max_lat;
123 double bounds[4];
124 int steps;
125
126 geohashBoundingBox(longitude, latitude, radius_meters, bounds);
127 min_lon = bounds[0];
128 min_lat = bounds[1];
129 max_lon = bounds[2];
130 max_lat = bounds[3];
131
132 steps = geohashEstimateStepsByRadius(radius_meters,latitude);
133
134 geohashGetCoordRange(&long_range,&lat_range);
135 geohashEncode(&long_range,&lat_range,longitude,latitude,steps,&hash);
136 geohashNeighbors(&hash,&neighbors);
137 geohashDecode(long_range,lat_range,hash,&area);
138
139 /* Check if the step is enough at the limits of the covered area.
140 * Sometimes when the search area is near an edge of the
141 * area, the estimated step is not small enough, since one of the
142 * north / south / west / east square is too near to the search area
143 * to cover everything. */
144 int decrease_step = 0;
145 {
146 GeoHashArea north, south, east, west;
147
148 geohashDecode(long_range, lat_range, neighbors.north, &north);
149 geohashDecode(long_range, lat_range, neighbors.south, &south);
150 geohashDecode(long_range, lat_range, neighbors.east, &east);
151 geohashDecode(long_range, lat_range, neighbors.west, &west);
152
153 if (geohashGetDistance(longitude,latitude,longitude,north.latitude.max)
154 < radius_meters) decrease_step = 1;
155 if (geohashGetDistance(longitude,latitude,longitude,south.latitude.min)
156 < radius_meters) decrease_step = 1;
157 if (geohashGetDistance(longitude,latitude,east.longitude.max,latitude)
158 < radius_meters) decrease_step = 1;
159 if (geohashGetDistance(longitude,latitude,west.longitude.min,latitude)
160 < radius_meters) decrease_step = 1;
161 }
162
163 if (steps > 1 && decrease_step) {
164 steps--;
165 geohashEncode(&long_range,&lat_range,longitude,latitude,steps,&hash);
166 geohashNeighbors(&hash,&neighbors);
167 geohashDecode(long_range,lat_range,hash,&area);
168 }
169
170 /* Exclude the search areas that are useless. */
171 if (steps >= 2) {
172 if (area.latitude.min < min_lat) {
173 GZERO(neighbors.south);
174 GZERO(neighbors.south_west);
175 GZERO(neighbors.south_east);
176 }
177 if (area.latitude.max > max_lat) {
178 GZERO(neighbors.north);
179 GZERO(neighbors.north_east);
180 GZERO(neighbors.north_west);
181 }
182 if (area.longitude.min < min_lon) {
183 GZERO(neighbors.west);
184 GZERO(neighbors.south_west);
185 GZERO(neighbors.north_west);
186 }
187 if (area.longitude.max > max_lon) {
188 GZERO(neighbors.east);
189 GZERO(neighbors.south_east);
190 GZERO(neighbors.north_east);
191 }
192 }
193 radius.hash = hash;
194 radius.neighbors = neighbors;
195 radius.area = area;
196 return radius;
197 }
198
geohashGetAreasByRadiusWGS84(double longitude,double latitude,double radius_meters)199 GeoHashRadius geohashGetAreasByRadiusWGS84(double longitude, double latitude,
200 double radius_meters) {
201 return geohashGetAreasByRadius(longitude, latitude, radius_meters);
202 }
203
geohashAlign52Bits(const GeoHashBits hash)204 GeoHashFix52Bits geohashAlign52Bits(const GeoHashBits hash) {
205 uint64_t bits = hash.bits;
206 bits <<= (52 - hash.step * 2);
207 return bits;
208 }
209
210 /* Calculate distance using haversin great circle distance formula. */
geohashGetDistance(double lon1d,double lat1d,double lon2d,double lat2d)211 double geohashGetDistance(double lon1d, double lat1d, double lon2d, double lat2d) {
212 double lat1r, lon1r, lat2r, lon2r, u, v;
213 lat1r = deg_rad(lat1d);
214 lon1r = deg_rad(lon1d);
215 lat2r = deg_rad(lat2d);
216 lon2r = deg_rad(lon2d);
217 u = sin((lat2r - lat1r) / 2);
218 v = sin((lon2r - lon1r) / 2);
219 return 2.0 * EARTH_RADIUS_IN_METERS *
220 asin(sqrt(u * u + cos(lat1r) * cos(lat2r) * v * v));
221 }
222
geohashGetDistanceIfInRadius(double x1,double y1,double x2,double y2,double radius,double * distance)223 int geohashGetDistanceIfInRadius(double x1, double y1,
224 double x2, double y2, double radius,
225 double *distance) {
226 *distance = geohashGetDistance(x1, y1, x2, y2);
227 if (*distance > radius) return 0;
228 return 1;
229 }
230
geohashGetDistanceIfInRadiusWGS84(double x1,double y1,double x2,double y2,double radius,double * distance)231 int geohashGetDistanceIfInRadiusWGS84(double x1, double y1, double x2,
232 double y2, double radius,
233 double *distance) {
234 return geohashGetDistanceIfInRadius(x1, y1, x2, y2, radius, distance);
235 }
236