Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
102 changes: 102 additions & 0 deletions x/h3go/area.go
Original file line number Diff line number Diff line change
@@ -0,0 +1,102 @@
/*
* Copyright 2026 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.
*/

package h3go

import "math"

// CellAreaRads2 returns the exact area of a cell in square radians.
func CellAreaRads2(c Cell) (float64, error) {
boundary, err := c.Boundary()
if err != nil {
return 0, err
}

return boundary.areaRads2(), nil
}

// CellAreaKm2 returns the exact area of a cell in square kilometers.
func CellAreaKm2(c Cell) (float64, error) {
rads2, err := CellAreaRads2(c)
if err != nil {
return 0, err
}

return rads2 * earthRadiusKm * earthRadiusKm, nil
}

// CellAreaM2 returns the exact area of a cell in square meters.
func CellAreaM2(c Cell) (float64, error) {
km2, err := CellAreaKm2(c)
if err != nil {
return 0, err
}

return km2 * 1000 * 1000, nil
}

// areaRads2 returns the area in square radians enclosed by the boundary loop. It
// sums the signed Cagnoli area contribution of each edge arc (assumed to be the
// shorter geodesic) with compensated summation, then normalizes a clockwise loop
// into [0, 4π] by adding the full-sphere area.
func (b CellBoundary) areaRads2() float64 {
var adder kahanAdder

verts := len(b)
for i := range verts {
next := (i + 1) % verts
adder.add(b[i].cagnoli(b[next]))
}

if adder.sum < 0 {
adder.add(2 * m2PI) // 4π, the area of the whole sphere
}

return adder.sum
}

// cagnoli returns the signed area contribution, in radians, of the boundary edge
// arc from ll to other (lat/lng in degrees), following the d3-geo spherical-area
// formulation.
func (ll LatLng) cagnoli(other LatLng) float64 {
lat := ll.Lat*degsToRads/2 + math.Pi/4
otherLat := other.Lat*degsToRads/2 + math.Pi/4

sa := math.Sin(lat) * math.Sin(otherLat)
ca := math.Cos(lat) * math.Cos(otherLat)

delta := (other.Lng - ll.Lng) * degsToRads
sinDelta := math.Sin(delta)
cosDelta := math.Cos(delta)

return -2 * math.Atan2(sa*sinDelta, sa*cosDelta+ca)
}

// kahanAdder accumulates a sum of floating-point terms using Kahan compensated
// summation, which preserves precision when adding many terms of varying
// magnitude (as a polygon area does).
type kahanAdder struct {
sum float64
comp float64 // running compensation for lost low-order bits
}

// add folds x into the running sum, carrying the rounding error forward.
func (a *kahanAdder) add(x float64) {
y := x - a.comp
t := a.sum + y
a.comp = (t - a.sum) - y
a.sum = t
}
156 changes: 156 additions & 0 deletions x/h3go/area_test.go
Original file line number Diff line number Diff line change
@@ -0,0 +1,156 @@
/*
* Copyright 2026 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.
*/

package h3go

import (
"math"
"testing"
)

// TestCellAreaKnownValues checks CellAreaRads2/Km2/M2 against reference values
// for hexagons and pentagons across several resolutions, including a base cell
// and cells whose boundaries span icosahedron faces.
func TestCellAreaKnownValues(t *testing.T) {
t.Parallel()

tests := map[string]struct {
giveCell string
wantRads2 float64
wantKm2 float64
wantM2 float64
}{
"hex_res9": {"8928308280fffff", 2.6952182709906241e-09, 0.1093981886467751, 109398.1886467751},
"base_res0": {"8001fffffffffff", 0.10116268528089556, 4106166.3344639186, 4106166334463.9185},
"hex_res5": {"85283473fffffff", 6.5310250106417195e-06, 265.09255812828178, 265092558.1282818},
"pent_res5": {"851c0003fffffff", 3.1482243104266963e-06, 127.78558260805931, 127785582.60805932},
"pent_res0": {"8009fffffffffff", 0.06312389871006796, 2562182.1629554993, 2562182162955.499},
}

for name, tt := range tests {
t.Run(name, func(t *testing.T) {
t.Parallel()

c := CellFromString(tt.giveCell)

rads2, err := CellAreaRads2(c)
if err != nil {
t.Fatalf("CellAreaRads2: %v", err)
}

assertRelClose(t, rads2, tt.wantRads2, "rads2")

km2, err := CellAreaKm2(c)
if err != nil {
t.Fatalf("CellAreaKm2: %v", err)
}

assertRelClose(t, km2, tt.wantKm2, "km2")

m2, err := CellAreaM2(c)
if err != nil {
t.Fatalf("CellAreaM2: %v", err)
}

assertRelClose(t, m2, tt.wantM2, "m2")
})
}
}

// TestCellAreaRes0SumsToSphere checks that the areas of all resolution-0 cells
// sum to the area of the unit sphere (4π), a global consistency property.
func TestCellAreaRes0SumsToSphere(t *testing.T) {
t.Parallel()

res0, err := Res0Cells()
if err != nil {
t.Fatalf("Res0Cells: %v", err)
}

var sum float64

for _, c := range res0 {
area, err := CellAreaRads2(c)
if err != nil {
t.Fatalf("CellAreaRads2(%015x): %v", uint64(c), err)
}

if area <= 0 {
t.Fatalf("CellAreaRads2(%015x) = %v, want positive", uint64(c), area)
}

sum += area
}

assertRelClose(t, sum, 4*math.Pi, "res0 area sum")
}

// TestCellAreaInvalidBaseCell covers the error branch of CellAreaRads2/Km2/M2,
// reached when the boundary cannot be computed for an out-of-range base cell.
func TestCellAreaInvalidBaseCell(t *testing.T) {
t.Parallel()

bad := Cell(h3Init) | Cell(cellMode)<<modeOffset | Cell(numBaseCells)<<baseCellOffset

if _, err := CellAreaRads2(bad); err == nil {
t.Fatal("CellAreaRads2: got nil error, want failure")
}

if _, err := CellAreaKm2(bad); err == nil {
t.Fatal("CellAreaKm2: got nil error, want failure")
}

if _, err := CellAreaM2(bad); err == nil {
t.Fatal("CellAreaM2: got nil error, want failure")
}
}

// TestBoundaryAreaClockwiseNormalizes covers the clockwise-loop branch of
// areaRads2: a loop wound clockwise has negative signed area and must be
// normalized into [0, 4π] by adding the full-sphere area. (Real cell boundaries
// are wound counterclockwise, so only a reversed loop reaches this branch.)
func TestBoundaryAreaClockwiseNormalizes(t *testing.T) {
t.Parallel()

ccw := CellBoundary{
{Lat: 0, Lng: 0},
{Lat: 0, Lng: 1},
{Lat: 1, Lng: 1},
{Lat: 1, Lng: 0},
}

cw := CellBoundary{ccw[3], ccw[2], ccw[1], ccw[0]}

small := ccw.areaRads2()
large := cw.areaRads2()

if small <= 0 || small >= 2*math.Pi {
t.Fatalf("ccw area = %v, want a small positive value", small)
}

assertRelClose(t, small+large, 4*math.Pi, "ccw + cw area")
}

// assertRelClose fails if got and want differ by more than a small relative
// tolerance, suitable for comparing values that span many orders of magnitude.
func assertRelClose(t *testing.T, got, want float64, label string) {
t.Helper()

const tol = 1e-9
if math.Abs(got-want) > tol*math.Max(1, math.Abs(want)) {
t.Fatalf("%s = %.17g, want %.17g", label, got, want)
}
}
4 changes: 4 additions & 0 deletions x/h3go/h3go.go
Original file line number Diff line number Diff line change
Expand Up @@ -118,6 +118,10 @@ const (
// radsToDegs converts radians to degrees by multiplying radians by this constant.
radsToDegs = 180.0 / math.Pi

// earthRadiusKm is the authalic (equal-area) radius of the Earth in
// kilometers, used to convert spherical measures to physical units.
earthRadiusKm = 6371.007180918475

// h3Init has all 15 digit slots set to 7 (invalid); mode/res/base cell are 0.
h3Init = 35184372088831

Expand Down
112 changes: 112 additions & 0 deletions x/h3go/latlng.go
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,118 @@ func (c Cell) LatLng() (LatLng, error) {
return fijk.toVec3(c.Resolution()).toLatLng(), nil
}

// --- Distances ---

// GreatCircleDistanceRads returns the great-circle distance between two points
// in radians, using the haversine formula. The points are given in degrees.
func GreatCircleDistanceRads(a, b LatLng) float64 {
aLat := a.Lat * degsToRads
aLng := a.Lng * degsToRads
bLat := b.Lat * degsToRads
bLng := b.Lng * degsToRads

sinLat := math.Sin((bLat - aLat) / 2)
sinLng := math.Sin((bLng - aLng) / 2)

h := sinLat*sinLat + math.Cos(aLat)*math.Cos(bLat)*sinLng*sinLng

return 2 * math.Atan2(math.Sqrt(h), math.Sqrt(1-h))
}

// GreatCircleDistanceKm returns the great-circle distance between two points in
// kilometers. The points are given in degrees.
func GreatCircleDistanceKm(a, b LatLng) float64 {
return GreatCircleDistanceRads(a, b) * earthRadiusKm
}

// GreatCircleDistanceM returns the great-circle distance between two points in
// meters. The points are given in degrees.
func GreatCircleDistanceM(a, b LatLng) float64 {
return GreatCircleDistanceKm(a, b) * 1000
}

// --- Average cell sizes ---

// hexAreaAvgKm2 holds the average hexagon area at each resolution in square
// kilometers, indexed by resolution.
var hexAreaAvgKm2 = [maxResolution + 1]float64{
4.357449416078383e+06, 6.097884417941332e+05, 8.680178039899720e+04,
1.239343465508816e+04, 1.770347654491307e+03, 2.529038581819449e+02,
3.612906216441245e+01, 5.161293359717191e+00, 7.373275975944177e-01,
1.053325134272067e-01, 1.504750190766435e-02, 2.149643129451879e-03,
3.070918756316060e-04, 4.387026794728296e-05, 6.267181135324313e-06,
8.953115907605790e-07,
}

// hexAreaAvgM2 holds the average hexagon area at each resolution in square
// meters, indexed by resolution.
var hexAreaAvgM2 = [maxResolution + 1]float64{
4.357449416078390e+12, 6.097884417941339e+11, 8.680178039899731e+10,
1.239343465508818e+10, 1.770347654491309e+09, 2.529038581819452e+08,
3.612906216441250e+07, 5.161293359717198e+06, 7.373275975944188e+05,
1.053325134272069e+05, 1.504750190766437e+04, 2.149643129451882e+03,
3.070918756316063e+02, 4.387026794728301e+01, 6.267181135324322e+00,
8.953115907605802e-01,
}

// hexEdgeLenAvgKm holds the average hexagon edge length at each resolution in
// kilometers, indexed by resolution.
var hexEdgeLenAvgKm = [maxResolution + 1]float64{
1281.256011, 483.0568391, 182.5129565, 68.97922179,
26.07175968, 9.854090990, 3.724532667, 1.406475763,
0.531414010, 0.200786148, 0.075863783, 0.028663897,
0.010830188, 0.004092010, 0.001546100, 0.000584169,
}

// hexEdgeLenAvgM holds the average hexagon edge length at each resolution in
// meters, indexed by resolution.
var hexEdgeLenAvgM = [maxResolution + 1]float64{
1281256.011, 483056.8391, 182512.9565, 68979.22179,
26071.75968, 9854.090990, 3724.532667, 1406.475763,
531.4140101, 200.7861476, 75.86378287, 28.66389748,
10.83018784, 4.092010473, 1.546099657, 0.584168630,
}

// HexagonAreaAvgKm2 returns the average area of a hexagon at the given
// resolution in square kilometers.
func HexagonAreaAvgKm2(res int) (float64, error) {
if res < 0 || res > maxResolution {
return 0, ErrResolutionDomain
}

return hexAreaAvgKm2[res], nil
}

// HexagonAreaAvgM2 returns the average area of a hexagon at the given resolution
// in square meters.
func HexagonAreaAvgM2(res int) (float64, error) {
if res < 0 || res > maxResolution {
return 0, ErrResolutionDomain
}

return hexAreaAvgM2[res], nil
}

// HexagonEdgeLengthAvgKm returns the average edge length of a hexagon at the
// given resolution in kilometers.
func HexagonEdgeLengthAvgKm(res int) (float64, error) {
if res < 0 || res > maxResolution {
return 0, ErrResolutionDomain
}

return hexEdgeLenAvgKm[res], nil
}

// HexagonEdgeLengthAvgM returns the average edge length of a hexagon at the given
// resolution in meters.
func HexagonEdgeLengthAvgM(res int) (float64, error) {
if res < 0 || res > maxResolution {
return 0, ErrResolutionDomain
}

return hexEdgeLenAvgM[res], nil
}

// --- Vec3d math ---

func latLngToVec3(lat, lng float64) vec3d {
Expand Down
Loading
Loading