OmniSciDB  a5dc49c757
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
geoCoord.hpp
Go to the documentation of this file.
1 /*
2  * Copyright 2016-2020 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  */
21 
22 // #include <math.h>
23 // #include <stdbool.h>
24 
25 // #include "constants.h"
26 // #include "h3api.h"
27 // #include "mathExtensions.h"
28 
35 EXTENSION_NOINLINE double _posAngleRads(double rads) {
36  double tmp = ((rads < 0.0) ? rads + M_2PI : rads);
37  if (rads >= M_2PI)
38  tmp -= M_2PI;
39  return tmp;
40 }
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 GeoCoord *p1, const GeoCoord *p2,
53 // double threshold) {
54 // return fabs(p1->lat - p2->lat) < threshold &&
55 // fabs(p1->lon - p2->lon) < threshold;
56 // }
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 // bool geoAlmostEqual(const GeoCoord *p1, const GeoCoord *p2) {
68 // return geoAlmostEqualThreshold(p1, p2, EPSILON_RAD);
69 // }
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 lonDegs The desired longitude in decimal degrees.
77 // */
78 // void setGeoDegs(GeoCoord *p, double latDegs, double lonDegs) {
79 // _setGeoRads(p, H3_EXPORT(degsToRads)(latDegs),
80 // H3_EXPORT(degsToRads)(lonDegs));
81 // }
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 lonRads The desired longitude in decimal radians.
89 // */
90 // void _setGeoRads(GeoCoord *p, double latRads, double lonRads) {
91 // p->lat = latRads;
92 // p->lon = lonRads;
93 // }
94 
102  return degrees * M_PI_180;
103 }
104 
112  return radians * M_180_PI;
113 }
114 
115 // /**
116 // * constrainLat makes sure latitudes are in the proper bounds
117 // *
118 // * @param lat The original lat value
119 // * @return The corrected lat value
120 // */
121 // double constrainLat(double lat) {
122 // while (lat > M_PI_2) {
123 // lat = lat - M_PI;
124 // }
125 // return lat;
126 // }
127 
134 EXTENSION_INLINE double constrainLng(double lng) {
135  while (lng > M_PI) {
136  lng = lng - (2 * M_PI);
137  }
138  while (lng < -M_PI) {
139  lng = lng + (2 * M_PI);
140  }
141  return lng;
142 }
143 
144 // /**
145 // * The great circle distance in radians between two spherical coordinates.
146 // *
147 // * This function uses the Haversine formula.
148 // * For math details, see:
149 // * https://en.wikipedia.org/wiki/Haversine_formula
150 // * https://www.movable-type.co.uk/scripts/latlong.html
151 // *
152 // * @param a the first lat/lng pair (in radians)
153 // * @param b the second lat/lng pair (in radians)
154 // *
155 // * @return the great circle distance in radians between a and b
156 // */
157 // double H3_EXPORT(pointDistRads)(const GeoCoord *a, const GeoCoord *b) {
158 // double sinLat = sin((b->lat - a->lat) / 2.0);
159 // double sinLng = sin((b->lon - a->lon) / 2.0);
160 
161 // double A = sinLat * sinLat + cos(a->lat) * cos(b->lat) * sinLng * sinLng;
162 
163 // return 2 * atan2(sqrt(A), sqrt(1 - A));
164 // }
165 
166 // /**
167 // * The great circle distance in kilometers between two spherical coordinates.
168 // */
169 // double H3_EXPORT(pointDistKm)(const GeoCoord *a, const GeoCoord *b) {
170 // return H3_EXPORT(pointDistRads)(a, b) * EARTH_RADIUS_KM;
171 // }
172 
173 // /**
174 // * The great circle distance in meters between two spherical coordinates.
175 // */
176 // double H3_EXPORT(pointDistM)(const GeoCoord *a, const GeoCoord *b) {
177 // return H3_EXPORT(pointDistKm)(a, b) * 1000;
178 // }
179 
187 EXTENSION_NOINLINE double _geoAzimuthRads(const GeoCoord(p1), const GeoCoord(p2)) {
188  return atan2(
189  cos(p2[LAT_INDEX]) * sin(p2[LON_INDEX] - p1[LON_INDEX]),
190  cos(p1[LAT_INDEX]) * sin(p2[LAT_INDEX]) -
191  sin(p1[LAT_INDEX]) * cos(p2[LAT_INDEX]) * cos(p2[LON_INDEX] - p1[LON_INDEX]));
192 }
193 
205  double az,
206  double distance,
207  GeoCoord(p2)) {
208  if (distance < EPSILON) {
209  GeoCoordCopy(p2, p1);
210  return true;
211  }
212 
213  double sinlat, sinlon, coslon;
214 
215  az = _posAngleRads(az);
216 
217  // check for due north/south azimuth
218  if (az < EPSILON || fabs(az - M_PI) < EPSILON) {
219  if (az < EPSILON) // due north
220  p2[LAT_INDEX] = p1[LAT_INDEX] + distance;
221  else // due south
222  p2[LAT_INDEX] = p1[LAT_INDEX] - distance;
223 
224  if (fabs(p2[LAT_INDEX] - M_PI_2) < EPSILON) // north pole
225  {
226  p2[LAT_INDEX] = M_PI_2;
227  p2[LON_INDEX] = 0.0;
228  } else if (fabs(p2[LAT_INDEX] + M_PI_2) < EPSILON) // south pole
229  {
230  p2[LAT_INDEX] = -M_PI_2;
231  p2[LON_INDEX] = 0.0;
232  } else
233  p2[LON_INDEX] = constrainLng(p1[LON_INDEX]);
234  } else // not due north or south
235  {
236  sinlat =
237  sin(p1[LAT_INDEX]) * cos(distance) + cos(p1[LAT_INDEX]) * sin(distance) * cos(az);
238  if (sinlat > 1.0)
239  sinlat = 1.0;
240  if (sinlat < -1.0)
241  sinlat = -1.0;
242  p2[LAT_INDEX] = asin(sinlat);
243  if (fabs(p2[LAT_INDEX] - M_PI_2) < EPSILON) // north pole
244  {
245  p2[LAT_INDEX] = M_PI_2;
246  p2[LON_INDEX] = 0.0;
247  } else if (fabs(p2[LAT_INDEX] + M_PI_2) < EPSILON) // south pole
248  {
249  p2[LAT_INDEX] = -M_PI_2;
250  p2[LON_INDEX] = 0.0;
251  } else {
252  sinlon = sin(az) * sin(distance) / cos(p2[LAT_INDEX]);
253  coslon = (cos(distance) - sin(p1[LAT_INDEX]) * sin(p2[LAT_INDEX])) /
254  cos(p1[LAT_INDEX]) / cos(p2[LAT_INDEX]);
255  if (sinlon > 1.0)
256  sinlon = 1.0;
257  if (sinlon < -1.0)
258  sinlon = -1.0;
259  if (coslon > 1.0)
260  coslon = 1.0;
261  if (coslon < -1.0)
262  coslon = -1.0;
263  p2[LON_INDEX] = constrainLng(p1[LON_INDEX] + atan2(sinlon, coslon));
264  }
265  }
266  return true;
267 }
268 
269 // /*
270 // * The following functions provide meta information about the H3 hexagons at
271 // * each zoom level. Since there are only 16 total levels, these are current
272 // * handled with hardwired static values, but it may be worthwhile to put these
273 // * static values into another file that can be autogenerated by source code in
274 // * the future.
275 // */
276 
277 // double H3_EXPORT(hexAreaKm2)(int res) {
278 // static const double areas[] = {
279 // 4250546.848, 607220.9782, 86745.85403, 12392.26486,
280 // 1770.323552, 252.9033645, 36.1290521, 5.1612932,
281 // 0.7373276, 0.1053325, 0.0150475, 0.0021496,
282 // 0.0003071, 0.0000439, 0.0000063, 0.0000009};
283 // return areas[res];
284 // }
285 
286 // double H3_EXPORT(hexAreaM2)(int res) {
287 // static const double areas[] = {
288 // 4.25055E+12, 6.07221E+11, 86745854035, 12392264862,
289 // 1770323552, 252903364.5, 36129052.1, 5161293.2,
290 // 737327.6, 105332.5, 15047.5, 2149.6,
291 // 307.1, 43.9, 6.3, 0.9};
292 // return areas[res];
293 // }
294 
295 // double H3_EXPORT(edgeLengthKm)(int res) {
296 // static const double lens[] = {
297 // 1107.712591, 418.6760055, 158.2446558, 59.81085794,
298 // 22.6063794, 8.544408276, 3.229482772, 1.220629759,
299 // 0.461354684, 0.174375668, 0.065907807, 0.024910561,
300 // 0.009415526, 0.003559893, 0.001348575, 0.000509713};
301 // return lens[res];
302 // }
303 
304 // double H3_EXPORT(edgeLengthM)(int res) {
305 // static const double lens[] = {
306 // 1107712.591, 418676.0055, 158244.6558, 59810.85794,
307 // 22606.3794, 8544.408276, 3229.482772, 1220.629759,
308 // 461.3546837, 174.3756681, 65.90780749, 24.9105614,
309 // 9.415526211, 3.559893033, 1.348574562, 0.509713273};
310 // return lens[res];
311 // }
312 
313 // int64_t H3_EXPORT(numHexagons)(int res) { return 2 + 120 * _ipow(7, res); }
314 
315 // /**
316 // * Surface area in radians^2 of spherical triangle on unit sphere.
317 // *
318 // * For the math, see:
319 // * https://en.wikipedia.org/wiki/Spherical_trigonometry#Area_and_spherical_excess
320 // *
321 // * @param a length of triangle side A in radians
322 // * @param b length of triangle side B in radians
323 // * @param c length of triangle side C in radians
324 // *
325 // * @return area in radians^2 of triangle on unit sphere
326 // */
327 // double triangleEdgeLengthsToArea(double a, double b, double c) {
328 // double s = (a + b + c) / 2;
329 
330 // a = (s - a) / 2;
331 // b = (s - b) / 2;
332 // c = (s - c) / 2;
333 // s = s / 2;
334 
335 // return 4 * atan(sqrt(tan(s) * tan(a) * tan(b) * tan(c)));
336 // }
337 
338 // /**
339 // * Compute area in radians^2 of a spherical triangle, given its vertices.
340 // *
341 // * @param a vertex lat/lng in radians
342 // * @param b vertex lat/lng in radians
343 // * @param c vertex lat/lng in radians
344 // *
345 // * @return area of triangle on unit sphere, in radians^2
346 // */
347 // double triangleArea(const GeoCoord *a, const GeoCoord *b, const GeoCoord *c) {
348 // return triangleEdgeLengthsToArea(H3_EXPORT(pointDistRads)(a, b),
349 // H3_EXPORT(pointDistRads)(b, c),
350 // H3_EXPORT(pointDistRads)(c, a));
351 // }
352 
353 // /**
354 // * Area of H3 cell in radians^2.
355 // *
356 // * The area is calculated by breaking the cell into spherical triangles and
357 // * summing up their areas. Note that some H3 cells (hexagons and pentagons)
358 // * are irregular, and have more than 6 or 5 sides.
359 // *
360 // * todo: optimize the computation by re-using the edges shared between triangles
361 // *
362 // * @param cell H3 cell
363 // *
364 // * @return cell area in radians^2
365 // */
366 // double H3_EXPORT(cellAreaRads2)(H3Index cell) {
367 // GeoCoord c;
368 // GeoBoundary gb;
369 // H3_EXPORT(h3ToGeo)(cell, &c);
370 // H3_EXPORT(h3ToGeoBoundary)(cell, &gb);
371 
372 // double area = 0.0;
373 // for (int i = 0; i < gb.numVerts; i++) {
374 // int j = (i + 1) % gb.numVerts;
375 // area += triangleArea(&gb.verts[i], &gb.verts[j], &c);
376 // }
377 
378 // return area;
379 // }
380 
381 // /**
382 // * Area of H3 cell in kilometers^2.
383 // */
384 // double H3_EXPORT(cellAreaKm2)(H3Index h) {
385 // return H3_EXPORT(cellAreaRads2)(h) * EARTH_RADIUS_KM * EARTH_RADIUS_KM;
386 // }
387 
388 // /**
389 // * Area of H3 cell in meters^2.
390 // */
391 // double H3_EXPORT(cellAreaM2)(H3Index h) {
392 // return H3_EXPORT(cellAreaKm2)(h) * 1000 * 1000;
393 // }
394 
395 // /**
396 // * Length of a unidirectional edge in radians.
397 // *
398 // * @param edge H3 unidirectional edge
399 // *
400 // * @return length in radians
401 // */
402 // double H3_EXPORT(exactEdgeLengthRads)(H3Index edge) {
403 // GeoBoundary gb;
404 
405 // H3_EXPORT(getH3UnidirectionalEdgeBoundary)(edge, &gb);
406 
407 // double length = 0.0;
408 // for (int i = 0; i < gb.numVerts - 1; i++) {
409 // length += H3_EXPORT(pointDistRads)(&gb.verts[i], &gb.verts[i + 1]);
410 // }
411 
412 // return length;
413 // }
414 
415 // /**
416 // * Length of a unidirectional edge in kilometers.
417 // */
418 // double H3_EXPORT(exactEdgeLengthKm)(H3Index edge) {
419 // return H3_EXPORT(exactEdgeLengthRads)(edge) * EARTH_RADIUS_KM;
420 // }
421 
422 // /**
423 // * Length of a unidirectional edge in meters.
424 // */
425 // double H3_EXPORT(exactEdgeLengthM)(H3Index edge) {
426 // return H3_EXPORT(exactEdgeLengthKm)(edge) * 1000;
427 // }
#define M_PI_180
Definition: constants.h:37
#define GeoCoord(variable_name)
Definition: h3api.h:94
#define EXTENSION_NOINLINE
Definition: heavydbTypes.h:58
#define H3_EXPORT(name)
Definition: h3api.h:43
EXTENSION_INLINE double H3_EXPORT() radsToDegs(double radians)
converts radians to degrees
Definition: geoCoord.hpp:111
EXTENSION_INLINE double constrainLng(double lng)
Definition: geoCoord.hpp:134
#define M_PI
Definition: constants.h:25
#define LAT_INDEX
Definition: h3api.h:92
EXTENSION_NOINLINE double _posAngleRads(double rads)
Definition: geoCoord.hpp:35
#define M_2PI
Definition: constants.h:34
#define EXTENSION_INLINE
Definition: heavydbTypes.h:57
EXTENSION_NOINLINE double radians(const double x)
Geodetic (lat/lon) functions.
EXTENSION_NOINLINE bool _geoAzDistanceRads(const GeoCoord(p1), double az, double distance, GeoCoord(p2))
Definition: geoCoord.hpp:204
#define M_PI_2
Definition: constants.h:30
EXTENSION_NOINLINE double degrees(double x)
#define EPSILON
Definition: constants.h:42
EXTENSION_NOINLINE double _geoAzimuthRads(const GeoCoord(p1), const GeoCoord(p2))
Definition: geoCoord.hpp:187
#define LON_INDEX
Definition: h3api.h:93
EXTENSION_INLINE double H3_EXPORT() degsToRads(double degrees)
converts degrees to radians
Definition: geoCoord.hpp:101
#define M_180_PI
Definition: constants.h:39
#define GeoCoordCopy(dest_coord, src_coord)
Definition: h3api.h:96