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 "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
50const double DEG_TO_RAD = 0.017453292519943295769236907684886;
51/// @brief Earth's quatratic mean radius for WGS-84
52const double EARTH_RADIUS_IN_METERS = 6372797.560856;
53
54const double MERCATOR_MAX = 20037726.37;
55const double MERCATOR_MIN = -20037726.37;
56
57static inline double deg_rad(double ang) { return ang * D_R; }
58static 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. */
62uint8_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 */
98int 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 */
121GeoHashRadius 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 (north.latitude.max < max_lat)
168 decrease_step = 1;
169 if (south.latitude.min > min_lat)
170 decrease_step = 1;
171 if (east.longitude.max < max_lon)
172 decrease_step = 1;
173 if (west.longitude.min > min_lon)
174 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
213GeoHashFix52Bits 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. */
220double 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
232int 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
240int 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 */
254int 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