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 |
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 | |
57 | static inline double deg_rad(double ang) { return ang * D_R; } |
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. */ |
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 | */ |
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 */ |
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 (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 | |
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. */ |
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 | |
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 | |
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 | */ |
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 | |