/src/h3/src/h3lib/lib/latLng.c
Line | Count | Source (jump to first uncovered line) |
1 | | /* |
2 | | * Copyright 2016-2021 Uber Technologies, Inc. |
3 | | * |
4 | | * Licensed under the Apache License, Version 2.0 (the "License"); |
5 | | * you may not use this file except in compliance with the License. |
6 | | * You may obtain a copy of the License at |
7 | | * |
8 | | * http://www.apache.org/licenses/LICENSE-2.0 |
9 | | * |
10 | | * Unless required by applicable law or agreed to in writing, software |
11 | | * distributed under the License is distributed on an "AS IS" BASIS, |
12 | | * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
13 | | * See the License for the specific language governing permissions and |
14 | | * limitations under the License. |
15 | | */ |
16 | | /** @file latLng.c |
17 | | * @brief Functions for working with lat/lng coordinates. |
18 | | */ |
19 | | |
20 | | #include "latLng.h" |
21 | | |
22 | | #include <math.h> |
23 | | #include <stdbool.h> |
24 | | |
25 | | #include "constants.h" |
26 | | #include "h3Assert.h" |
27 | | #include "h3api.h" |
28 | | #include "mathExtensions.h" |
29 | | |
30 | | /** |
31 | | * Normalizes radians to a value between 0.0 and two PI. |
32 | | * |
33 | | * @param rads The input radians value. |
34 | | * @return The normalized radians value. |
35 | | */ |
36 | 159M | double _posAngleRads(double rads) { |
37 | 159M | double tmp = ((rads < 0.0L) ? rads + M_2PI : rads); |
38 | 159M | if (rads >= M_2PI) tmp -= M_2PI; |
39 | 159M | return tmp; |
40 | 159M | } |
41 | | |
42 | | /** |
43 | | * Determines if the components of two spherical coordinates are within some |
44 | | * threshold distance of each other. |
45 | | * |
46 | | * @param p1 The first spherical coordinates. |
47 | | * @param p2 The second spherical coordinates. |
48 | | * @param threshold The threshold distance. |
49 | | * @return Whether or not the two coordinates are within the threshold distance |
50 | | * of each other. |
51 | | */ |
52 | | bool geoAlmostEqualThreshold(const LatLng *p1, const LatLng *p2, |
53 | 0 | double threshold) { |
54 | 0 | return fabs(p1->lat - p2->lat) < threshold && |
55 | 0 | fabs(p1->lng - p2->lng) < threshold; |
56 | 0 | } |
57 | | |
58 | | /** |
59 | | * Determines if the components of two spherical coordinates are within our |
60 | | * standard epsilon distance of each other. |
61 | | * |
62 | | * @param p1 The first spherical coordinates. |
63 | | * @param p2 The second spherical coordinates. |
64 | | * @return Whether or not the two coordinates are within the epsilon distance |
65 | | * of each other. |
66 | | */ |
67 | 0 | bool geoAlmostEqual(const LatLng *p1, const LatLng *p2) { |
68 | 0 | return geoAlmostEqualThreshold(p1, p2, EPSILON_RAD); |
69 | 0 | } |
70 | | |
71 | | /** |
72 | | * Set the components of spherical coordinates in decimal degrees. |
73 | | * |
74 | | * @param p The spherical coordinates. |
75 | | * @param latDegs The desired latitude in decimal degrees. |
76 | | * @param lngDegs The desired longitude in decimal degrees. |
77 | | */ |
78 | 0 | void setGeoDegs(LatLng *p, double latDegs, double lngDegs) { |
79 | 0 | _setGeoRads(p, H3_EXPORT(degsToRads)(latDegs), |
80 | 0 | H3_EXPORT(degsToRads)(lngDegs)); |
81 | 0 | } |
82 | | |
83 | | /** |
84 | | * Set the components of spherical coordinates in radians. |
85 | | * |
86 | | * @param p The spherical coordinates. |
87 | | * @param latRads The desired latitude in decimal radians. |
88 | | * @param lngRads The desired longitude in decimal radians. |
89 | | */ |
90 | 0 | void _setGeoRads(LatLng *p, double latRads, double lngRads) { |
91 | 0 | p->lat = latRads; |
92 | 0 | p->lng = lngRads; |
93 | 0 | } |
94 | | |
95 | | /** |
96 | | * Convert from decimal degrees to radians. |
97 | | * |
98 | | * @param degrees The decimal degrees. |
99 | | * @return The corresponding radians. |
100 | | */ |
101 | 0 | double H3_EXPORT(degsToRads)(double degrees) { return degrees * M_PI_180; } |
102 | | |
103 | | /** |
104 | | * Convert from radians to decimal degrees. |
105 | | * |
106 | | * @param radians The radians. |
107 | | * @return The corresponding decimal degrees. |
108 | | */ |
109 | 0 | double H3_EXPORT(radsToDegs)(double radians) { return radians * M_180_PI; } |
110 | | |
111 | | /** |
112 | | * constrainLat makes sure latitudes are in the proper bounds |
113 | | * |
114 | | * @param lat The original lat value |
115 | | * @return The corrected lat value |
116 | | */ |
117 | 0 | double constrainLat(double lat) { |
118 | 0 | while (lat > M_PI_2) { |
119 | 0 | lat = lat - M_PI; |
120 | 0 | } |
121 | 0 | return lat; |
122 | 0 | } |
123 | | |
124 | | /** |
125 | | * constrainLng makes sure longitudes are in the proper bounds |
126 | | * |
127 | | * @param lng The origin lng value |
128 | | * @return The corrected lng value |
129 | | */ |
130 | 44.0M | double constrainLng(double lng) { |
131 | 45.8M | while (lng > M_PI) { |
132 | 1.83M | lng = lng - (2 * M_PI); |
133 | 1.83M | } |
134 | 44.6M | while (lng < -M_PI) { |
135 | 592k | lng = lng + (2 * M_PI); |
136 | 592k | } |
137 | 44.0M | return lng; |
138 | 44.0M | } |
139 | | |
140 | | /** |
141 | | * The great circle distance in radians between two spherical coordinates. |
142 | | * |
143 | | * This function uses the Haversine formula. |
144 | | * For math details, see: |
145 | | * https://en.wikipedia.org/wiki/Haversine_formula |
146 | | * https://www.movable-type.co.uk/scripts/latlong.html |
147 | | * |
148 | | * @param a the first lat/lng pair (in radians) |
149 | | * @param b the second lat/lng pair (in radians) |
150 | | * |
151 | | * @return the great circle distance in radians between a and b |
152 | | */ |
153 | 401k | double H3_EXPORT(greatCircleDistanceRads)(const LatLng *a, const LatLng *b) { |
154 | 401k | double sinLat = sin((b->lat - a->lat) / 2.0); |
155 | 401k | double sinLng = sin((b->lng - a->lng) / 2.0); |
156 | | |
157 | 401k | double A = sinLat * sinLat + cos(a->lat) * cos(b->lat) * sinLng * sinLng; |
158 | | |
159 | 401k | return 2 * atan2(sqrt(A), sqrt(1 - A)); |
160 | 401k | } |
161 | | |
162 | | /** |
163 | | * The great circle distance in kilometers between two spherical coordinates. |
164 | | */ |
165 | 401k | double H3_EXPORT(greatCircleDistanceKm)(const LatLng *a, const LatLng *b) { |
166 | 401k | return H3_EXPORT(greatCircleDistanceRads)(a, b) * EARTH_RADIUS_KM; |
167 | 401k | } |
168 | | |
169 | | /** |
170 | | * The great circle distance in meters between two spherical coordinates. |
171 | | */ |
172 | 0 | double H3_EXPORT(greatCircleDistanceM)(const LatLng *a, const LatLng *b) { |
173 | 0 | return H3_EXPORT(greatCircleDistanceKm)(a, b) * 1000; |
174 | 0 | } |
175 | | |
176 | | /** |
177 | | * Determines the azimuth to p2 from p1 in radians. |
178 | | * |
179 | | * @param p1 The first spherical coordinates. |
180 | | * @param p2 The second spherical coordinates. |
181 | | * @return The azimuth in radians from p1 to p2. |
182 | | */ |
183 | 13.9M | double _geoAzimuthRads(const LatLng *p1, const LatLng *p2) { |
184 | 13.9M | return atan2(cos(p2->lat) * sin(p2->lng - p1->lng), |
185 | 13.9M | cos(p1->lat) * sin(p2->lat) - |
186 | 13.9M | sin(p1->lat) * cos(p2->lat) * cos(p2->lng - p1->lng)); |
187 | 13.9M | } |
188 | | |
189 | | /** |
190 | | * Computes the point on the sphere a specified azimuth and distance from |
191 | | * another point. |
192 | | * |
193 | | * @param p1 The first spherical coordinates. |
194 | | * @param az The desired azimuth from p1. |
195 | | * @param distance The desired distance from p1, must be non-negative. |
196 | | * @param p2 The spherical coordinates at the desired azimuth and distance from |
197 | | * p1. |
198 | | */ |
199 | | void _geoAzDistanceRads(const LatLng *p1, double az, double distance, |
200 | 44.0M | LatLng *p2) { |
201 | 44.0M | if (distance < EPSILON) { |
202 | 0 | *p2 = *p1; |
203 | 0 | return; |
204 | 0 | } |
205 | | |
206 | 44.0M | double sinlat, sinlng, coslng; |
207 | | |
208 | 44.0M | az = _posAngleRads(az); |
209 | | |
210 | | // check for due north/south azimuth |
211 | 44.0M | if (az < EPSILON || fabs(az - M_PI) < EPSILON) { |
212 | 0 | if (az < EPSILON) // due north |
213 | 0 | p2->lat = p1->lat + distance; |
214 | 0 | else // due south |
215 | 0 | p2->lat = p1->lat - distance; |
216 | |
|
217 | 0 | if (fabs(p2->lat - M_PI_2) < EPSILON) // north pole |
218 | 0 | { |
219 | 0 | p2->lat = M_PI_2; |
220 | 0 | p2->lng = 0.0; |
221 | 0 | } else if (fabs(p2->lat + M_PI_2) < EPSILON) // south pole |
222 | 0 | { |
223 | 0 | p2->lat = -M_PI_2; |
224 | 0 | p2->lng = 0.0; |
225 | 0 | } else |
226 | 0 | p2->lng = constrainLng(p1->lng); |
227 | 0 | } else // not due north or south |
228 | 44.0M | { |
229 | 44.0M | sinlat = sin(p1->lat) * cos(distance) + |
230 | 44.0M | cos(p1->lat) * sin(distance) * cos(az); |
231 | 44.0M | if (sinlat > 1.0) sinlat = 1.0; |
232 | 44.0M | if (sinlat < -1.0) sinlat = -1.0; |
233 | 44.0M | p2->lat = asin(sinlat); |
234 | 44.0M | if (fabs(p2->lat - M_PI_2) < EPSILON) // north pole |
235 | 0 | { |
236 | 0 | p2->lat = M_PI_2; |
237 | 0 | p2->lng = 0.0; |
238 | 44.0M | } else if (fabs(p2->lat + M_PI_2) < EPSILON) // south pole |
239 | 0 | { |
240 | 0 | p2->lat = -M_PI_2; |
241 | 0 | p2->lng = 0.0; |
242 | 44.0M | } else { |
243 | 44.0M | sinlng = sin(az) * sin(distance) / cos(p2->lat); |
244 | 44.0M | coslng = (cos(distance) - sin(p1->lat) * sin(p2->lat)) / |
245 | 44.0M | cos(p1->lat) / cos(p2->lat); |
246 | 44.0M | if (sinlng > 1.0) sinlng = 1.0; |
247 | 44.0M | if (sinlng < -1.0) sinlng = -1.0; |
248 | 44.0M | if (coslng > 1.0) coslng = 1.0; |
249 | 44.0M | if (coslng < -1.0) coslng = -1.0; |
250 | 44.0M | p2->lng = constrainLng(p1->lng + atan2(sinlng, coslng)); |
251 | 44.0M | } |
252 | 44.0M | } |
253 | 44.0M | } |
254 | | |
255 | | /* |
256 | | * The following functions provide meta information about the H3 hexagons at |
257 | | * each zoom level. Since there are only 16 total levels, these are current |
258 | | * handled with hardwired static values, but it may be worthwhile to put these |
259 | | * static values into another file that can be autogenerated by source code in |
260 | | * the future. |
261 | | */ |
262 | | |
263 | 0 | H3Error H3_EXPORT(getHexagonAreaAvgKm2)(int res, double *out) { |
264 | 0 | static const double areas[] = { |
265 | 0 | 4.357449416078383e+06, 6.097884417941332e+05, 8.680178039899720e+04, |
266 | 0 | 1.239343465508816e+04, 1.770347654491307e+03, 2.529038581819449e+02, |
267 | 0 | 3.612906216441245e+01, 5.161293359717191e+00, 7.373275975944177e-01, |
268 | 0 | 1.053325134272067e-01, 1.504750190766435e-02, 2.149643129451879e-03, |
269 | 0 | 3.070918756316060e-04, 4.387026794728296e-05, 6.267181135324313e-06, |
270 | 0 | 8.953115907605790e-07}; |
271 | 0 | if (res < 0 || res > MAX_H3_RES) { |
272 | 0 | return E_RES_DOMAIN; |
273 | 0 | } |
274 | 0 | *out = areas[res]; |
275 | 0 | return E_SUCCESS; |
276 | 0 | } |
277 | | |
278 | 0 | H3Error H3_EXPORT(getHexagonAreaAvgM2)(int res, double *out) { |
279 | 0 | static const double areas[] = { |
280 | 0 | 4.357449416078390e+12, 6.097884417941339e+11, 8.680178039899731e+10, |
281 | 0 | 1.239343465508818e+10, 1.770347654491309e+09, 2.529038581819452e+08, |
282 | 0 | 3.612906216441250e+07, 5.161293359717198e+06, 7.373275975944188e+05, |
283 | 0 | 1.053325134272069e+05, 1.504750190766437e+04, 2.149643129451882e+03, |
284 | 0 | 3.070918756316063e+02, 4.387026794728301e+01, 6.267181135324322e+00, |
285 | 0 | 8.953115907605802e-01}; |
286 | 0 | if (res < 0 || res > MAX_H3_RES) { |
287 | 0 | return E_RES_DOMAIN; |
288 | 0 | } |
289 | 0 | *out = areas[res]; |
290 | 0 | return E_SUCCESS; |
291 | 0 | } |
292 | | |
293 | 0 | H3Error H3_EXPORT(getHexagonEdgeLengthAvgKm)(int res, double *out) { |
294 | 0 | static const double lens[] = { |
295 | 0 | 1281.256011, 483.0568391, 182.5129565, 68.97922179, |
296 | 0 | 26.07175968, 9.854090990, 3.724532667, 1.406475763, |
297 | 0 | 0.531414010, 0.200786148, 0.075863783, 0.028663897, |
298 | 0 | 0.010830188, 0.004092010, 0.001546100, 0.000584169}; |
299 | 0 | if (res < 0 || res > MAX_H3_RES) { |
300 | 0 | return E_RES_DOMAIN; |
301 | 0 | } |
302 | 0 | *out = lens[res]; |
303 | 0 | return E_SUCCESS; |
304 | 0 | } |
305 | | |
306 | 0 | H3Error H3_EXPORT(getHexagonEdgeLengthAvgM)(int res, double *out) { |
307 | 0 | static const double lens[] = { |
308 | 0 | 1281256.011, 483056.8391, 182512.9565, 68979.22179, |
309 | 0 | 26071.75968, 9854.090990, 3724.532667, 1406.475763, |
310 | 0 | 531.4140101, 200.7861476, 75.86378287, 28.66389748, |
311 | 0 | 10.83018784, 4.092010473, 1.546099657, 0.584168630}; |
312 | 0 | if (res < 0 || res > MAX_H3_RES) { |
313 | 0 | return E_RES_DOMAIN; |
314 | 0 | } |
315 | 0 | *out = lens[res]; |
316 | 0 | return E_SUCCESS; |
317 | 0 | } |
318 | | |
319 | 0 | H3Error H3_EXPORT(getNumCells)(int res, int64_t *out) { |
320 | 0 | if (res < 0 || res > MAX_H3_RES) { |
321 | 0 | return E_RES_DOMAIN; |
322 | 0 | } |
323 | 0 | *out = 2 + 120 * _ipow(7, res); |
324 | 0 | return E_SUCCESS; |
325 | 0 | } |
326 | | |
327 | | /** |
328 | | * Surface area in radians^2 of spherical triangle on unit sphere. |
329 | | * |
330 | | * For the math, see: |
331 | | * https://en.wikipedia.org/wiki/Spherical_trigonometry#Area_and_spherical_excess |
332 | | * |
333 | | * @param a length of triangle side A in radians |
334 | | * @param b length of triangle side B in radians |
335 | | * @param c length of triangle side C in radians |
336 | | * |
337 | | * @return area in radians^2 of triangle on unit sphere |
338 | | */ |
339 | 0 | double triangleEdgeLengthsToArea(double a, double b, double c) { |
340 | 0 | double s = (a + b + c) / 2; |
341 | |
|
342 | 0 | a = (s - a) / 2; |
343 | 0 | b = (s - b) / 2; |
344 | 0 | c = (s - c) / 2; |
345 | 0 | s = s / 2; |
346 | |
|
347 | 0 | return 4 * atan(sqrt(tan(s) * tan(a) * tan(b) * tan(c))); |
348 | 0 | } |
349 | | |
350 | | /** |
351 | | * Compute area in radians^2 of a spherical triangle, given its vertices. |
352 | | * |
353 | | * @param a vertex lat/lng in radians |
354 | | * @param b vertex lat/lng in radians |
355 | | * @param c vertex lat/lng in radians |
356 | | * |
357 | | * @return area of triangle on unit sphere, in radians^2 |
358 | | */ |
359 | 0 | double triangleArea(const LatLng *a, const LatLng *b, const LatLng *c) { |
360 | 0 | return triangleEdgeLengthsToArea(H3_EXPORT(greatCircleDistanceRads)(a, b), |
361 | 0 | H3_EXPORT(greatCircleDistanceRads)(b, c), |
362 | 0 | H3_EXPORT(greatCircleDistanceRads)(c, a)); |
363 | 0 | } |
364 | | |
365 | | /** |
366 | | * Area of H3 cell in radians^2. |
367 | | * |
368 | | * The area is calculated by breaking the cell into spherical triangles and |
369 | | * summing up their areas. Note that some H3 cells (hexagons and pentagons) |
370 | | * are irregular, and have more than 6 or 5 sides. |
371 | | * |
372 | | * todo: optimize the computation by re-using the edges shared between triangles |
373 | | * |
374 | | * @param cell H3 cell |
375 | | * @param out cell area in radians^2 |
376 | | * @return E_SUCCESS on success, or an error code otherwise |
377 | | */ |
378 | 0 | H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) { |
379 | 0 | LatLng c; |
380 | 0 | CellBoundary cb; |
381 | 0 | H3Error err = H3_EXPORT(cellToLatLng)(cell, &c); |
382 | 0 | if (err) { |
383 | 0 | return err; |
384 | 0 | } |
385 | 0 | err = H3_EXPORT(cellToBoundary)(cell, &cb); |
386 | 0 | if (NEVER(err)) { |
387 | | // Uncoverable because cellToLatLng will have returned an error already |
388 | 0 | return err; |
389 | 0 | } |
390 | | |
391 | 0 | double area = 0.0; |
392 | 0 | for (int i = 0; i < cb.numVerts; i++) { |
393 | 0 | int j = (i + 1) % cb.numVerts; |
394 | 0 | area += triangleArea(&cb.verts[i], &cb.verts[j], &c); |
395 | 0 | } |
396 | |
|
397 | 0 | *out = area; |
398 | 0 | return E_SUCCESS; |
399 | 0 | } |
400 | | |
401 | | /** |
402 | | * Area of H3 cell in kilometers^2. |
403 | | */ |
404 | 0 | H3Error H3_EXPORT(cellAreaKm2)(H3Index cell, double *out) { |
405 | 0 | H3Error err = H3_EXPORT(cellAreaRads2)(cell, out); |
406 | 0 | if (!err) { |
407 | 0 | *out = *out * EARTH_RADIUS_KM * EARTH_RADIUS_KM; |
408 | 0 | } |
409 | 0 | return err; |
410 | 0 | } |
411 | | |
412 | | /** |
413 | | * Area of H3 cell in meters^2. |
414 | | */ |
415 | 0 | H3Error H3_EXPORT(cellAreaM2)(H3Index cell, double *out) { |
416 | 0 | H3Error err = H3_EXPORT(cellAreaKm2)(cell, out); |
417 | 0 | if (!err) { |
418 | 0 | *out = *out * 1000 * 1000; |
419 | 0 | } |
420 | 0 | return err; |
421 | 0 | } |
422 | | |
423 | | /** |
424 | | * Length of a directed edge in radians. |
425 | | * |
426 | | * @param edge H3 directed edge |
427 | | * |
428 | | * @return length in radians |
429 | | */ |
430 | 0 | H3Error H3_EXPORT(edgeLengthRads)(H3Index edge, double *length) { |
431 | 0 | CellBoundary cb; |
432 | |
|
433 | 0 | H3Error err = H3_EXPORT(directedEdgeToBoundary)(edge, &cb); |
434 | 0 | if (err) { |
435 | 0 | return err; |
436 | 0 | } |
437 | | |
438 | 0 | *length = 0.0; |
439 | 0 | for (int i = 0; i < cb.numVerts - 1; i++) { |
440 | 0 | *length += |
441 | 0 | H3_EXPORT(greatCircleDistanceRads)(&cb.verts[i], &cb.verts[i + 1]); |
442 | 0 | } |
443 | |
|
444 | 0 | return E_SUCCESS; |
445 | 0 | } |
446 | | |
447 | | /** |
448 | | * Length of a directed edge in kilometers. |
449 | | */ |
450 | 0 | H3Error H3_EXPORT(edgeLengthKm)(H3Index edge, double *length) { |
451 | 0 | H3Error err = H3_EXPORT(edgeLengthRads)(edge, length); |
452 | 0 | *length = *length * EARTH_RADIUS_KM; |
453 | 0 | return err; |
454 | 0 | } |
455 | | |
456 | | /** |
457 | | * Length of a directed edge in meters. |
458 | | */ |
459 | 0 | H3Error H3_EXPORT(edgeLengthM)(H3Index edge, double *length) { |
460 | 0 | H3Error err = H3_EXPORT(edgeLengthKm)(edge, length); |
461 | 0 | *length = *length * 1000; |
462 | 0 | return err; |
463 | 0 | } |