Add S2Earth

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°)
+}