1/*
2 * Copyright 2016-2017 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 geoCoord.c
17 * @brief Functions for working with lat/lon coordinates.
18 */
19
20#include "geoCoord.h"
21#include <math.h>
22#include <stdbool.h>
23#include "constants.h"
24#include "h3api.h"
25
26/**
27 * Normalizes radians to a value between 0.0 and two PI.
28 *
29 * @param rads The input radians value.
30 * @return The normalized radians value.
31 */
32double _posAngleRads(double rads) {
33 double tmp = ((rads < 0.0L) ? rads + M_2PI : rads);
34 if (rads >= M_2PI) tmp -= M_2PI;
35 return tmp;
36}
37
38/**
39 * Determines if the components of two spherical coordinates are within some
40 * threshold distance of each other.
41 *
42 * @param p1 The first spherical coordinates.
43 * @param p2 The second spherical coordinates.
44 * @param threshold The threshold distance.
45 * @return Whether or not the two coordinates are within the threshold distance
46 * of each other.
47 */
48bool geoAlmostEqualThreshold(const GeoCoord* p1, const GeoCoord* p2,
49 double threshold) {
50 return fabs(p1->lat - p2->lat) < threshold &&
51 fabs(p1->lon - p2->lon) < threshold;
52}
53
54/**
55 * Determines if the components of two spherical coordinates are within our
56 * standard epsilon distance of each other.
57 *
58 * @param p1 The first spherical coordinates.
59 * @param p2 The second spherical coordinates.
60 * @return Whether or not the two coordinates are within the epsilon distance
61 * of each other.
62 */
63bool geoAlmostEqual(const GeoCoord* p1, const GeoCoord* p2) {
64 return geoAlmostEqualThreshold(p1, p2, EPSILON_RAD);
65}
66
67/**
68 * Set the components of spherical coordinates in decimal degrees.
69 *
70 * @param p The spherical coodinates.
71 * @param latDegs The desired latitidue in decimal degrees.
72 * @param lonDegs The desired longitude in decimal degrees.
73 */
74void setGeoDegs(GeoCoord* p, double latDegs, double lonDegs) {
75 _setGeoRads(p, H3_EXPORT(degsToRads)(latDegs),
76 H3_EXPORT(degsToRads)(lonDegs));
77}
78
79/**
80 * Set the components of spherical coordinates in radians.
81 *
82 * @param p The spherical coodinates.
83 * @param latRads The desired latitidue in decimal radians.
84 * @param lonRads The desired longitude in decimal radians.
85 */
86void _setGeoRads(GeoCoord* p, double latRads, double lonRads) {
87 p->lat = latRads;
88 p->lon = lonRads;
89}
90
91/**
92 * Convert from decimal degrees to radians.
93 *
94 * @param degrees The decimal degrees.
95 * @return The corresponding radians.
96 */
97double H3_EXPORT(degsToRads)(double degrees) { return degrees * M_PI_180; }
98
99/**
100 * Convert from radians to decimal degrees.
101 *
102 * @param radians The radians.
103 * @return The corresponding decimal degrees.
104 */
105double H3_EXPORT(radsToDegs)(double radians) { return radians * M_180_PI; }
106
107/**
108 * constrainLat makes sure latitudes are in the proper bounds
109 *
110 * @param lat The original lat value
111 * @return The corrected lat value
112 */
113double constrainLat(double lat) {
114 while (lat > M_PI_2) {
115 lat = lat - M_PI;
116 }
117 return lat;
118}
119
120/**
121 * constrainLng makes sure longitudes are in the proper bounds
122 *
123 * @param lng The origin lng value
124 * @return The corrected lng value
125 */
126double constrainLng(double lng) {
127 while (lng > M_PI) {
128 lng = lng - (2 * M_PI);
129 }
130 while (lng < -M_PI) {
131 lng = lng + (2 * M_PI);
132 }
133 return lng;
134}
135
136/**
137 * Find the great circle distance in radians between two spherical coordinates.
138 *
139 * @param p1 The first spherical coordinates.
140 * @param p2 The second spherical coordinates.
141 * @return The great circle distance in radians between p1 and p2.
142 */
143double _geoDistRads(const GeoCoord* p1, const GeoCoord* p2) {
144 // use spherical triangle with p1 at A, p2 at B, and north pole at C
145 double bigC = fabs(p2->lon - p1->lon);
146 if (bigC > M_PI) // assume we want the complement
147 {
148 // note that in this case they can't both be negative
149 double lon1 = p1->lon;
150 if (lon1 < 0.0L) lon1 += 2.0L * M_PI;
151 double lon2 = p2->lon;
152 if (lon2 < 0.0L) lon2 += 2.0L * M_PI;
153
154 bigC = fabs(lon2 - lon1);
155 }
156
157 double b = M_PI_2 - p1->lat;
158 double a = M_PI_2 - p2->lat;
159
160 // use law of cosines to find c
161 double cosc = cos(a) * cos(b) + sin(a) * sin(b) * cos(bigC);
162 if (cosc > 1.0L) cosc = 1.0L;
163 if (cosc < -1.0L) cosc = -1.0L;
164
165 return acos(cosc);
166}
167
168/**
169 * Find the great circle distance in kilometers between two spherical
170 * coordinates.
171 *
172 * @param p1 The first spherical coordinates.
173 * @param p2 The second spherical coordinates.
174 * @return The distance in kilometers between p1 and p2.
175 */
176double _geoDistKm(const GeoCoord* p1, const GeoCoord* p2) {
177 return EARTH_RADIUS_KM * _geoDistRads(p1, p2);
178}
179
180/**
181 * Determines the azimuth to p2 from p1 in radians.
182 *
183 * @param p1 The first spherical coordinates.
184 * @param p2 The second spherical coordinates.
185 * @return The azimuth in radians from p1 to p2.
186 */
187double _geoAzimuthRads(const GeoCoord* p1, const GeoCoord* p2) {
188 return atan2(cos(p2->lat) * sin(p2->lon - p1->lon),
189 cos(p1->lat) * sin(p2->lat) -
190 sin(p1->lat) * cos(p2->lat) * cos(p2->lon - p1->lon));
191}
192
193/**
194 * Computes the point on the sphere a specified azimuth and distance from
195 * another point.
196 *
197 * @param p1 The first spherical coordinates.
198 * @param az The desired azimuth from p1.
199 * @param distance The desired distance from p1, must be non-negative.
200 * @param p2 The spherical coordinates at the desired azimuth and distance from
201 * p1.
202 */
203void _geoAzDistanceRads(const GeoCoord* p1, double az, double distance,
204 GeoCoord* p2) {
205 if (distance < EPSILON) {
206 *p2 = *p1;
207 return;
208 }
209
210 double sinlat, sinlon, coslon;
211
212 az = _posAngleRads(az);
213
214 // check for due north/south azimuth
215 if (az < EPSILON || fabs(az - M_PI) < EPSILON) {
216 if (az < EPSILON) // due north
217 p2->lat = p1->lat + distance;
218 else // due south
219 p2->lat = p1->lat - distance;
220
221 if (fabs(p2->lat - M_PI_2) < EPSILON) // north pole
222 {
223 p2->lat = M_PI_2;
224 p2->lon = 0.0L;
225 } else if (fabs(p2->lat + M_PI_2) < EPSILON) // south pole
226 {
227 p2->lat = -M_PI_2;
228 p2->lon = 0.0L;
229 } else
230 p2->lon = constrainLng(p1->lon);
231 } else // not due north or south
232 {
233 sinlat = sin(p1->lat) * cos(distance) +
234 cos(p1->lat) * sin(distance) * cos(az);
235 if (sinlat > 1.0L) sinlat = 1.0L;
236 if (sinlat < -1.0L) sinlat = -1.0L;
237 p2->lat = asin(sinlat);
238 if (fabs(p2->lat - M_PI_2) < EPSILON) // north pole
239 {
240 p2->lat = M_PI_2;
241 p2->lon = 0.0L;
242 } else if (fabs(p2->lat + M_PI_2) < EPSILON) // south pole
243 {
244 p2->lat = -M_PI_2;
245 p2->lon = 0.0L;
246 } else {
247 sinlon = sin(az) * sin(distance) / cos(p2->lat);
248 coslon = (cos(distance) - sin(p1->lat) * sin(p2->lat)) /
249 cos(p1->lat) / cos(p2->lat);
250 if (sinlon > 1.0L) sinlon = 1.0L;
251 if (sinlon < -1.0L) sinlon = -1.0L;
252 if (coslon > 1.0L) sinlon = 1.0L;
253 if (coslon < -1.0L) sinlon = -1.0L;
254 p2->lon = constrainLng(p1->lon + atan2(sinlon, coslon));
255 }
256 }
257}
258
259/*
260 * The following functions provide meta information about the H3 hexagons at
261 * each zoom level. Since there are only 16 total levels, these are current
262 * handled with hardwired static values, but it may be worthwhile to put these
263 * static values into another file that can be autogenerated by source code in
264 * the future.
265 */
266
267double H3_EXPORT(hexAreaKm2)(int res) {
268 static const double areas[] = {
269 4250546.848, 607220.9782, 86745.85403, 12392.26486,
270 1770.323552, 252.9033645, 36.1290521, 5.1612932,
271 0.7373276, 0.1053325, 0.0150475, 0.0021496,
272 0.0003071, 0.0000439, 0.0000063, 0.0000009};
273 return areas[res];
274}
275
276double H3_EXPORT(hexAreaM2)(int res) {
277 static const double areas[] = {
278 4.25055E+12, 6.07221E+11, 86745854035, 12392264862,
279 1770323552, 252903364.5, 36129052.1, 5161293.2,
280 737327.6, 105332.5, 15047.5, 2149.6,
281 307.1, 43.9, 6.3, 0.9};
282 return areas[res];
283}
284
285double H3_EXPORT(edgeLengthKm)(int res) {
286 static const double lens[] = {
287 1107.712591, 418.6760055, 158.2446558, 59.81085794,
288 22.6063794, 8.544408276, 3.229482772, 1.220629759,
289 0.461354684, 0.174375668, 0.065907807, 0.024910561,
290 0.009415526, 0.003559893, 0.001348575, 0.000509713};
291 return lens[res];
292}
293
294double H3_EXPORT(edgeLengthM)(int res) {
295 static const double lens[] = {
296 1107712.591, 418676.0055, 158244.6558, 59810.85794,
297 22606.3794, 8544.408276, 3229.482772, 1220.629759,
298 461.3546837, 174.3756681, 65.90780749, 24.9105614,
299 9.415526211, 3.559893033, 1.348574562, 0.509713273};
300 return lens[res];
301}
302
303/** @brief Number of unique valid H3Indexes at given resolution. */
304int64_t H3_EXPORT(numHexagons)(int res) {
305 static const int64_t nums[] = {122L,
306 842L,
307 5882L,
308 41162L,
309 288122L,
310 2016842L,
311 14117882L,
312 98825162L,
313 691776122L,
314 4842432842L,
315 33897029882L,
316 237279209162L,
317 1660954464122L,
318 11626681248842L,
319 81386768741882L,
320 569707381193162L};
321 return nums[res];
322}
323