h3: keep cell areas exact across the antimeridian

area.go

@@ -106,7 +106,18 @@ func (v cagnoliVertex) cagnoli(other cagnoliVertex) float64 {
 	sa := v.sinHalfLat * other.sinHalfLat
 	ca := v.cosHalfLat * other.cosHalfLat
 
-	delta := (other.lng - v.lng) * DegsToRads
+	// For an edge crossing the antimeridian the raw difference is close to
+	// ±360°, and the true (tiny) longitude step is lost to the rounding of a
+	// value that large. Shift both longitudes to the same side of the
+	// antimeridian first so the subtraction is between two nearby values; sin
+	// and cos are 2π-periodic, so the result is mathematically unchanged.
+	deltaDegs := other.lng - v.lng
+	if deltaDegs > piDeg || deltaDegs < -piDeg {
+		deltaDegs = (other.lng - math.Copysign(piDeg, other.lng)) -
+			(v.lng - math.Copysign(piDeg, v.lng))
+	}
+
+	delta := deltaDegs * DegsToRads
 	sinDelta := math.Sin(delta)
 	cosDelta := math.Cos(delta)
 

area_test.go

@@ -18,6 +18,7 @@ package h3
 
 import (
 	"math"
+	"strconv"
 	"testing"
 )
 
@@ -185,6 +186,123 @@ func TestBoundaryAreaEmpty(t *testing.T) {
 	}
 }
 
+// TestCellAreaAntimeridian checks that a cell straddling the antimeridian has
+// the same area as a neighbouring cell at the same resolution that does not.
+// Subtracting raw longitudes across ±180° loses the tiny true step to the
+// rounding of a value near 360°, which made straddling cells too small.
+func TestCellAreaAntimeridian(t *testing.T) {
+	t.Parallel()
+
+	tests := map[string]struct {
+		giveStraddling string
+		giveNeighbor   string
+	}{
+		"res13": {"8d056898db628bf", "8d056898db62d7f"},
+		"res14": {"8e32b605c51e067", "8e32b605c51e15f"},
+		"res15": {"8f7f809893b5b9d", "8f7f809893b5b9c"},
+	}
+
+	for name, tt := range tests {
+		t.Run(name, func(t *testing.T) {
+			t.Parallel()
+
+			straddling, err := CellAreaRads2(CellFromString(tt.giveStraddling))
+			if err != nil {
+				t.Fatalf("CellAreaRads2(%s): %v", tt.giveStraddling, err)
+			}
+
+			neighbor, err := CellAreaRads2(CellFromString(tt.giveNeighbor))
+			if err != nil {
+				t.Fatalf("CellAreaRads2(%s): %v", tt.giveNeighbor, err)
+			}
+
+			const tol = 1e-5
+			if rel := math.Abs(straddling-neighbor) / neighbor; rel > tol {
+				t.Fatalf("straddling area %.17g, neighbor %.17g: relative difference %.2e > %.0e",
+					straddling, neighbor, rel, tol)
+			}
+		})
+	}
+}
+
+// TestCellAreaAntimeridianRotation checks, at every resolution, that a cell on
+// the antimeridian keeps its area when its boundary is rotated 180° about the
+// polar axis, which moves it onto the prime meridian where no edge crosses
+// ±180°. A rotation cannot change an area, so no reference values are needed.
+func TestCellAreaAntimeridianRotation(t *testing.T) {
+	t.Parallel()
+
+	for res := 0; res <= MaxResolution; res++ {
+		t.Run(strconv.Itoa(res), func(t *testing.T) {
+			t.Parallel()
+
+			for _, lat := range []float64{-60, 0, 45} {
+				cell, err := LatLngToCell(NewLatLng(lat, 180), res)
+				if err != nil {
+					t.Fatalf("LatLngToCell: %v", err)
+				}
+
+				boundary, err := cell.Boundary()
+				if err != nil {
+					t.Fatalf("Boundary: %v", err)
+				}
+
+				rotated := make(CellBoundary, len(boundary))
+				for i := range boundary {
+					rotated[i] = LatLng{Lat: boundary[i].Lat, Lng: rotateHalfTurn(boundary[i].Lng)}
+				}
+
+				const tol = 1e-6
+
+				got, want := boundary.areaRads2(), rotated.areaRads2()
+				if rel := math.Abs(got-want) / want; rel > tol {
+					t.Fatalf("%s: area %.17g, rotated %.17g: relative difference %.2e > %.0e",
+						cell, got, want, rel, tol)
+				}
+			}
+		})
+	}
+}
+
+// TestBoundaryAreaAntimeridianTriangle checks that a small triangle across
+// ±180° has the same area as the same triangle across 0°. A triangle rather
+// than a rectangle: a rectangle's two crossing edges have identical longitude
+// differences, so their rounding errors cancel and hide the defect.
+func TestBoundaryAreaAntimeridianTriangle(t *testing.T) {
+	t.Parallel()
+
+	const halfWidth = 1e-6
+
+	prime := CellBoundary{
+		{Lat: 0, Lng: -halfWidth},
+		{Lat: 0, Lng: halfWidth},
+		{Lat: 2 * halfWidth, Lng: 0},
+	}
+
+	anti := make(CellBoundary, len(prime))
+	for i := range prime {
+		anti[i] = LatLng{Lat: prime[i].Lat, Lng: rotateHalfTurn(prime[i].Lng)}
+	}
+
+	const tol = 1e-9
+
+	got, want := anti.areaRads2(), prime.areaRads2()
+	if rel := math.Abs(got-want) / want; rel > tol {
+		t.Fatalf("antimeridian triangle area %.17g, prime meridian %.17g: relative difference %.2e > %.0e",
+			got, want, rel, tol)
+	}
+}
+
+// rotateHalfTurn returns the longitude 180° east of lng, wrapped to (-180, 180].
+func rotateHalfTurn(lng float64) float64 {
+	rotated := lng + piDeg
+	if rotated > piDeg {
+		rotated -= twoPiDeg
+	}
+
+	return rotated
+}
+
 // 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) {

bbox.go

@@ -18,20 +18,6 @@ package h3
 
 import "math"
 
-// Longitude span constants in degrees, the units of the public API. They stand
-// in for the H3 C library's radian M_PI / M_2PI when detecting and normalizing
-// antimeridian-crossing geometry.
-const (
-	piDeg     = 180.0
-	twoPiDeg  = 360.0
-	halfPiDeg = 90.0
-
-	// twoPiRad is a full longitude turn in radians, used to normalize
-	// antimeridian-crossing geometry where the math must run at the radian scale
-	// of the H3 C library.
-	twoPiRad = 2 * math.Pi
-)
-
 // longitudeNormalization selects how a longitude is shifted so two bounding
 // boxes, either of which may cross the antimeridian, can be compared in one
 // frame of reference.

conformance/README.md

@@ -124,6 +124,7 @@ floats/edges.jsonl
 floats/vertexes.jsonl
 floats/distances.jsonl
 floats/resolutions.jsonl
+floats/curated.jsonl
 digests/resolutions.jsonl
 digests/baseCells.jsonl
 digests/patterns.jsonl
@@ -676,10 +677,11 @@ runner stops at the one-million tier unless
 ## Curated files
 
 Seeded sampling gives coverage density but is blind to the particular inputs
-that have broken implementations. Six files hold such inputs by hand:
+that have broken implementations. Seven files hold such inputs by hand:
 `inspection/curated.jsonl`, `hierarchy/curated.jsonl`,
-`traversal/curated.jsonl`, `localij/curated.jsonl`, `sets/curated.jsonl` and
-`regions/curated.jsonl`. Each has exactly the record type of its group, so a
+`traversal/curated.jsonl`, `localij/curated.jsonl`, `sets/curated.jsonl`,
+`regions/curated.jsonl` and `floats/curated.jsonl`. Each has exactly the
+record type of its group (for `floats`, that of `floats/cells.jsonl`), so a
 consumer runs it with the group's checks and nothing else. The expected
 values are the reference's answers, produced by the same generator as every
 other file; only the subjects are chosen. The tables below are the
@@ -796,6 +798,29 @@ degrees as it does for the sampled polygons.
 | `cliTriangle:7` | `polygon_test1.txt`, three vertices in San Francisco | inputfiles |
 | `cliDecagon:7` | `polygon_test2.txt`, ten vertices in San Francisco | inputfiles |
 
+### `floats/curated.jsonl`
+
+Cells on the antimeridian, with the record type of `floats/cells.jsonl`.
+The area of a cell is a sum over its edges of a term in the difference of
+the endpoints' longitudes; for an edge crossing ±180° the raw difference is
+close to ±360°, and the true, tiny step is lost to the rounding of a value
+that large unless both longitudes are first brought to the same side of the
+antimeridian. The seeded sample has straddling cells only at resolutions 0
+and 1, where the loss is far below tolerance; the defect is 0.1% at res 13
+and 9% at res 15.
+
+| Index | Source | Why |
+|---|---|---|
+| `8d056898db628bf`, `8d056898db62d7f` | tests | a res 13 straddling cell and a neighbour that does not straddle |
+| `8e32b605c51e067`, `8e32b605c51e15f` | tests | the same at res 14 |
+| `8f7f809893b5b9d`, `8f7f809893b5b9c` | tests | the same at res 15 |
+| `807ffffffffffff` … `8f7eb57221a2bb0` | hand | the cell containing (0°, 180°) at each resolution 0 through 15, in order |
+| `8023fffffffffff` … `8f2261a9d5b3701` | hand | the cell containing (45°, 180°) at each resolution |
+| `80ebfffffffffff` … `8feab634872d2b4` | hand | the cell containing (−60°, 180°) at each resolution |
+
+The *hand* rows are `latLngToCell` of the stated point at each resolution,
+54 subjects in all; the two res 0 cells are pentagons.
+
 ## Sampling
 
 Inputs are chosen by a deterministic procedure from `generator.seed`, so that

conformance/floats_test.go

@@ -56,13 +56,13 @@ func measure(value float64, err error) measurement {
 	return measurement{value: value, err: err}
 }
 
-// TestFloatCells runs every floats/cells record against h3.
+// TestFloatCells runs every floats/cells and floats/curated record against h3.
 func TestFloatCells(t *testing.T) {
 	t.Parallel()
 
 	tolerances := loadSuite(t).Tolerances
 
-	runGroup(t, "floats/cells", func(t *testing.T, record FloatCellRecord) {
+	check := func(t *testing.T, record FloatCellRecord) {
 		t.Helper()
 
 		cell := h3.CellFromString(record.Index)
@@ -91,7 +91,10 @@ func TestFloatCells(t *testing.T) {
 				t.Errorf("LatLngToCell(%v, %d) = %s (%v), want %s", point, record.InteriorRes, got, err, record.InteriorCells[i])
 			}
 		}
-	})
+	}
+
+	runGroup(t, "floats/cells", check)
+	runGroup(t, "floats/curated", check)
 }
 
 // TestFloatEdges runs every floats/edges record against h3.

conformance/internal/gen/gen.c

@@ -2175,7 +2175,7 @@ static void writeManifest(const char *outDir, const char *version, uint64_t seed
 /* ---- curated files (README, "Curated files") ---- */
 
 #define COUNT(array) (sizeof(array) / sizeof((array)[0]))
-#define NUM_CURATED_FILES 6
+#define NUM_CURATED_FILES 7
 #define MAX_CURATED_VERTS 16
 
 /* Indexes from the reference's own tests and command-line fixtures, plus
@@ -2189,6 +2189,31 @@ static const uint64_t CURATED_INSPECTION[] = {
     UINT64_C(0x8029fffffffffff),  UINT64_C(0x85283472fffffff),
 };
 
+/* Cells on the antimeridian: the three straddling and non-straddling pairs
+ * from the reference's cellArea tests, then the cell containing (lat, 180)
+ * at every resolution for lat 0, 45 and -60. Their areas depend on the
+ * reference subtracting edge longitudes on one side of +-180. */
+static const uint64_t CURATED_FLOATS[] = {
+    UINT64_C(0x8d056898db628bf), UINT64_C(0x8d056898db62d7f), UINT64_C(0x8e32b605c51e067),
+    UINT64_C(0x8e32b605c51e15f), UINT64_C(0x8f7f809893b5b9d), UINT64_C(0x8f7f809893b5b9c),
+    UINT64_C(0x807ffffffffffff),  UINT64_C(0x817ebffffffffff), UINT64_C(0x827eb7fffffffff),
+    UINT64_C(0x837eb5fffffffff),  UINT64_C(0x847eb57ffffffff), UINT64_C(0x857eb573fffffff),
+    UINT64_C(0x867eb5727ffffff),  UINT64_C(0x877eb5722ffffff), UINT64_C(0x887eb57221fffff),
+    UINT64_C(0x897eb57221bffff),  UINT64_C(0x8a7eb57221a7fff), UINT64_C(0x8b7eb57221a2fff),
+    UINT64_C(0x8c7eb57221a2bff),  UINT64_C(0x8d7eb57221a2bbf), UINT64_C(0x8e7eb57221a2bb7),
+    UINT64_C(0x8f7eb57221a2bb0),  UINT64_C(0x8023fffffffffff), UINT64_C(0x81227ffffffffff),
+    UINT64_C(0x822267fffffffff),  UINT64_C(0x832261fffffffff), UINT64_C(0x842261bffffffff),
+    UINT64_C(0x852261abfffffff),  UINT64_C(0x862261a9fffffff), UINT64_C(0x872261a9dffffff),
+    UINT64_C(0x882261a9d5fffff),  UINT64_C(0x892261a9d5bffff), UINT64_C(0x8a2261a9d5b7fff),
+    UINT64_C(0x8b2261a9d5b3fff),  UINT64_C(0x8c2261a9d5b37ff), UINT64_C(0x8d2261a9d5b373f),
+    UINT64_C(0x8e2261a9d5b3707),  UINT64_C(0x8f2261a9d5b3701), UINT64_C(0x80ebfffffffffff),
+    UINT64_C(0x81eabffffffffff),  UINT64_C(0x82eab7fffffffff), UINT64_C(0x83eab6fffffffff),
+    UINT64_C(0x84eab63ffffffff),  UINT64_C(0x85eab637fffffff), UINT64_C(0x86eab634fffffff),
+    UINT64_C(0x87eab6348ffffff),  UINT64_C(0x88eab63487fffff), UINT64_C(0x89eab634873ffff),
+    UINT64_C(0x8aeab634872ffff),  UINT64_C(0x8beab634872dfff), UINT64_C(0x8ceab634872d3ff),
+    UINT64_C(0x8deab634872d2bf),  UINT64_C(0x8eeab634872d2b7), UINT64_C(0x8feab634872d2b4),
+};
+
 static const uint64_t CURATED_HIERARCHY[] = {
     UINT64_C(0x88283080ddfffff), UINT64_C(0x8f283080dcb0ae2), UINT64_C(0x81083ffffffffff),
     UINT64_C(0x820807fffffffff), UINT64_C(0x85283473fffffff), UINT64_C(0x80c3fffffffffff),
@@ -2486,7 +2511,7 @@ static void writeCuratedPolygons(const char *outDir, FileInfo *info) {
     finishFile(out, path, info);
 }
 
-/* writeCuratedFiles writes the six curated files and their manifest entries. */
+/* writeCuratedFiles writes the seven curated files and their manifest entries. */
 static void writeCuratedFiles(const char *outDir, uint64_t seed, FileInfo *files) {
     files[0] = (FileInfo){"inspection/curated.jsonl", 0, ""};
     writeCuratedIndexes(outDir, seed, &files[0], CURATED_INSPECTION, COUNT(CURATED_INSPECTION),
@@ -2502,6 +2527,9 @@ static void writeCuratedFiles(const char *outDir, uint64_t seed, FileInfo *files
     writeCuratedSets(outDir, &files[4]);
     files[5] = (FileInfo){"regions/curated.jsonl", 0, ""};
     writeCuratedPolygons(outDir, &files[5]);
+    files[6] = (FileInfo){"floats/curated.jsonl", 0, ""};
+    writeCuratedIndexes(outDir, seed, &files[6], CURATED_FLOATS, COUNT(CURATED_FLOATS),
+                        writeFloatCellRecord);
 }
 
 /* Every record file of the suite, in manifest order. */
@@ -2514,7 +2542,7 @@ static const char *const SUITE_FILES[] = {
     "floats/resolutions.jsonl", "digests/resolutions.jsonl", "digests/baseCells.jsonl",
     "digests/patterns.jsonl",  "inspection/curated.jsonl",  "hierarchy/curated.jsonl",
     "traversal/curated.jsonl", "localij/curated.jsonl",     "sets/curated.jsonl",
-    "regions/curated.jsonl",
+    "regions/curated.jsonl",   "floats/curated.jsonl",
 };
 
 /* writeManifestFromDisk describes every suite file as it is on disk and

conformance/internal/gen/patches/0002-cellArea-antimeridian.patch

@@ -0,0 +1,33 @@
+Keep cell areas exact for cells crossing the antimeridian.
+
+Upstream change uber/h3#1274, not yet in a release. cagnoli() subtracts
+the raw longitudes of an edge's endpoints; for an edge crossing +-180 the
+difference is close to +-2*pi and the true (tiny) longitude step is lost
+to the rounding of a value that large. Straddling cells came out too
+small: about 0.1% at res 13 and 9% at res 15. Shifting both longitudes to
+the same side of the antimeridian before subtracting keeps the result
+exact, and sin and cos are 2*pi-periodic, so nothing else changes. The
+suite's floats/curated.jsonl holds straddling cells at every resolution
+that depend on both the reference and the port using this form.
+
+Delete this file once the release pinned in third_party/h3 contains the
+change; build.sh fails when a patch is already applied.
+
+--- a/src/h3lib/lib/area.c
++++ b/src/h3lib/lib/area.c
+@@ -37,7 +37,15 @@ static inline double cagnoli(LatLng x, LatLng y) {
+     double sa = sin(x.lat) * sin(y.lat);
+     double ca = cos(x.lat) * cos(y.lat);
+
++    // For an edge crossing the antimeridian, the raw difference is close to
++    // +-2*pi, and the true (tiny) longitude step is lost to the ulp of 2*pi.
++    // Shift both longitudes to the same side of the antimeridian first so
++    // that the subtraction is between two nearby values; sin and cos are
++    // 2*pi-periodic, so the result is mathematically unchanged.
+     double d = y.lng - x.lng;
++    if (d > M_PI || d < -M_PI) {
++        d = (y.lng - copysign(M_PI, y.lng)) - (x.lng - copysign(M_PI, x.lng));
++    }
+     double sd = sin(d);
+     double cd = cos(d);
+