Add S2Earth (#191)

This is a copy of the google3 implementation, which was not included in
previously because it lived elsewhere in the source tree. References to
other files in google3 have been removed from comments.

Fixes #151

Author: flwyd

earth/earth.go

@@ -0,0 +1,119 @@
+// Copyright 2025 Google LLC
+//
+// 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 earth implements functions for working with the planet Earth modeled as
+a sphere.
+*/
+package earth
+
+import (
+	"math"
+
+	"github.com/golang/geo/s1"
+	"github.com/golang/geo/s2"
+	"github.com/google/go-units/unit"
+)
+
+const (
+	// Radius is the Earth's mean radius, which is the radius of the
+	// equivalent sphere with the same surface area. According to NASA,
+	// this value is 6371.01 +/- 0.02 km. The equatorial radius is 6378.136
+	// km, and the polar radius is 6356.752 km. They differ by one part in
+	// 298.257.
+	//
+	// Reference: http://ssd.jpl.nasa.gov/phys_props_earth.html, which
+	// quotes Yoder, C.F. 1995. "Astrometric and Geodetic Properties of
+	// Earth and the Solar System" in Global Earth Physics, A Handbook of
+	// Physical Constants, AGU Reference Shelf 1, American Geophysical
+	// Union, Table 2.
+	//
+	// This value is the same as in s2earth.h and S2Earth.java in order to be
+	// able to make consistent conversions across programming languages.
+	Radius = 6371.01 * unit.Kilometer
+
+	// LowestAltitude is the altitude of the lowest known point on Earth,
+	// the Challenger Deep, below the surface of the spherical Earth. This value
+	// is the same as the C++ and Java implementations of S2Earth. The value may
+	// change as more precise measurements are made.
+
+	LowestAltitude = -10898 * unit.Meter
+
+	// HighestAltitude is the altitude of the highest known point on Earth,
+	// Mount Everest, above the surface of the spherical Earth. This value is the
+	// same as the C++ and Java implementations of S2Earth. The value may change
+	// as more precise measurements are made.
+
+	HighestAltitude = 8848 * unit.Meter
+)
+
+// AngleFromLength returns the angle from a given distance on the spherical
+// earth's surface.
+func AngleFromLength(d unit.Length) s1.Angle {
+	return s1.Angle(float64(d/Radius)) * s1.Radian
+}
+
+// LengthFromAngle returns the distance on the spherical earth's surface from
+// a given angle.
+func LengthFromAngle(a s1.Angle) unit.Length {
+	return unit.Length(a.Radians()) * Radius
+}
+
+// LengthFromPoints returns the distance between two points on the spherical
+// earth's surface.
+func LengthFromPoints(a, b s2.Point) unit.Length {
+	return LengthFromAngle(a.Distance(b))
+}
+
+// LengthFromLatLngs returns the distance on the spherical earth's surface
+// between two LatLngs.
+func LengthFromLatLngs(a, b s2.LatLng) unit.Length {
+	return LengthFromAngle(a.Distance(b))
+}
+
+// AreaFromSteradians returns the area on the spherical Earth's surface covered
+// by s steradians, as returned by Area() methods on s2 geometry types.
+func AreaFromSteradians(s float64) unit.Area {
+	return unit.Area(s * Radius.Meters() * Radius.Meters())
+}
+
+// SteradiansFromArea returns the number of steradians covered by an area on the
+// spherical Earth's surface.  The value will be between 0 and 4 * math.Pi if a
+// does not exceed the area of the Earth.
+func SteradiansFromArea(a unit.Area) float64 {
+	return a.SquareMeters() / (Radius.Meters() * Radius.Meters())
+}
+
+// InitialBearingFromLatLngs computes the initial bearing from a to b.
+//
+// This is the bearing an observer at point a has when facing point b. A bearing
+// of 0 degrees is north, and it increases clockwise (90 degrees is east, etc).
+//
+// If a == b, a == -b, or a is one of the Earth's poles, the return value is
+// undefined.
+func InitialBearingFromLatLngs(a, b s2.LatLng) s1.Angle {
+	lat1 := a.Lat.Radians()
+	cosLat2 := math.Cos(b.Lat.Radians())
+	latDiff := b.Lat.Radians() - a.Lat.Radians()
+	lngDiff := b.Lng.Radians() - a.Lng.Radians()
+
+	x := math.Sin(latDiff) + math.Sin(lat1)*cosLat2*2*haversine(lngDiff)
+	y := math.Sin(lngDiff) * cosLat2
+	return s1.Angle(math.Atan2(y, x)) * s1.Radian
+}
+
+func haversine(radians float64) float64 {
+	sinHalf := math.Sin(radians / 2)
+	return sinHalf * sinHalf
+}

earth/earth_example_test.go

@@ -0,0 +1,82 @@
+// Copyright 2025 Google LLC
+//
+// 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 earth_test
+
+import (
+	"fmt"
+	"math"
+
+	"github.com/golang/geo/earth"
+	"github.com/golang/geo/s1"
+	"github.com/golang/geo/s2"
+	"github.com/google/go-units/unit"
+)
+
+func ExampleAngleFromLength() {
+	length := 500 * unit.Mile
+	angle := earth.AngleFromLength(length)
+	fmt.Printf("I would walk 500 miles (%.4f rad)", angle.Radians())
+	// Output: I would walk 500 miles (0.1263 rad)
+}
+
+func ExampleLengthFromAngle() {
+	angle := 2 * s1.Degree
+	length := earth.LengthFromAngle(angle)
+	fmt.Printf("2° is %.0f miles", length.Miles())
+	// Output: 2° is 138 miles
+}
+
+func ExampleLengthFromPoints() {
+	equator := s2.PointFromCoords(1, 0, 0)
+	pole := s2.PointFromCoords(0, 1, 0)
+	length := earth.LengthFromPoints(equator, pole)
+	fmt.Printf("Equator to pole is %.2f km, π/2*Earth radius is %.2f km",
+		length.Kilometers(), math.Pi/2*earth.Radius.Kilometers())
+	// Output: Equator to pole is 10007.56 km, π/2*Earth radius is 10007.56 km
+}
+
+func ExampleLengthFromLatLngs() {
+	chukchi := s2.LatLngFromDegrees(66.025893, -169.699684)
+	seward := s2.LatLngFromDegrees(65.609727, -168.093694)
+	length := earth.LengthFromLatLngs(chukchi, seward)
+	fmt.Printf("Bering Strait is %.0f feet", length.Feet())
+	// Output: Bering Strait is 283979 feet
+}
+
+func ExampleAreaFromSteradians() {
+	bermuda := s2.PointFromLatLng(s2.LatLngFromDegrees(32.361457, -64.663495))
+	culebra := s2.PointFromLatLng(s2.LatLngFromDegrees(18.311199, -65.300765))
+	miami := s2.PointFromLatLng(s2.LatLngFromDegrees(25.802018, -80.269892))
+	triangle := s2.PolygonFromLoops(
+		[]*s2.Loop{s2.LoopFromPoints([]s2.Point{bermuda, miami, culebra})})
+	area := earth.AreaFromSteradians(triangle.Area())
+	fmt.Printf("Bermuda Triangle is %.2f square miles", area.SquareMiles())
+	// Output: Bermuda Triangle is 464541.15 square miles
+}
+
+func ExampleSteradiansFromArea() {
+	steradians := earth.SteradiansFromArea(unit.SquareCentimeter)
+	fmt.Printf("1 square centimeter is %g steradians, close to a level %d cell",
+		steradians, s2.AvgAreaMetric.ClosestLevel(steradians))
+	// Output: 1 square centimeter is 2.4636750563592804e-18 steradians, close to a level 30 cell
+}
+
+func ExampleInitialBearingFromLatLngs() {
+	christchurch := s2.LatLngFromDegrees(-43.491402, 172.522275)
+	riogrande := s2.LatLngFromDegrees(-53.777156, -67.734719)
+	bearing := earth.InitialBearingFromLatLngs(christchurch, riogrande)
+	fmt.Printf("Head southeast (%.2f°)", bearing.Degrees())
+	// Output: Head southeast (146.90°)
+}