diff --git a/CHANGELOG.md b/CHANGELOG.md index 8476c8f78b..73f71e5dee 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,7 +7,11 @@ The public API of this library consists of the functions declared in file ## [Unreleased] ### Added -- `reverseDirectedEdge` function +- `reverseDirectedEdge` function (#1098) +- (internal) `geoLoopArea` function (#1101) + +### Changed +- `cellAreaRads2` uses `geoLoopArea` (#1101) ## [4.4.1] - 2025-11-11 ### Fixed diff --git a/CMakeLists.txt b/CMakeLists.txt index 7eb6dbddb5..6b4c8983ec 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -170,6 +170,8 @@ set(LIB_SOURCE_FILES src/h3lib/include/constants.h src/h3lib/include/coordijk.h src/h3lib/include/algos.h + src/h3lib/include/adder.h + src/h3lib/include/area.h src/h3lib/lib/h3Assert.c src/h3lib/lib/algos.c src/h3lib/lib/coordijk.c @@ -188,7 +190,8 @@ set(LIB_SOURCE_FILES src/h3lib/lib/iterators.c src/h3lib/lib/vertexGraph.c src/h3lib/lib/faceijk.c - src/h3lib/lib/baseCells.c) + src/h3lib/lib/baseCells.c + src/h3lib/lib/area.c) set(APP_SOURCE_FILES src/apps/applib/include/kml.h src/apps/applib/include/benchmark.h @@ -278,6 +281,7 @@ set(OTHER_SOURCE_FILES src/apps/testapps/testH3IteratorsInternal.c src/apps/testapps/testMathExtensionsInternal.c src/apps/testapps/testDescribeH3Error.c + src/apps/testapps/testGeoLoopArea.c src/apps/miscapps/cellToBoundaryHier.c src/apps/miscapps/cellToLatLngHier.c src/apps/miscapps/generateBaseCellNeighbors.c @@ -317,7 +321,8 @@ set(OTHER_SOURCE_FILES src/apps/benchmarks/benchmarkDirectedEdge.c src/apps/benchmarks/benchmarkVertex.c src/apps/benchmarks/benchmarkIsValidCell.c - src/apps/benchmarks/benchmarkH3Api.c) + src/apps/benchmarks/benchmarkH3Api.c + src/apps/benchmarks/benchmarkArea.c) set(ALL_SOURCE_FILES ${LIB_SOURCE_FILES} ${APP_SOURCE_FILES} ${TEST_APP_SOURCE_FILES} @@ -672,6 +677,8 @@ if(BUILD_BENCHMARKS) src/apps/benchmarks/benchmarkPolygonToCells.c) add_h3_benchmark(benchmarkPolygonToCellsExperimental src/apps/benchmarks/benchmarkPolygonToCellsExperimental.c) + add_h3_benchmark(benchmarkArea + src/apps/benchmarks/benchmarkArea.c) if(ENABLE_REQUIRES_ALL_SYMBOLS) add_h3_benchmark(benchmarkPolygon src/apps/benchmarks/benchmarkPolygon.c) diff --git a/CMakeTests.cmake b/CMakeTests.cmake index 1901a61fac..3508ff23a3 100644 --- a/CMakeTests.cmake +++ b/CMakeTests.cmake @@ -259,6 +259,7 @@ add_h3_test(testH3IteratorsInternal src/apps/testapps/testH3IteratorsInternal.c) add_h3_test(testMathExtensionsInternal src/apps/testapps/testMathExtensionsInternal.c) add_h3_test(testDescribeH3Error src/apps/testapps/testDescribeH3Error.c) +add_h3_test(testGeoLoopArea src/apps/testapps/testGeoLoopArea.c) add_h3_test_with_arg(testH3NeighborRotations src/apps/testapps/testH3NeighborRotations.c 0) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c new file mode 100644 index 0000000000..fa1b29e75e --- /dev/null +++ b/src/apps/benchmarks/benchmarkArea.c @@ -0,0 +1,59 @@ +/* + * Copyright 2025 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +#include +#include +#include + +#include "adder.h" +#include "benchmark.h" +#include "constants.h" +#include "h3api.h" +#include "iterators.h" + +static void doResSum(int res, bool print) { + Adder adder = {}; + double cellArea; + IterCellsResolution iter = iterInitRes(res); + + for (; iter.h; iterStepRes(&iter)) { + H3_EXPORT(cellAreaRads2)(iter.h, &cellArea); + kadd(&adder, cellArea); + } + + double total_area = adder.sum; + double diff = fabs(total_area - 4 * M_PI); + if (print) { + printf("res: %d, diff: %e\n", res, diff); + } +} + +BEGIN_BENCHMARKS(); + +static int MAX_RES = 3; + +BENCHMARK(allCellsAtRes_print, 1, { + for (int i = 0; i <= MAX_RES; i++) { + doResSum(i, true); + } +}); + +BENCHMARK(allCellsAtRes_noprint, 10, { + for (int i = 0; i <= MAX_RES; i++) { + doResSum(i, false); + } +}); + +END_BENCHMARKS(); diff --git a/src/apps/testapps/testGeoLoopArea.c b/src/apps/testapps/testGeoLoopArea.c new file mode 100644 index 0000000000..49643ca245 --- /dev/null +++ b/src/apps/testapps/testGeoLoopArea.c @@ -0,0 +1,228 @@ +/* + * Copyright 2025 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +/** @file + * @brief tests H3 GeoLoop area calculation + * + * usage: `testGeoLoopArea` + */ + +#include + +#include "area.h" +#include "h3api.h" +#include "test.h" +#include "utility.h" + +static void _compareArea(LatLng *verts, int numVerts, double target_area) { + double tol = 1e-14; + GeoLoop loop = {.verts = verts, .numVerts = numVerts}; + + double out; + t_assertSuccess(geoLoopAreaRads2(loop, &out)); + t_assert(fabs(out - target_area) < tol, "area should match"); +} + +SUITE(geoLoopArea) { + TEST(triangle_basic) { + /* + GeoLoop representing a triangle covering 1/8 of the globe, with + points ordered according to right-hand rule (counter-clockwise). + + The triangle starts at the north pole, moves down 90 degrees to the + equator, and then sweeps out 90 degrees along the equator + before returning to the north pole. + + The globe has an area of 4*pi radians, so this 1/8 triangle piece of + the globe should have area pi/2. + */ + LatLng verts[] = { + {M_PI_2, 0.0}, + {0.0, 0.0}, + {0.0, M_PI_2}, + }; + + _compareArea(verts, ARRAY_SIZE(verts), M_PI / 2); + } + + TEST(triangle_reversed) { + /* + Reverse the order of the points in the triangle from the previous test, + so that they are in clockwise order. + + Since the points are in clockwise order, GeoLoop represents the whole + globe minus the triangle above. + */ + LatLng verts[] = { + {0.0, M_PI_2}, + {0.0, 0.0}, + {M_PI_2, 0.0}, + }; + + _compareArea(verts, ARRAY_SIZE(verts), 7 * M_PI / 2); + } + + TEST(slice) { + /* + Stitch two 1/8 triangles together, sharing an edge along the equator + to create a 1/4 slice of the globe, with vertices at the north and south + pole. + */ + LatLng verts[] = { + {M_PI_2, 0.0}, + {0.0, 0.0}, + {-M_PI_2, 0.0}, + {0.0, M_PI_2}, + }; + + _compareArea(verts, ARRAY_SIZE(verts), M_PI); + } + + TEST(slice_reversed) { + /* + 3/4 slice of the globe, from north to south pole, formed by reversing + order of points from example above. + */ + LatLng verts[] = { + {M_PI_2, 0.0}, + {0.0, M_PI_2}, + {-M_PI_2, 0.0}, + {0.0, 0.0}, + }; + + _compareArea(verts, ARRAY_SIZE(verts), 3 * M_PI); + } + + TEST(hemisphere_east) { + /* + Stitch two 1/4 triangles together to cover the eastern hemisphere. + */ + LatLng verts[] = { + {M_PI_2, 0.0}, + {0.0, 0.0}, + {-M_PI_2, 0.0}, + {0.0, M_PI}, + }; + + _compareArea(verts, ARRAY_SIZE(verts), 2 * M_PI); + } + + TEST(hemisphere_north) { + /* + Stitch 4 1/8 triangles together to cover the northern hemisphere. + */ + LatLng verts[] = { + {0.0, -M_PI}, + {0.0, -M_PI_2}, + {0.0, 0.0}, + {0.0, M_PI_2}, + }; + + _compareArea(verts, ARRAY_SIZE(verts), 2 * M_PI); + } + + TEST(percentageSlice) { + /* + Demonstrate that edge arcs between points in a polygon or geoloop + should be less than 180 degrees (M_PI radians). + + Create a triangle from north pole to equator and back to the north pole + that sweeps out an edge arc of t * M_PI radians along the equator, + so it should have an area of t*M_PI for t \in [0,1]. + + However, there is a discontinuity at t = 1 (i.e., M_PI radians or 180 + degrees), where expected area goes to (2 + t) * M_PI for 1 < t < 2. + + Recall that the area in stradians of the entire globe is 4*M_PI. + */ + for (double t = 0; t <= 1.2; t += 0.01) { + LatLng verts[] = { + {M_PI_2, 0.0}, + {0.0, -M_PI_2}, + {0.0, t * M_PI - M_PI_2}, + }; + GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; + + double out; + t_assertSuccess(geoLoopAreaRads2(loop, &out)); + + double tol = 1e-13; + if (t < 0.99) { + // When t < 1, the largest angle in the triangle is less than + // 180 degrees + t_assert(fabs(out - t * M_PI) <= tol, "expected area"); + } else if (t > 1.01) { + /* + Discontinuity at t == 1. For t > 1, the triangle "flips", + because the shortest geodesic path is on the other side of + the globe. The triangle is now oriented in clockwise order, + and the area computed is the area *outside* of the triangle, + which starts at 3*pi. + */ + t_assert(fabs(out - (2 + t) * M_PI) <= tol, "expected area"); + } + /* + Note that we avoid testing t == 1, since the triangle + isn't well defined because there are many possible geodesic + shortest paths when consecutive points are antipodal + (180 degrees apart). + */ + } + } + + TEST(percentageSlice_large) { + /* + Continuing from the test above, note that a large polygon with + t > 1 is *still* representable and we can compute its area accurately; + we just need to add intermediate vertices so that + no edge arc is greater than 180 degrees. + */ + double t = 1.2; + LatLng verts[] = { + {M_PI_2, 0.0}, + {0.0, -M_PI_2}, + {0.0, 0.0}, // Extra vertex so every angle is < 180 degrees + {0.0, t * M_PI - M_PI_2}, + }; + + _compareArea(verts, ARRAY_SIZE(verts), t * M_PI); + } + + TEST(degenerateLoop2) { + // Note that `geoLoopAreaRads2()` works without error on degenerate + // loops, returning 0 area. + LatLng verts[] = { + {M_PI_2, 0.0}, + {0.0, -M_PI_2}, + }; + _compareArea(verts, ARRAY_SIZE(verts), 0.0); + } + + TEST(degenerateLoop1) { + // Note that `geoLoopAreaRads2()` works without error on degenerate + // loops, returning 0 area. + LatLng verts[] = { + {0.0, 0.0}, + }; + _compareArea(verts, ARRAY_SIZE(verts), 0.0); + } + + TEST(degenerateLoop0) { + // Note that `geoLoopAreaRads2()` works without error on degenerate + // loops, returning 0 area. + _compareArea(NULL, 0, 0.0); + } +} diff --git a/src/apps/testapps/testH3CellArea.c b/src/apps/testapps/testH3CellArea.c index 060e14d169..4189dfc7de 100644 --- a/src/apps/testapps/testH3CellArea.c +++ b/src/apps/testapps/testH3CellArea.c @@ -22,7 +22,7 @@ #include -#include "h3Index.h" +#include "h3api.h" #include "test.h" static const double areasKm2[] = { diff --git a/src/apps/testapps/testH3CellAreaExhaustive.c b/src/apps/testapps/testH3CellAreaExhaustive.c index 13261dae15..784aa51bde 100644 --- a/src/apps/testapps/testH3CellAreaExhaustive.c +++ b/src/apps/testapps/testH3CellAreaExhaustive.c @@ -23,6 +23,7 @@ #include +#include "adder.h" #include "constants.h" #include "h3api.h" #include "iterators.h" @@ -136,15 +137,16 @@ static void cell_area_assert(H3Index cell) { */ static void earth_area_test(int res, H3Error (*cell_area)(H3Index, double *), double target, double tol) { - double area = 0.0; - for (IterCellsResolution iter = iterInitRes(res); iter.h; - iterStepRes(&iter)) { + Adder adder = {}; + IterCellsResolution iter = iterInitRes(res); + + for (; iter.h; iterStepRes(&iter)) { double cellArea; t_assertSuccess((*cell_area)(iter.h, &cellArea)); - area += cellArea; + kadd(&adder, cellArea); } - t_assert(fabs(area - target) < tol, + t_assert(fabs(adder.sum - target) < tol, "sum of all cells should give earth area"); } @@ -176,28 +178,24 @@ SUITE(h3CellAreaExhaustive) { double km2 = rads2 * EARTH_RADIUS_KM * EARTH_RADIUS_KM; double m2 = km2 * 1000 * 1000; - // Notice the drop in accuracy at resolution 1. - // I think this has something to do with Class II vs Class III - // resolutions. - earth_area_test(0, H3_EXPORT(cellAreaRads2), rads2, 1e-14); earth_area_test(0, H3_EXPORT(cellAreaKm2), km2, 1e-6); earth_area_test(0, H3_EXPORT(cellAreaM2), m2, 1e0); - earth_area_test(1, H3_EXPORT(cellAreaRads2), rads2, 1e-9); - earth_area_test(1, H3_EXPORT(cellAreaKm2), km2, 1e-1); - earth_area_test(1, H3_EXPORT(cellAreaM2), m2, 1e5); + earth_area_test(1, H3_EXPORT(cellAreaRads2), rads2, 1e-14); + earth_area_test(1, H3_EXPORT(cellAreaKm2), km2, 1e-6); + earth_area_test(1, H3_EXPORT(cellAreaM2), m2, 1e0); - earth_area_test(2, H3_EXPORT(cellAreaRads2), rads2, 1e-12); - earth_area_test(2, H3_EXPORT(cellAreaKm2), km2, 1e-5); + earth_area_test(2, H3_EXPORT(cellAreaRads2), rads2, 1e-14); + earth_area_test(2, H3_EXPORT(cellAreaKm2), km2, 1e-6); earth_area_test(2, H3_EXPORT(cellAreaM2), m2, 1e0); - earth_area_test(3, H3_EXPORT(cellAreaRads2), rads2, 1e-11); - earth_area_test(3, H3_EXPORT(cellAreaKm2), km2, 1e-3); - earth_area_test(3, H3_EXPORT(cellAreaM2), m2, 1e3); + earth_area_test(3, H3_EXPORT(cellAreaRads2), rads2, 1e-14); + earth_area_test(3, H3_EXPORT(cellAreaKm2), km2, 1e-6); + earth_area_test(3, H3_EXPORT(cellAreaM2), m2, 1e0); - earth_area_test(4, H3_EXPORT(cellAreaRads2), rads2, 1e-11); - earth_area_test(4, H3_EXPORT(cellAreaKm2), km2, 1e-3); - earth_area_test(4, H3_EXPORT(cellAreaM2), m2, 1e2); + earth_area_test(4, H3_EXPORT(cellAreaRads2), rads2, 1e-14); + earth_area_test(4, H3_EXPORT(cellAreaKm2), km2, 1e-6); + earth_area_test(4, H3_EXPORT(cellAreaM2), m2, 1e0); } } diff --git a/src/h3lib/include/adder.h b/src/h3lib/include/adder.h new file mode 100644 index 0000000000..715470dd67 --- /dev/null +++ b/src/h3lib/include/adder.h @@ -0,0 +1,67 @@ +/* + * Copyright 2025 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#ifndef ADDER_H +#define ADDER_H + +/* +Header-only implementation of a "compensated summation" algorithm (Kahan +summation), which allows us to add up sequences of floating-point numbers +with better precision than naive summation, especially when the terms in +the sum vary significantly in magnitude. +See: https://en.wikipedia.org/wiki/Kahan_summation_algorithm + +This is useful when computing the area of (multi)polygons, which +often involves adding many small terms to a large aggregate. For example, +D3 uses an improved accuracy summation when computing polygonal areas via its +`Adder` class: +https://github.com/d3/d3-geo/blob/main/src/area.js + +There are a few potential algorithms we might consider for summation: + +1. Naive sum +2. Kahan summation +3. Neumaier summation +4. Other approaches like pairwise summation, or Python's `fsum` + +We considered the first three for simplicity, and settled on Kahan summation: +it achieves noticeably better accuracy than naive summation and almost as +good accuracy as Neumaier, while being only slightly slower than naive and +slightly faster and simpler than Neumaier. + +See also: https://github.com/python/cpython/issues/100425 for discussion of +tradeoffs between Kahan, Neumaier, and `fsum`. + +## Usage + +Adder adder = {}; // initialize +kadd(&adder, x); // add x to summation +out = adder.sum; // extract final sum +*/ + +typedef struct { + double sum; // running total. public + double _c; // compensation term. private +} Adder; + +static inline void kadd(Adder *adder, double x) { + double y = x - adder->_c; + double t = adder->sum + y; + adder->_c = (t - adder->sum) - y; + adder->sum = t; +} + +#endif diff --git a/src/h3lib/include/area.h b/src/h3lib/include/area.h new file mode 100644 index 0000000000..c561f53773 --- /dev/null +++ b/src/h3lib/include/area.h @@ -0,0 +1,27 @@ +/* + * Copyright 2025 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +/** @file area.h + * @brief Area computation functions + */ + +#ifndef AREA_H +#define AREA_H + +#include "h3api.h" + +H3Error geoLoopAreaRads2(GeoLoop loop, double *out); + +#endif diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c new file mode 100644 index 0000000000..166059fcd9 --- /dev/null +++ b/src/h3lib/lib/area.c @@ -0,0 +1,156 @@ +/* + * Copyright 2025 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +/** @file area.c + * @brief Algorithms for computing areas of regions on a sphere (GeoLoop, + * cells, polygons, multipolygons, etc.) + */ + +#include "area.h" + +#include + +#include "adder.h" +#include "constants.h" +#include "h3Assert.h" +#include "h3api.h" + +// Cagnoli contribution for edge arc x to y, following d3-geo’s +// area implementation: +// https://github.com/d3/d3-geo/blob/8c53a90ae70c94bace73ecb02f2c792c649c86ba/src/area.js#L51-L70 +static inline double cagnoli(LatLng x, LatLng y) { + x.lat = x.lat / 2.0 + M_PI / 4.0; + y.lat = y.lat / 2.0 + M_PI / 4.0; + + double sa = sin(x.lat) * sin(y.lat); + double ca = cos(x.lat) * cos(y.lat); + + double d = y.lng - x.lng; + double sd = sin(d); + double cd = cos(d); + + return -2.0 * atan2(sa * sd, sa * cd + ca); +} + +/** + * Area in radians^2 enclosed by vertices in GeoLoop. + * + * The GeoLoop should represent a simple curve with no self-intersections. + * Vertices should be ordered according to the "right hand rule". + * That is, if you are looking from outer space at a spherical polygon loop + * on the surface of the earth whose interior is contained within a hemisphere, + * then the vertices should be ordered counter-clockwise. The interior of the + * loop is to the left of a person walking along the boundary of the polygon + * in the counter-clockwise direction. + * + * Note that GeoLoops do not need to repeat the first vertex at the end of the + * array to close the loop; this is done automatically. + * + * The edge arcs between adjacent vertices are assumed to be the shortest + * geodesic path between them; that is, all arcs are interpreted to be less + * than 180 degrees or pi radians. + * Avoid arcs that are exactly pi (i.e, two antipodal vertices). + * "Large" polygon loops (e.g., that cannot be contained in a hemisphere) can + * still be constructed by using intermediate vertices with arcs less than + * 180 degrees, and the loop area will still be computed correctly. + * + * The area of the entire globe is 4*pi radians^2. If, for example, you have a + * small GeoLoop with area `a << 4*pi` and then reverse the order of the + * vertices, you produce GeoLoop with area `4*pi - a`, since, by the right hand + * rule, the new loop's interior is the majority of the globe, + * or "everything except the original polygon". + * Note that the area enclosed by the loop is determined by the vertex order; + * this function does **not** return `min(a, 4*pi - a)`. + * + * @param loop GeoLoop of boundary vertices in counter-clockwise order + * @param out loop area in radians^2, in interval [0, 4*pi] + * @return E_SUCCESS on success, or an error code otherwise + */ +H3Error geoLoopAreaRads2(GeoLoop loop, double *out) { + // Use `Adder` to improve numerical accuracy of the sum of many Cagnoli + // terms + Adder adder = {}; + + for (int i = 0; i < loop.numVerts; i++) { + int j = (i + 1) % loop.numVerts; + kadd(&adder, cagnoli(loop.verts[i], loop.verts[j])); + } + + // The Cagnoli sum above yields a signed area, with the sign switching + // with the orientation of the vertices. Since we want our area to always be + // positive, we normalize into [0, 4*pi] by adding 4*pi when the signed + // area is negative. + if (adder.sum < 0.0) { + kadd(&adder, 4.0 * M_PI); + } + + *out = adder.sum; + return E_SUCCESS; +} + +/** + * Area of H3 cell in radians^2. + * + * Uses `geoLoopAreaRads2` to compute cell area. + * + * @param cell H3 cell + * @param out cell area in radians^2 + * @return E_SUCCESS on success, or an error code otherwise + */ +H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) { + CellBoundary cb; + H3Error err = H3_EXPORT(cellToBoundary)(cell, &cb); + if (err) { + return err; + } + + GeoLoop loop = {.verts = cb.verts, .numVerts = cb.numVerts}; + err = geoLoopAreaRads2(loop, out); + if (NEVER(err)) { + return err; + } + + return E_SUCCESS; +} + +/** + * Area of H3 cell in kilometers^2. + * + * @param cell H3 cell + * @param out cell area in kilometers^2 + * @return E_SUCCESS on success, or an error code otherwise + */ +H3Error H3_EXPORT(cellAreaKm2)(H3Index cell, double *out) { + H3Error err = H3_EXPORT(cellAreaRads2)(cell, out); + if (!err) { + *out *= EARTH_RADIUS_KM * EARTH_RADIUS_KM; + } + return err; +} + +/** + * Area of H3 cell in meters^2. + * + * @param cell H3 cell + * @param out cell area in meters^2 + * @return E_SUCCESS on success, or an error code otherwise + */ +H3Error H3_EXPORT(cellAreaM2)(H3Index cell, double *out) { + H3Error err = H3_EXPORT(cellAreaKm2)(cell, out); + if (!err) { + *out *= 1000 * 1000; + } + return err; +} diff --git a/src/h3lib/lib/latLng.c b/src/h3lib/lib/latLng.c index f90d9ee2b4..082db77460 100644 --- a/src/h3lib/lib/latLng.c +++ b/src/h3lib/lib/latLng.c @@ -353,110 +353,6 @@ H3Error H3_EXPORT(getNumCells)(int res, int64_t *out) { return E_SUCCESS; } -/** - * Surface area in radians^2 of spherical triangle on unit sphere. - * - * For the math, see: - * https://en.wikipedia.org/wiki/Spherical_trigonometry#Area_and_spherical_excess - * - * @param a length of triangle side A in radians - * @param b length of triangle side B in radians - * @param c length of triangle side C in radians - * - * @return area in radians^2 of triangle on unit sphere - */ -double triangleEdgeLengthsToArea(double a, double b, double c) { - double s = (a + b + c) * 0.5; - - a = (s - a) * 0.5; - b = (s - b) * 0.5; - c = (s - c) * 0.5; - s = s * 0.5; - - return 4 * atan(sqrt(tan(s) * tan(a) * tan(b) * tan(c))); -} - -/** - * Compute area in radians^2 of a spherical triangle, given its vertices. - * - * @param a vertex lat/lng in radians - * @param b vertex lat/lng in radians - * @param c vertex lat/lng in radians - * - * @return area of triangle on unit sphere, in radians^2 - */ -double triangleArea(const LatLng *a, const LatLng *b, const LatLng *c) { - return triangleEdgeLengthsToArea(H3_EXPORT(greatCircleDistanceRads)(a, b), - H3_EXPORT(greatCircleDistanceRads)(b, c), - H3_EXPORT(greatCircleDistanceRads)(c, a)); -} - -/** - * Area of H3 cell in radians^2. - * - * The area is calculated by breaking the cell into spherical triangles and - * summing up their areas. Note that some H3 cells (hexagons and pentagons) - * are irregular, and have more than 6 or 5 sides. - * - * todo: optimize the computation by re-using the edges shared between triangles - * - * @param cell H3 cell - * @param out cell area in radians^2 - * @return E_SUCCESS on success, or an error code otherwise - */ -H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) { - LatLng c; - CellBoundary cb; - H3Error err = H3_EXPORT(cellToLatLng)(cell, &c); - if (err) { - return err; - } - err = H3_EXPORT(cellToBoundary)(cell, &cb); - if (NEVER(err)) { - // Uncoverable because cellToLatLng will have returned an error already - return err; - } - - double area = 0.0; - for (int i = 0; i < cb.numVerts; i++) { - int j = (i + 1) % cb.numVerts; - area += triangleArea(&cb.verts[i], &cb.verts[j], &c); - } - - *out = area; - return E_SUCCESS; -} - -/** - * Area of H3 cell in kilometers^2. - * - * @param cell H3 cell - * @param out cell area in kilometers^2 - * @return E_SUCCESS on success, or an error code otherwise - */ -H3Error H3_EXPORT(cellAreaKm2)(H3Index cell, double *out) { - H3Error err = H3_EXPORT(cellAreaRads2)(cell, out); - if (!err) { - *out = *out * EARTH_RADIUS_KM * EARTH_RADIUS_KM; - } - return err; -} - -/** - * Area of H3 cell in meters^2. - * - * @param cell H3 cell - * @param out cell area in meters^2 - * @return E_SUCCESS on success, or an error code otherwise - */ -H3Error H3_EXPORT(cellAreaM2)(H3Index cell, double *out) { - H3Error err = H3_EXPORT(cellAreaKm2)(cell, out); - if (!err) { - *out = *out * 1000 * 1000; - } - return err; -} - /** * Length of a directed edge in radians. *