/src/h3/src/h3lib/lib/area.c

Source
/*
 * 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 <math.h>

#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 GeoPolygon in radians^2.
 *
 * Outer GeoLoop vertices should be in counter-clockwise order.
 * Hole GeoLoop vertices should be in clockwise order.
 * Returned area is the area contained by the outer loop, minus
 * the areas of the holes. See `geoLoopAreaRads2` for the expected
 * form of GeoLoops.
 *
 * No check is made to ensure holes are disjoint or are contained
 * within the outer GeoLoop.
 *
 * @param   poly  GeoPolygon
 * @param    out  GeoPolygon area in radians^2
 * @return  E_SUCCESS on success; error code otherwise
 */
H3Error geoPolygonAreaRads2(GeoPolygon poly, double *out) {
    H3Error err;
    Adder adder = {};
    double term;

    err = geoLoopAreaRads2(poly.geoloop, &term);
    if (err) return err;
    kadd(&adder, term);

    for (int i = 0; i < poly.numHoles; i++) {
        err = geoLoopAreaRads2(poly.holes[i], &term);
        if (err) return err;

        // Due to clockwise order, holes will contribute area
        // of "everything except the hole", so adjust with -4*pi term.
        kadd(&adder, term);
        kadd(&adder, -4.0 * M_PI);
    }

    *out = adder.sum;

    return E_SUCCESS;
}

/**
 * Area of GeoMultiPolygon in radians^2.
 *
 * Area is the sum of the areas of the polygons contained by `mpoly`.
 * See `geoPolygonAreaRads2` for expected polygon format.
 * No check is made to ensure polygons are disjoint.
 *
 * @param   mpoly  GeoMultiPolygon
 * @param     out  GeoMultiPolygon area in radians^2
 * @return  E_SUCCESS on success; error code otherwise
 */
H3Error geoMultiPolygonAreaRads2(GeoMultiPolygon mpoly, double *out) {
    H3Error err;
    Adder adder = {};
    double term;

    for (int i = 0; i < mpoly.numPolygons; i++) {
        err = geoPolygonAreaRads2(mpoly.polygons[i], &term);
        if (err) return err;
        kadd(&adder, term);
    }

    *out = adder.sum;

    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;
}

Coverage Report

Created: 2026-01-09 06:51

Line	Count	Source
1		/*
2		* Copyright 2025 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 area.c
17		* @brief Algorithms for computing areas of regions on a sphere (GeoLoop,
18		* cells, polygons, multipolygons, etc.)
19		*/
20
21		#include "area.h"
22
23		#include <math.h>
24
25		#include "adder.h"
26		#include "constants.h"
27		#include "h3Assert.h"
28		#include "h3api.h"
29
30		// Cagnoli contribution for edge arc x to y, following d3-geo’s
31		// area implementation:
32		// https://github.com/d3/d3-geo/blob/8c53a90ae70c94bace73ecb02f2c792c649c86ba/src/area.js#L51-L70
33	4.96k	static inline double cagnoli(LatLng x, LatLng y) {
34	4.96k	x.lat = x.lat / 2.0 + M_PI / 4.0;
35	4.96k	y.lat = y.lat / 2.0 + M_PI / 4.0;
36
37	4.96k	double sa = sin(x.lat) * sin(y.lat);
38	4.96k	double ca = cos(x.lat) * cos(y.lat);
39
40	4.96k	double d = y.lng - x.lng;
41	4.96k	double sd = sin(d);
42	4.96k	double cd = cos(d);
43
44	4.96k	return -2.0 * atan2(sa * sd, sa * cd + ca);
45	4.96k	}
46
47		/**
48		* Area in radians^2 enclosed by vertices in GeoLoop.
49		*
50		* The GeoLoop should represent a simple curve with no self-intersections.
51		* Vertices should be ordered according to the "right hand rule".
52		* That is, if you are looking from outer space at a spherical polygon loop
53		* on the surface of the earth whose interior is contained within a hemisphere,
54		* then the vertices should be ordered counter-clockwise. The interior of the
55		* loop is to the left of a person walking along the boundary of the polygon
56		* in the counter-clockwise direction.
57		*
58		* Note that GeoLoops do not need to repeat the first vertex at the end of the
59		* array to close the loop; this is done automatically.
60		*
61		* The edge arcs between adjacent vertices are assumed to be the shortest
62		* geodesic path between them; that is, all arcs are interpreted to be less
63		* than 180 degrees or pi radians.
64		* Avoid arcs that are exactly pi (i.e, two antipodal vertices).
65		* "Large" polygon loops (e.g., that cannot be contained in a hemisphere) can
66		* still be constructed by using intermediate vertices with arcs less than
67		* 180 degrees, and the loop area will still be computed correctly.
68		*
69		* The area of the entire globe is 4*pi radians^2. If, for example, you have a
70		* small GeoLoop with area `a << 4*pi` and then reverse the order of the
71		* vertices, you produce GeoLoop with area `4*pi - a`, since, by the right hand
72		* rule, the new loop's interior is the majority of the globe,
73		* or "everything except the original polygon".
74		* Note that the area enclosed by the loop is determined by the vertex order;
75		* this function does not return `min(a, 4*pi - a)`.
76		*
77		* @param loop GeoLoop of boundary vertices in counter-clockwise order
78		* @param out loop area in radians^2, in interval [0, 4*pi]
79		* @return E_SUCCESS on success, or an error code otherwise
80		*/
81	801	H3Error geoLoopAreaRads2(GeoLoop loop, double *out) {
82		// Use `Adder` to improve numerical accuracy of the sum of many Cagnoli
83		// terms
84	801	Adder adder = {};
85
86	5.76k	for (int i = 0; i < loop.numVerts; i++) {
87	4.96k	int j = (i + 1) % loop.numVerts;
88	4.96k	kadd(&adder, cagnoli(loop.verts[i], loop.verts[j]));
89	4.96k	}
90
91		// The Cagnoli sum above yields a signed area, with the sign switching
92		// with the orientation of the vertices. Since we want our area to always be
93		// positive, we normalize into [0, 4pi] by adding 4pi when the signed
94		// area is negative.
95	801	if (adder.sum < 0.0) {
96	24	kadd(&adder, 4.0 * M_PI);
97	24	}
98
99	801	*out = adder.sum;
100	801	return E_SUCCESS;
101	801	}
102
103		/**
104		* Area of H3 cell in radians^2.
105		*
106		* Uses `geoLoopAreaRads2` to compute cell area.
107		*
108		* @param cell H3 cell
109		* @param out cell area in radians^2
110		* @return E_SUCCESS on success, or an error code otherwise
111		*/
112	807	H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) {
113	807	CellBoundary cb;
114	807	H3Error err = H3_EXPORT(cellToBoundary)(cell, &cb);
115	807	if (err) {
116	6	return err;
117	6	}
118
119	801	GeoLoop loop = {.verts = cb.verts, .numVerts = cb.numVerts};
120	801	err = geoLoopAreaRads2(loop, out);
121	801	if (NEVER(err)) {
122	0	return err;
123	0	}
124
125	801	return E_SUCCESS;
126	801	}
127
128		/**
129		* Area of GeoPolygon in radians^2.
130		*
131		* Outer GeoLoop vertices should be in counter-clockwise order.
132		* Hole GeoLoop vertices should be in clockwise order.
133		* Returned area is the area contained by the outer loop, minus
134		* the areas of the holes. See `geoLoopAreaRads2` for the expected
135		* form of GeoLoops.
136		*
137		* No check is made to ensure holes are disjoint or are contained
138		* within the outer GeoLoop.
139		*
140		* @param poly GeoPolygon
141		* @param out GeoPolygon area in radians^2
142		* @return E_SUCCESS on success; error code otherwise
143		*/
144	0	H3Error geoPolygonAreaRads2(GeoPolygon poly, double *out) {
145	0	H3Error err;
146	0	Adder adder = {};
147	0	double term;
148
149	0	err = geoLoopAreaRads2(poly.geoloop, &term);
150	0	if (err) return err;
151	0	kadd(&adder, term);
152
153	0	for (int i = 0; i < poly.numHoles; i++) {
154	0	err = geoLoopAreaRads2(poly.holes[i], &term);
155	0	if (err) return err;
156
157		// Due to clockwise order, holes will contribute area
158		// of "everything except the hole", so adjust with -4*pi term.
159	0	kadd(&adder, term);
160	0	kadd(&adder, -4.0 * M_PI);
161	0	}
162
163	0	*out = adder.sum;
164
165	0	return E_SUCCESS;
166	0	}
167
168		/**
169		* Area of GeoMultiPolygon in radians^2.
170		*
171		* Area is the sum of the areas of the polygons contained by `mpoly`.
172		* See `geoPolygonAreaRads2` for expected polygon format.
173		* No check is made to ensure polygons are disjoint.
174		*
175		* @param mpoly GeoMultiPolygon
176		* @param out GeoMultiPolygon area in radians^2
177		* @return E_SUCCESS on success; error code otherwise
178		*/
179	0	H3Error geoMultiPolygonAreaRads2(GeoMultiPolygon mpoly, double *out) {
180	0	H3Error err;
181	0	Adder adder = {};
182	0	double term;
183
184	0	for (int i = 0; i < mpoly.numPolygons; i++) {
185	0	err = geoPolygonAreaRads2(mpoly.polygons[i], &term);
186	0	if (err) return err;
187	0	kadd(&adder, term);
188	0	}
189
190	0	*out = adder.sum;
191
192	0	return E_SUCCESS;
193	0	}
194
195		/**
196		* Area of H3 cell in kilometers^2.
197		*
198		* @param cell H3 cell
199		* @param out cell area in kilometers^2
200		* @return E_SUCCESS on success, or an error code otherwise
201		*/
202	538	H3Error H3_EXPORT(cellAreaKm2)(H3Index cell, double *out) {
203	538	H3Error err = H3_EXPORT(cellAreaRads2)(cell, out);
204	538	if (!err) {
205	534	out = EARTH_RADIUS_KM * EARTH_RADIUS_KM;
206	534	}
207	538	return err;
208	538	}
209
210		/**
211		* Area of H3 cell in meters^2.
212		*
213		* @param cell H3 cell
214		* @param out cell area in meters^2
215		* @return E_SUCCESS on success, or an error code otherwise
216		*/
217	269	H3Error H3_EXPORT(cellAreaM2)(H3Index cell, double *out) {
218	269	H3Error err = H3_EXPORT(cellAreaKm2)(cell, out);
219	269	if (!err) {
220	267	out = 1000 * 1000;
221	267	}
222	269	return err;
223	269	}