1 /*
2 * Copyright (c) 2013-2014, yinqiwen <[email protected]>
3 * Copyright (c) 2014, Matt Stancliff <[email protected]>.
4 * Copyright (c) 2015-2016, Salvatore Sanfilippo <[email protected]>.
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 "geohash_helper.h"
38 #include <math.h>
39
40 #define D_R (M_PI / 180.0)
41 #define R_MAJOR 6378137.0
42 #define R_MINOR 6356752.3142
43 #define RATIO (R_MINOR / R_MAJOR)
44 #define ECCENT (sqrt(1.0 - (RATIO *RATIO)))
45 #define COM (0.5 * ECCENT)
46
47 /// @brief The usual PI/180 constant
48 const double DEG_TO_RAD = 0.017453292519943295769236907684886;
49 /// @brief Earth's quatratic mean radius for WGS-84
50 const double EARTH_RADIUS_IN_METERS = 6372797.560856;
51
52 const double MERCATOR_MAX = 20037726.37;
53 const double MERCATOR_MIN = -20037726.37;
54
deg_rad(double ang)55 static inline double deg_rad(double ang) { return ang * D_R; }
rad_deg(double ang)56 static inline double rad_deg(double ang) { return ang / D_R; }
57
58 /* This function is used in order to estimate the step (bits precision)
59 * of the 9 search area boxes during radius queries. */
geohashEstimateStepsByRadius(double range_meters,double lat)60 uint8_t geohashEstimateStepsByRadius(double range_meters, double lat) {
61 if (range_meters == 0) return 26;
62 int step = 1;
63 while (range_meters < MERCATOR_MAX) {
64 range_meters *= 2;
65 step++;
66 }
67 step -= 2; /* Make sure range is included in most of the base cases. */
68
69 /* Wider range torwards the poles... Note: it is possible to do better
70 * than this approximation by computing the distance between meridians
71 * at this latitude, but this does the trick for now. */
72 if (lat > 66 || lat < -66) {
73 step--;
74 if (lat > 80 || lat < -80) step--;
75 }
76
77 /* Frame to valid range. */
78 if (step < 1) step = 1;
79 if (step > 26) step = 26;
80 return step;
81 }
82
geohashBitsComparator(const GeoHashBits * a,const GeoHashBits * b)83 int geohashBitsComparator(const GeoHashBits *a, const GeoHashBits *b) {
84 /* If step not equal, compare on step. Else, compare on bits. */
85 return a->step != b->step ? a->step - b->step : a->bits - b->bits;
86 }
87
geohashBoundingBox(double longitude,double latitude,double radius_meters,double * bounds)88 int geohashBoundingBox(double longitude, double latitude, double radius_meters,
89 double *bounds) {
90 if (!bounds) return 0;
91
92 double lonr, latr;
93 lonr = deg_rad(longitude);
94 latr = deg_rad(latitude);
95
96 if (radius_meters > EARTH_RADIUS_IN_METERS)
97 radius_meters = EARTH_RADIUS_IN_METERS;
98 double distance = radius_meters / EARTH_RADIUS_IN_METERS;
99 double min_latitude = latr - distance;
100 double max_latitude = latr + distance;
101
102 /* Note: we're being lazy and not accounting for coordinates near poles */
103 double min_longitude, max_longitude;
104 double difference_longitude = asin(sin(distance) / cos(latr));
105 min_longitude = lonr - difference_longitude;
106 max_longitude = lonr + difference_longitude;
107
108 bounds[0] = rad_deg(min_longitude);
109 bounds[1] = rad_deg(min_latitude);
110 bounds[2] = rad_deg(max_longitude);
111 bounds[3] = rad_deg(max_latitude);
112 return 1;
113 }
114
115 /* Return a set of areas (center + 8) that are able to cover a range query
116 * for the specified position and radius. */
geohashGetAreasByRadius(double longitude,double latitude,double radius_meters)117 GeoHashRadius geohashGetAreasByRadius(double longitude, double latitude, double radius_meters) {
118 GeoHashRange long_range, lat_range;
119 GeoHashRadius radius;
120 GeoHashBits hash;
121 GeoHashNeighbors neighbors;
122 GeoHashArea area;
123 double min_lon, max_lon, min_lat, max_lat;
124 double bounds[4];
125 int steps;
126
127 geohashBoundingBox(longitude, latitude, radius_meters, bounds);
128 min_lon = bounds[0];
129 min_lat = bounds[1];
130 max_lon = bounds[2];
131 max_lat = bounds[3];
132
133 steps = geohashEstimateStepsByRadius(radius_meters,latitude);
134
135 geohashGetCoordRange(&long_range,&lat_range);
136 geohashEncode(&long_range,&lat_range,longitude,latitude,steps,&hash);
137 geohashNeighbors(&hash,&neighbors);
138 geohashDecode(long_range,lat_range,hash,&area);
139
140 /* Check if the step is enough at the limits of the covered area.
141 * Sometimes when the search area is near an edge of the
142 * area, the estimated step is not small enough, since one of the
143 * north / south / west / east square is too near to the search area
144 * to cover everything. */
145 int decrease_step = 0;
146 {
147 GeoHashArea north, south, east, west;
148
149 geohashDecode(long_range, lat_range, neighbors.north, &north);
150 geohashDecode(long_range, lat_range, neighbors.south, &south);
151 geohashDecode(long_range, lat_range, neighbors.east, &east);
152 geohashDecode(long_range, lat_range, neighbors.west, &west);
153
154 if (geohashGetDistance(longitude,latitude,longitude,north.latitude.max)
155 < radius_meters) decrease_step = 1;
156 if (geohashGetDistance(longitude,latitude,longitude,south.latitude.min)
157 < radius_meters) decrease_step = 1;
158 if (geohashGetDistance(longitude,latitude,east.longitude.max,latitude)
159 < radius_meters) decrease_step = 1;
160 if (geohashGetDistance(longitude,latitude,west.longitude.min,latitude)
161 < radius_meters) decrease_step = 1;
162 }
163
164 if (decrease_step) {
165 steps--;
166 geohashEncode(&long_range,&lat_range,longitude,latitude,steps,&hash);
167 geohashNeighbors(&hash,&neighbors);
168 geohashDecode(long_range,lat_range,hash,&area);
169 }
170
171 /* Exclude the search areas that are useless. */
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 radius.hash = hash;
193 radius.neighbors = neighbors;
194 radius.area = area;
195 return radius;
196 }
197
geohashGetAreasByRadiusWGS84(double longitude,double latitude,double radius_meters)198 GeoHashRadius geohashGetAreasByRadiusWGS84(double longitude, double latitude,
199 double radius_meters) {
200 return geohashGetAreasByRadius(longitude, latitude, radius_meters);
201 }
202
geohashAlign52Bits(const GeoHashBits hash)203 GeoHashFix52Bits geohashAlign52Bits(const GeoHashBits hash) {
204 uint64_t bits = hash.bits;
205 bits <<= (52 - hash.step * 2);
206 return bits;
207 }
208
209 /* Calculate distance using haversin great circle distance formula. */
geohashGetDistance(double lon1d,double lat1d,double lon2d,double lat2d)210 double geohashGetDistance(double lon1d, double lat1d, double lon2d, double lat2d) {
211 double lat1r, lon1r, lat2r, lon2r, u, v;
212 lat1r = deg_rad(lat1d);
213 lon1r = deg_rad(lon1d);
214 lat2r = deg_rad(lat2d);
215 lon2r = deg_rad(lon2d);
216 u = sin((lat2r - lat1r) / 2);
217 v = sin((lon2r - lon1r) / 2);
218 return 2.0 * EARTH_RADIUS_IN_METERS *
219 asin(sqrt(u * u + cos(lat1r) * cos(lat2r) * v * v));
220 }
221
geohashGetDistanceIfInRadius(double x1,double y1,double x2,double y2,double radius,double * distance)222 int geohashGetDistanceIfInRadius(double x1, double y1,
223 double x2, double y2, double radius,
224 double *distance) {
225 *distance = geohashGetDistance(x1, y1, x2, y2);
226 if (*distance > radius) return 0;
227 return 1;
228 }
229
geohashGetDistanceIfInRadiusWGS84(double x1,double y1,double x2,double y2,double radius,double * distance)230 int geohashGetDistanceIfInRadiusWGS84(double x1, double y1, double x2,
231 double y2, double radius,
232 double *distance) {
233 return geohashGetDistanceIfInRadius(x1, y1, x2, y2, radius, distance);
234 }
235