diff --git a/x/h3go/area.go b/x/h3go/area.go new file mode 100644 index 0000000..faaf50f --- /dev/null +++ b/x/h3go/area.go @@ -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 +} diff --git a/x/h3go/area_test.go b/x/h3go/area_test.go new file mode 100644 index 0000000..096d364 --- /dev/null +++ b/x/h3go/area_test.go @@ -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)<= 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) + } +} diff --git a/x/h3go/h3go.go b/x/h3go/h3go.go index 8aafb62..585cb11 100644 --- a/x/h3go/h3go.go +++ b/x/h3go/h3go.go @@ -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 diff --git a/x/h3go/latlng.go b/x/h3go/latlng.go index d24fbf0..eed0dd4 100644 --- a/x/h3go/latlng.go +++ b/x/h3go/latlng.go @@ -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 { diff --git a/x/h3go/latlng_test.go b/x/h3go/latlng_test.go index c70f763..d4f6b77 100644 --- a/x/h3go/latlng_test.go +++ b/x/h3go/latlng_test.go @@ -219,3 +219,141 @@ func TestCellToLatLngInvalidIndex(t *testing.T) { t.Fatalf("CellToLatLng(0x7fffffffffffffff): got %v, want %v", err, ErrCellInvalid) } } + +// TestGreatCircleDistance checks the haversine distance in radians, kilometers, +// and meters against a known city pair, plus the zero-distance and symmetry +// properties. +func TestGreatCircleDistance(t *testing.T) { + t.Parallel() + + sf := LatLng{Lat: 37.7749, Lng: -122.4194} + ny := LatLng{Lat: 40.7128, Lng: -74.0060} + + t.Run("known_pair", func(t *testing.T) { + t.Parallel() + + assertRelClose(t, GreatCircleDistanceRads(sf, ny), 0.64810644562192898, "rads") + assertRelClose(t, GreatCircleDistanceKm(sf, ny), 4129.090819109696, "km") + assertRelClose(t, GreatCircleDistanceM(sf, ny), 4129090.819109696, "m") + }) + + t.Run("zero_distance", func(t *testing.T) { + t.Parallel() + + if got := GreatCircleDistanceRads(sf, sf); math.Abs(got) > 1e-12 { + t.Fatalf("GreatCircleDistanceRads(sf, sf) = %v, want ~0", got) + } + }) + + t.Run("symmetric", func(t *testing.T) { + t.Parallel() + + if a, b := GreatCircleDistanceRads(sf, ny), GreatCircleDistanceRads(ny, sf); a != b { + t.Fatalf("distance not symmetric: %v vs %v", a, b) + } + }) +} + +// TestHexagonAreaAvg checks the average-area getters: a known resolution-0 value, +// strictly decreasing area with finer resolution, and the resolution-domain error. +func TestHexagonAreaAvg(t *testing.T) { + t.Parallel() + + t.Run("known_and_monotonic", func(t *testing.T) { + t.Parallel() + + km0, err := HexagonAreaAvgKm2(0) + if err != nil { + t.Fatalf("HexagonAreaAvgKm2(0): %v", err) + } + + assertRelClose(t, km0, 4.357449416078383e+06, "area km2 res0") + + prevKm, prevM := math.Inf(1), math.Inf(1) + + for res := 0; res <= maxResolution; res++ { + km, err := HexagonAreaAvgKm2(res) + if err != nil { + t.Fatalf("HexagonAreaAvgKm2(%d): %v", res, err) + } + + m2, err := HexagonAreaAvgM2(res) + if err != nil { + t.Fatalf("HexagonAreaAvgM2(%d): %v", res, err) + } + + if km >= prevKm || m2 >= prevM { + t.Fatalf("res %d: area not decreasing (km=%v prevKm=%v)", res, km, prevKm) + } + + assertRelClose(t, m2, km*1e6, "m2 vs km2*1e6") + prevKm, prevM = km, m2 + } + }) + + t.Run("out_of_range", func(t *testing.T) { + t.Parallel() + + for _, res := range []int{-1, maxResolution + 1} { + if _, err := HexagonAreaAvgKm2(res); !errors.Is(err, ErrResolutionDomain) { + t.Fatalf("HexagonAreaAvgKm2(%d): got %v, want %v", res, err, ErrResolutionDomain) + } + + if _, err := HexagonAreaAvgM2(res); !errors.Is(err, ErrResolutionDomain) { + t.Fatalf("HexagonAreaAvgM2(%d): got %v, want %v", res, err, ErrResolutionDomain) + } + } + }) +} + +// TestHexagonEdgeLengthAvg checks the average-edge-length getters: a known +// resolution-0 value, strictly decreasing length with finer resolution, and the +// resolution-domain error. +func TestHexagonEdgeLengthAvg(t *testing.T) { + t.Parallel() + + t.Run("known_and_monotonic", func(t *testing.T) { + t.Parallel() + + km0, err := HexagonEdgeLengthAvgKm(0) + if err != nil { + t.Fatalf("HexagonEdgeLengthAvgKm(0): %v", err) + } + + assertRelClose(t, km0, 1281.256011, "edge km res0") + + prevKm, prevM := math.Inf(1), math.Inf(1) + + for res := 0; res <= maxResolution; res++ { + km, err := HexagonEdgeLengthAvgKm(res) + if err != nil { + t.Fatalf("HexagonEdgeLengthAvgKm(%d): %v", res, err) + } + + m, err := HexagonEdgeLengthAvgM(res) + if err != nil { + t.Fatalf("HexagonEdgeLengthAvgM(%d): %v", res, err) + } + + if km >= prevKm || m >= prevM { + t.Fatalf("res %d: length not decreasing (km=%v prevKm=%v)", res, km, prevKm) + } + + prevKm, prevM = km, m + } + }) + + t.Run("out_of_range", func(t *testing.T) { + t.Parallel() + + for _, res := range []int{-1, maxResolution + 1} { + if _, err := HexagonEdgeLengthAvgKm(res); !errors.Is(err, ErrResolutionDomain) { + t.Fatalf("HexagonEdgeLengthAvgKm(%d): got %v, want %v", res, err, ErrResolutionDomain) + } + + if _, err := HexagonEdgeLengthAvgM(res); !errors.Is(err, ErrResolutionDomain) { + t.Fatalf("HexagonEdgeLengthAvgM(%d): got %v, want %v", res, err, ErrResolutionDomain) + } + } + }) +} diff --git a/x/h3go/paritytest/measures_test.go b/x/h3go/paritytest/measures_test.go new file mode 100644 index 0000000..ee5579a --- /dev/null +++ b/x/h3go/paritytest/measures_test.go @@ -0,0 +1,133 @@ +/* + * 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 paritytest + +import ( + "math" + "testing" + + "github.com/uber/h3-go/v4" + "github.com/uber/h3-go/v4/x/h3go" +) + +// measureRelTolerance is the maximum allowed relative difference between the +// pure-Go and cgo measure results. The two implementations run the same formulas, +// but Go's math package and C's libm differ by a few ULPs in the transcendental +// functions; that fixed-magnitude error grows in relative terms as cells shrink, +// reaching ~1e-6 at the finest resolution — so we compare to a small relative +// tolerance rather than bit-for-bit. measureAbsFloor handles the want==0 case +// (the distance from a point to itself). +const ( + measureRelTolerance = 1e-5 + measureAbsFloor = 1e-9 +) + +// TestCellAreaMatchesCgo asserts the pure-Go cell area matches the cgo-backed +// reference in radians², km², and m² for every cell in the shared corpus. +func TestCellAreaMatchesCgo(t *testing.T) { + t.Parallel() + + for _, ref := range referenceCorpus(t) { + goCell := h3goCell(ref) + + wantRads2, wErr := h3.CellAreaRads2(ref) + gotRads2, gErr := h3go.CellAreaRads2(goCell) + + if !bothErr(wErr, gErr) { + t.Fatalf("CellAreaRads2(%015x) error mismatch: cgo=%v h3go=%v", uint64(ref), wErr, gErr) + } + + assertRelClose(t, gotRads2, wantRads2, ref, "rads2") + + wantKm2, _ := h3.CellAreaKm2(ref) + gotKm2, _ := h3go.CellAreaKm2(goCell) + assertRelClose(t, gotKm2, wantKm2, ref, "km2") + + wantM2, _ := h3.CellAreaM2(ref) + gotM2, _ := h3go.CellAreaM2(goCell) + assertRelClose(t, gotM2, wantM2, ref, "m2") + } +} + +// TestGreatCircleDistanceMatchesCgo asserts the pure-Go haversine distance +// matches the cgo-backed reference (radians/km/m) for every pair of corpus +// points. +func TestGreatCircleDistanceMatchesCgo(t *testing.T) { + t.Parallel() + + for _, a := range corpusPoints { + for _, b := range corpusPoints { + goA := h3go.LatLng{Lat: a.Lat, Lng: a.Lng} + goB := h3go.LatLng{Lat: b.Lat, Lng: b.Lng} + + assertRelClose(t, h3go.GreatCircleDistanceRads(goA, goB), h3.GreatCircleDistanceRads(a, b), 0, "dist rads") + assertRelClose(t, h3go.GreatCircleDistanceKm(goA, goB), h3.GreatCircleDistanceKm(a, b), 0, "dist km") + assertRelClose(t, h3go.GreatCircleDistanceM(goA, goB), h3.GreatCircleDistanceM(a, b), 0, "dist m") + } + } +} + +// TestHexagonAvgMatchesCgo asserts the average area and edge-length getters match +// the cgo-backed reference for every resolution, including the out-of-range +// error. +func TestHexagonAvgMatchesCgo(t *testing.T) { + t.Parallel() + + for res := -1; res <= h3.MaxResolution+1; res++ { + areaKmWant, areaKmErr := h3.HexagonAreaAvgKm2(res) + areaKmGot, areaKmGotErr := h3go.HexagonAreaAvgKm2(res) + + if !bothErr(areaKmErr, areaKmGotErr) { + t.Fatalf("HexagonAreaAvgKm2(%d) error mismatch: cgo=%v h3go=%v", res, areaKmErr, areaKmGotErr) + } + + assertRelClose(t, areaKmGot, areaKmWant, 0, "areaAvgKm2") + + areaMWant, _ := h3.HexagonAreaAvgM2(res) + areaMGot, _ := h3go.HexagonAreaAvgM2(res) + assertRelClose(t, areaMGot, areaMWant, 0, "areaAvgM2") + + lenKmWant, _ := h3.HexagonEdgeLengthAvgKm(res) + lenKmGot, _ := h3go.HexagonEdgeLengthAvgKm(res) + assertRelClose(t, lenKmGot, lenKmWant, 0, "edgeLenAvgKm") + + lenMWant, _ := h3.HexagonEdgeLengthAvgM(res) + lenMGot, _ := h3go.HexagonEdgeLengthAvgM(res) + assertRelClose(t, lenMGot, lenMWant, 0, "edgeLenAvgM") + } +} + +// assertRelClose fails if got and want differ by more than measureRelTolerance +// relative to want (or by more than measureAbsFloor when want is zero). cell is +// included in the failure message when non-zero. +func assertRelClose(t *testing.T, got, want float64, cell h3.Cell, label string) { + t.Helper() + + diff := math.Abs(got - want) + + if want == 0 { + if diff > measureAbsFloor { + t.Fatalf("%s (cell %015x): got %.17g, want %.17g", label, uint64(cell), got, want) + } + + return + } + + if diff/math.Abs(want) > measureRelTolerance { + t.Fatalf("%s (cell %015x): got %.17g, want %.17g", label, uint64(cell), got, want) + } +}