diff --git a/earth/earth.go b/earth/earth.go
index f101eeef..f55fa217 100644
--- a/earth/earth.go
+++ b/earth/earth.go
@@ -13,16 +13,13 @@
 // limitations under the License.
 
 /*
-Package earth implements functions for working with the planet Earth modeled as
-a sphere.
+Package earth provides constants for working with the planet Earth modeled as
+a sphere. Functions that operate on angles are in s1 (e.g., s1.EarthAngleFromLength),
+and functions that operate on s2 geometry types are in s2 (e.g., s2.EarthLengthFromLatLngs).
 */
 package earth
 
 import (
-	"math"
-
-	"github.com/golang/geo/s1"
-	"github.com/golang/geo/s2"
 	"github.com/google/go-units/unit"
 )
 
@@ -57,63 +54,3 @@ const (
 
 	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
-}
diff --git a/earth/earth_example_test.go b/earth/earth_example_test.go
deleted file mode 100644
index 7aadd6df..00000000
--- a/earth/earth_example_test.go
+++ /dev/null
@@ -1,82 +0,0 @@
-// 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°)
-}
diff --git a/earth/earth_test.go b/earth/earth_test.go
deleted file mode 100644
index 2cf7163b..00000000
--- a/earth/earth_test.go
+++ /dev/null
@@ -1,209 +0,0 @@
-// 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
-
-import (
-	"math"
-	"testing"
-
-	"github.com/golang/geo/s1"
-	"github.com/golang/geo/s2"
-	"github.com/google/go-units/unit"
-)
-
-func float64Eq(x, y float64) bool {
-	if x == y {
-		return true
-	}
-	if math.Abs(x) > math.Abs(y) {
-		return math.Abs(1-y/x) < 1e-14
-	}
-	return math.Abs(1-x/y) < 1e-14
-}
-
-var degreesToMeters = []struct {
-	angle  s1.Angle
-	length unit.Length
-}{
-	{-89.93201943346866 * s1.Degree, -1e7 * unit.Meter},
-	{-30 * s1.Degree, -3335853.035324518 * unit.Meter},
-	{0 * s1.Degree, 0 * unit.Meter},
-	{30 * s1.Degree, 3335853.035324518 * unit.Meter},
-	{89.93201943346866 * s1.Degree, 1e7 * unit.Meter},
-	{90 * s1.Degree, 10007559.105973555 * unit.Meter},
-	{179.86403886693734 * s1.Degree, 2e7 * unit.Meter},
-	{180 * s1.Degree, 20015118.21194711 * unit.Meter},
-	{359.72807773387467 * s1.Degree, 4e7 * unit.Meter},
-	{360 * s1.Degree, 40030236.42389422 * unit.Meter},
-	{899.3201943346867 * s1.Degree, 1e8 * unit.Meter},
-}
-
-func TestAngleFromLength(t *testing.T) {
-	for _, test := range degreesToMeters {
-		if got, want := AngleFromLength(test.length), test.angle; !float64Eq(got.Radians(), want.Radians()) {
-			t.Errorf("AngleFromLength(%v) = %v, want %v", test.length, got, want)
-		}
-	}
-}
-
-func TestLengthFromAngle(t *testing.T) {
-	for _, test := range degreesToMeters {
-		if got, want := LengthFromAngle(test.angle), test.length; !float64Eq(got.Meters(), want.Meters()) {
-			t.Errorf("LengthFromAngle(%v) = %v, want %v", test.angle, got, want)
-		}
-	}
-}
-
-func TestLengthFromPoints(t *testing.T) {
-	tests := []struct {
-		x1, y1, z1 float64
-		x2, y2, z2 float64
-		length     unit.Length
-	}{
-		{1, 0, 0, 1, 0, 0, 0 * unit.Meter},
-		{1, 0, 0, 0, 1, 0, 10007559.105973555 * unit.Meter},
-		{1, 0, 0, 0, 1, 1, 10007559.105973555 * unit.Meter},
-		{1, 0, 0, -1, 0, 0, 20015118.21194711 * unit.Meter},
-		{1, 2, 3, 2, 3, -1, 7680820.247060414 * unit.Meter},
-	}
-	for _, test := range tests {
-		p1 := s2.PointFromCoords(test.x1, test.y1, test.z1)
-		p2 := s2.PointFromCoords(test.x2, test.y2, test.z2)
-		if got, want := LengthFromPoints(p1, p2), test.length; !float64Eq(got.Meters(), want.Meters()) {
-			t.Errorf("LengthFromPoints(%v, %v) = %v, want %v", p1, p2, got, want)
-		}
-	}
-}
-
-func TestLengthFromLatLngs(t *testing.T) {
-	tests := []struct {
-		lat1, lng1, lat2, lng2 float64
-		length                 unit.Length
-	}{
-		{90, 0, 90, 0, 0 * unit.Meter},
-		{-37, 25, -66, -155, 8562022.790666264 * unit.Meter},
-		{0, 165, 0, -80, 12787436.635410652 * unit.Meter},
-		{47, -127, -47, 53, 20015118.077688109 * unit.Meter},
-		{51.961951, -180.227156, 51.782383, 181.126878, 95.0783566198074 * unit.Kilometer},
-	}
-	for _, test := range tests {
-		ll1 := s2.LatLngFromDegrees(test.lat1, test.lng1)
-		ll2 := s2.LatLngFromDegrees(test.lat2, test.lng2)
-		if got, want := LengthFromLatLngs(ll1, ll2), test.length; !float64Eq(got.Meters(), want.Meters()) {
-			t.Errorf("LengthFromLatLngs(%v, %v) = %v, want %v", ll1, ll2, got, want)
-		}
-	}
-}
-
-var (
-	earthArea        = unit.Area(Radius.Meters()*Radius.Meters()) * math.Pi * 4
-	steradiansToArea = []struct {
-		steradians float64
-		area       unit.Area
-	}{
-		{1, earthArea / 4 / math.Pi},
-		{4 * math.Pi, earthArea},
-		{s2.PolygonFromLoops([]*s2.Loop{s2.FullLoop()}).Area(), earthArea},
-		{s2.PolygonFromLoops([]*s2.Loop{s2.LoopFromPoints([]s2.Point{
-			s2.PointFromLatLng(s2.LatLngFromDegrees(-90, 0)),
-			s2.PointFromLatLng(s2.LatLngFromDegrees(0, 0)),
-			s2.PointFromLatLng(s2.LatLngFromDegrees(90, 0)),
-			s2.PointFromLatLng(s2.LatLngFromDegrees(0, -90)),
-		})}).Area(), earthArea / 4},
-		{s2.CellFromCellID(s2.CellIDFromFace(2)).ExactArea(), earthArea / 6},
-		{s2.AvgAreaMetric.Value(10), 81.07281893380302 * unit.SquareKilometer},  // average area of level 10 cells
-		{s2.AvgAreaMetric.Value(20), 77.31706517582228 * unit.SquareMeter},      // average area of level 20 cells
-		{s2.AvgAreaMetric.Value(30), 73.73529927808979 * unit.SquareMillimeter}, // average area of level 30 cells
-		{s2.PolygonFromLoops([]*s2.Loop{s2.EmptyLoop()}).Area(), 0 * unit.SquareMeter},
-	}
-)
-
-func TestAreaFromSteradians(t *testing.T) {
-	for _, test := range steradiansToArea {
-		if got, want := AreaFromSteradians(test.steradians), test.area; !float64Eq(got.SquareMeters(), want.SquareMeters()) {
-			t.Errorf("AreaFromSteradians(%v) = %v, want %v", test.steradians, got, want)
-		}
-	}
-}
-
-func TestSteradiansFromArea(t *testing.T) {
-	for _, test := range steradiansToArea {
-		if got, want := SteradiansFromArea(test.area), test.steradians; !float64Eq(got, want) {
-			t.Errorf("SteradiansFromArea(%v) = %v, want %v", test.area, got, want)
-		}
-	}
-}
-
-func TestInitialBearingFromLatLngs(t *testing.T) {
-	for _, tc := range []struct {
-		name string
-		a, b s2.LatLng
-		want s1.Angle
-	}{
-		{"Westward on equator", s2.LatLngFromDegrees(0, 50),
-			s2.LatLngFromDegrees(0, 100), s1.Degree * 90},
-		{"Eastward on equator", s2.LatLngFromDegrees(0, 50),
-			s2.LatLngFromDegrees(0, 0), s1.Degree * -90},
-		{"Northward on meridian", s2.LatLngFromDegrees(16, 28),
-			s2.LatLngFromDegrees(81, 28), s1.Degree * 0},
-		{"Southward on meridian", s2.LatLngFromDegrees(24, 64),
-			s2.LatLngFromDegrees(-27, 64), s1.Degree * 180},
-		{"Towards north pole", s2.LatLngFromDegrees(12, 76),
-			s2.LatLngFromDegrees(90, 50), s1.Degree * 0},
-		{"Towards south pole", s2.LatLngFromDegrees(-35, 105),
-			s2.LatLngFromDegrees(-90, -120), s1.Degree * 180},
-		{"Spain to Japan", s2.LatLngFromDegrees(40.4379332, -3.749576),
-			s2.LatLngFromDegrees(35.6733227, 139.6403486), s1.Degree * 29.2},
-		{"Japan to Spain", s2.LatLngFromDegrees(35.6733227, 139.6403486),
-			s2.LatLngFromDegrees(40.4379332, -3.749576), s1.Degree * -27.2},
-	} {
-		t.Run(tc.name, func(t *testing.T) {
-			got := InitialBearingFromLatLngs(tc.a, tc.b)
-			if diff := (got - tc.want).Abs(); diff > 0.01 {
-				t.Errorf("InitialBearingFromLatLngs(%s, %s): got %s, want %s, diff %s", tc.a, tc.b, got, tc.want, diff)
-			}
-		})
-	}
-}
-
-func TestInitialBearingFromLatLngsUndefinedResultDoesNotCrash(t *testing.T) {
-	// InitialBearingFromLatLngs says if a == b, a == -b, or a is one of Earth's
-	// poles, the return value is undefined.  Make sure it returns a real value
-	// (but don't assert what it is) rather than panicking or NaN.
-	// Bearing from a pole is undefined because 0° is north, but the observer
-	// can't face north from the north pole, so the calculation depends on the
-	// latitude value at the pole, even though 90°N 123°E and 90°N 45°W represent
-	// the same point.  Bearing is undefined when a == b because the observer can
-	// point any direction and still be present.  Bearing is undefined when
-	// a == -b (two antipodal points) because there are two possible paths.
-	for _, tc := range []struct {
-		name string
-		a, b s2.LatLng
-	}{
-		{"North pole prime meridian to Null Island", s2.LatLngFromDegrees(90, 0), s2.LatLngFromDegrees(0, 0)},
-		{"North pole facing east to Guatemala", s2.LatLngFromDegrees(90, 90), s2.LatLngFromDegrees(15, -90)},
-		{"South pole facing west to McMurdo", s2.LatLngFromDegrees(-90, -90), s2.LatLngFromDegrees(-78, 166)},
-		{"South pole anti-prime meridian to Null Island", s2.LatLngFromDegrees(-90, -180), s2.LatLngFromDegrees(0, 0)},
-		{"Jakarta and antipode", s2.LatLngFromDegrees(-6.109, 106.668), s2.LatLngFromDegrees(6.109, -180+106.668)},
-		{"Alert and antipode", s2.LatLngFromDegrees(82.499, -62.350), s2.LatLngFromDegrees(-82.499, 180-62.350)},
-	} {
-		t.Run(tc.name, func(t *testing.T) {
-			got := InitialBearingFromLatLngs(tc.a, tc.b)
-			if math.IsNaN(got.Radians()) || math.IsInf(got.Radians(), 0) {
-				t.Errorf("InitialBearingFromLatLngs(%s, %s): got %s, want a real value", tc.a, tc.b, got)
-			}
-		})
-	}
-}
diff --git a/s1/earth.go b/s1/earth.go
new file mode 100644
index 00000000..a40590f8
--- /dev/null
+++ b/s1/earth.go
@@ -0,0 +1,44 @@
+// 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 s1
+
+import (
+	"github.com/golang/geo/earth"
+	"github.com/google/go-units/unit"
+)
+
+// EarthAngleFromLength returns the angle subtended by a given arc length on
+// the spherical earth's surface.
+func EarthAngleFromLength(d unit.Length) Angle {
+	return Angle(float64(d/earth.Radius)) * Radian
+}
+
+// EarthLengthFromAngle returns the arc length on the spherical earth's surface
+// corresponding to a given angle.
+func EarthLengthFromAngle(a Angle) unit.Length {
+	return unit.Length(a.Radians()) * earth.Radius
+}
+
+// EarthAreaFromSteradians returns the surface area on the spherical earth
+// corresponding to a given solid angle in steradians.
+func EarthAreaFromSteradians(s float64) unit.Area {
+	return unit.Area(s * earth.Radius.Meters() * earth.Radius.Meters())
+}
+
+// EarthSteradiansFromArea returns the solid angle in steradians corresponding
+// to a given surface area on the spherical earth.
+func EarthSteradiansFromArea(a unit.Area) float64 {
+	return a.SquareMeters() / (earth.Radius.Meters() * earth.Radius.Meters())
+}
diff --git a/s1/earth_example_test.go b/s1/earth_example_test.go
new file mode 100644
index 00000000..586967b5
--- /dev/null
+++ b/s1/earth_example_test.go
@@ -0,0 +1,58 @@
+// 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 s1_test
+
+import (
+	"fmt"
+
+	"github.com/golang/geo/s1"
+	"github.com/golang/geo/s2"
+	"github.com/google/go-units/unit"
+)
+
+func ExampleEarthAngleFromLength() {
+	length := 500 * unit.Mile
+	angle := s1.EarthAngleFromLength(length)
+	fmt.Printf("I would walk 500 miles (%.4f rad)", angle.Radians())
+	// Output: I would walk 500 miles (0.1263 rad)
+}
+
+func ExampleEarthLengthFromAngle() {
+	angle := 2 * s1.Degree
+	length := s1.EarthLengthFromAngle(angle)
+	fmt.Printf("2° is %.0f miles", length.Miles())
+	// Output: 2° is 138 miles
+}
+
+func ExampleEarthAreaFromSteradians() {
+	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 := s1.EarthAreaFromSteradians(triangle.Area())
+	fmt.Printf("Bermuda Triangle is %.2f square miles", area.SquareMiles())
+	// Output: Bermuda Triangle is 464541.15 square miles
+}
+
+func ExampleEarthSteradiansFromArea() {
+	steradians := s1.EarthSteradiansFromArea(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
+}
diff --git a/s1/earth_test.go b/s1/earth_test.go
new file mode 100644
index 00000000..add269ef
--- /dev/null
+++ b/s1/earth_test.go
@@ -0,0 +1,106 @@
+// 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 s1
+
+import (
+	"math"
+	"testing"
+
+	"github.com/golang/geo/earth"
+	"github.com/google/go-units/unit"
+)
+
+func earthFloat64Eq(x, y float64) bool {
+	if x == y {
+		return true
+	}
+	if math.Abs(x) > math.Abs(y) {
+		return math.Abs(1-y/x) < 1e-14
+	}
+	return math.Abs(1-x/y) < 1e-14
+}
+
+var degreesToMeters = []struct {
+	angle  Angle
+	length unit.Length
+}{
+	{-89.93201943346866 * Degree, -1e7 * unit.Meter},
+	{-30 * Degree, -3335853.035324518 * unit.Meter},
+	{0 * Degree, 0 * unit.Meter},
+	{30 * Degree, 3335853.035324518 * unit.Meter},
+	{89.93201943346866 * Degree, 1e7 * unit.Meter},
+	{90 * Degree, 10007559.105973555 * unit.Meter},
+	{179.86403886693734 * Degree, 2e7 * unit.Meter},
+	{180 * Degree, 20015118.21194711 * unit.Meter},
+	{359.72807773387467 * Degree, 4e7 * unit.Meter},
+	{360 * Degree, 40030236.42389422 * unit.Meter},
+	{899.3201943346867 * Degree, 1e8 * unit.Meter},
+	{1 * Radian, earth.Radius},
+}
+
+func TestEarthAngleFromLength(t *testing.T) {
+	for _, test := range degreesToMeters {
+		got := EarthAngleFromLength(test.length)
+		want := test.angle
+		if !earthFloat64Eq(got.Radians(), want.Radians()) {
+			t.Errorf("EarthAngleFromLength(%v) = %v, want %v", test.length, got, want)
+		}
+	}
+}
+
+func TestEarthLengthFromAngle(t *testing.T) {
+	for _, test := range degreesToMeters {
+		got := EarthLengthFromAngle(test.angle)
+		want := test.length
+		if !earthFloat64Eq(got.Meters(), want.Meters()) {
+			t.Errorf("EarthLengthFromAngle(%v) = %v, want %v", test.angle, got, want)
+		}
+	}
+}
+
+var (
+	earthArea = unit.Area(earth.Radius.Meters()*earth.Radius.Meters()) * math.Pi * 4
+
+	steradiansToArea = []struct {
+		steradians float64
+		area       unit.Area
+	}{
+		{0, 0 * unit.SquareMeter},
+		{1, earthArea / 4 / math.Pi},
+		{math.Pi, earthArea / 4},
+		{2 * math.Pi, earthArea / 2},
+		{4 * math.Pi, earthArea},
+	}
+)
+
+func TestEarthAreaFromSteradians(t *testing.T) {
+	for _, test := range steradiansToArea {
+		got := EarthAreaFromSteradians(test.steradians)
+		want := test.area
+		if !earthFloat64Eq(got.SquareMeters(), want.SquareMeters()) {
+			t.Errorf("EarthAreaFromSteradians(%v) = %v, want %v", test.steradians, got, want)
+		}
+	}
+}
+
+func TestEarthSteradiansFromArea(t *testing.T) {
+	for _, test := range steradiansToArea {
+		got := EarthSteradiansFromArea(test.area)
+		want := test.steradians
+		if !earthFloat64Eq(got, want) {
+			t.Errorf("EarthSteradiansFromArea(%v) = %v, want %v", test.area, got, want)
+		}
+	}
+}
diff --git a/s2/contains_point_query_test.go b/s2/contains_point_query_test.go
index 0a3d82d3..5d6c4471 100644
--- a/s2/contains_point_query_test.go
+++ b/s2/contains_point_query_test.go
@@ -19,6 +19,7 @@ import (
 	"testing"
 
 	"github.com/golang/geo/s1"
+	"github.com/google/go-units/unit"
 )
 
 func TestContainsPointQueryVertexModelOpen(t *testing.T) {
@@ -138,7 +139,7 @@ func TestContainsPointQueryVertexModelClosed(t *testing.T) {
 
 func TestContainsPointQueryContainingShapes(t *testing.T) {
 	const numVerticesPerLoop = 10
-	maxLoopRadius := kmToAngle(10)
+	maxLoopRadius := s1.EarthAngleFromLength(10 * unit.Kilometer)
 	centerCap := CapFromCenterAngle(randomPoint(), maxLoopRadius)
 	index := NewShapeIndex()
 
diff --git a/s2/earth.go b/s2/earth.go
new file mode 100644
index 00000000..587c8906
--- /dev/null
+++ b/s2/earth.go
@@ -0,0 +1,57 @@
+// 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 s2
+
+import (
+	"math"
+
+	"github.com/golang/geo/s1"
+	"github.com/google/go-units/unit"
+)
+
+// EarthLengthFromPoints returns the distance between two points on the
+// spherical earth's surface.
+func EarthLengthFromPoints(a, b Point) unit.Length {
+	return s1.EarthLengthFromAngle(a.Distance(b))
+}
+
+// EarthLengthFromLatLngs returns the distance on the spherical earth's surface
+// between two LatLngs.
+func EarthLengthFromLatLngs(a, b LatLng) unit.Length {
+	return s1.EarthLengthFromAngle(a.Distance(b))
+}
+
+// EarthInitialBearingFromLatLngs 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 EarthInitialBearingFromLatLngs(a, b LatLng) s1.Angle {
+	latitudeA := a.Lat.Radians()
+	cosLatitudeB := math.Cos(b.Lat.Radians())
+	latitudeDiff := b.Lat.Radians() - a.Lat.Radians()
+	longitudeDiff := b.Lng.Radians() - a.Lng.Radians()
+
+	x := math.Sin(latitudeDiff) + math.Sin(latitudeA)*cosLatitudeB*2*haversine(longitudeDiff)
+	y := math.Sin(longitudeDiff) * cosLatitudeB
+	return s1.Angle(math.Atan2(y, x)) * s1.Radian
+}
+
+func haversine(radians float64) float64 {
+	sinHalf := math.Sin(radians / 2)
+	return sinHalf * sinHalf
+}
diff --git a/s2/earth_example_test.go b/s2/earth_example_test.go
new file mode 100644
index 00000000..a50bfab5
--- /dev/null
+++ b/s2/earth_example_test.go
@@ -0,0 +1,51 @@
+// 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 s2_test
+
+import (
+	"fmt"
+	"math"
+
+	"github.com/golang/geo/earth"
+	"github.com/golang/geo/s2"
+)
+
+func ExampleEarthLengthFromPoints() {
+	equator := s2.PointFromCoords(1, 0, 0)
+	pole := s2.PointFromCoords(0, 1, 0)
+	length := s2.EarthLengthFromPoints(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 ExampleEarthLengthFromLatLngs() {
+	chukchi := s2.LatLngFromDegrees(66.025893, -169.699684)
+	seward := s2.LatLngFromDegrees(65.609727, -168.093694)
+	length := s2.EarthLengthFromLatLngs(chukchi, seward)
+	fmt.Printf("Bering Strait is %.0f feet", length.Feet())
+	// Output: Bering Strait is 283979 feet
+}
+
+func ExampleEarthInitialBearingFromLatLngs() {
+	christchurch := s2.LatLngFromDegrees(-43.491402, 172.522275)
+	riogrande := s2.LatLngFromDegrees(-53.777156, -67.734719)
+	bearing := s2.EarthInitialBearingFromLatLngs(christchurch, riogrande)
+	fmt.Printf("Head southeast (%.2f°)", bearing.Degrees())
+	// Output: Head southeast (146.90°)
+}
diff --git a/s2/earth_test.go b/s2/earth_test.go
new file mode 100644
index 00000000..b4e5792c
--- /dev/null
+++ b/s2/earth_test.go
@@ -0,0 +1,202 @@
+// 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 s2
+
+import (
+	"math"
+	"testing"
+
+	"github.com/golang/geo/s1"
+	"github.com/google/go-units/unit"
+)
+
+func earthFloat64Eq(x, y float64) bool {
+	if x == y {
+		return true
+	}
+	if math.Abs(x) > math.Abs(y) {
+		return math.Abs(1-y/x) < 1e-14
+	}
+	return math.Abs(1-x/y) < 1e-14
+}
+
+func TestEarthLengthFromPoints(t *testing.T) {
+	tests := []struct {
+		x1, y1, z1 float64
+		x2, y2, z2 float64
+		length     unit.Length
+	}{
+		{1, 0, 0, 1, 0, 0, 0 * unit.Meter},
+		{1, 0, 0, 0, 1, 0, 10007559.105973555 * unit.Meter},
+		{1, 0, 0, 0, 1, 1, 10007559.105973555 * unit.Meter},
+		{1, 0, 0, -1, 0, 0, 20015118.21194711 * unit.Meter},
+		{1, 2, 3, 2, 3, -1, 7680820.247060414 * unit.Meter},
+	}
+	for _, test := range tests {
+		p1 := PointFromCoords(test.x1, test.y1, test.z1)
+		p2 := PointFromCoords(test.x2, test.y2, test.z2)
+		got := EarthLengthFromPoints(p1, p2)
+		want := test.length
+		if !earthFloat64Eq(got.Meters(), want.Meters()) {
+			t.Errorf("EarthLengthFromPoints(%v, %v) = %v, want %v", p1, p2, got, want)
+		}
+	}
+}
+
+func TestEarthLengthFromLatLngs(t *testing.T) {
+	tests := []struct {
+		lat1, lng1, lat2, lng2 float64
+		length                 unit.Length
+	}{
+		{90, 0, 90, 0, 0 * unit.Meter},
+		{-37, 25, -66, -155, 8562022.790666264 * unit.Meter},
+		{0, 165, 0, -80, 12787436.635410652 * unit.Meter},
+		{47, -127, -47, 53, 20015118.077688109 * unit.Meter},
+		{51.961951, -180.227156, 51.782383, 181.126878, 95.0783566198074 * unit.Kilometer},
+	}
+	for _, test := range tests {
+		ll1 := LatLngFromDegrees(test.lat1, test.lng1)
+		ll2 := LatLngFromDegrees(test.lat2, test.lng2)
+		got := EarthLengthFromLatLngs(ll1, ll2)
+		want := test.length
+		if !earthFloat64Eq(got.Meters(), want.Meters()) {
+			t.Errorf("EarthLengthFromLatLngs(%v, %v) = %v, want %v", ll1, ll2, got, want)
+		}
+	}
+}
+
+func TestEarthInitialBearingFromLatLngs(t *testing.T) {
+	for _, tc := range []struct {
+		name string
+		a, b LatLng
+		want s1.Angle
+	}{
+		{
+			"Westward on equator",
+			LatLngFromDegrees(0, 50),
+			LatLngFromDegrees(0, 100),
+			s1.Degree * 90,
+		},
+		{
+			"Eastward on equator",
+			LatLngFromDegrees(0, 50),
+			LatLngFromDegrees(0, 0),
+			s1.Degree * -90,
+		},
+		{
+			"Northward on meridian",
+			LatLngFromDegrees(16, 28),
+			LatLngFromDegrees(81, 28),
+			s1.Degree * 0,
+		},
+		{
+			"Southward on meridian",
+			LatLngFromDegrees(24, 64),
+			LatLngFromDegrees(-27, 64),
+			s1.Degree * 180,
+		},
+		{
+			"Towards north pole",
+			LatLngFromDegrees(12, 76),
+			LatLngFromDegrees(90, 50),
+			s1.Degree * 0,
+		},
+		{
+			"Towards south pole",
+			LatLngFromDegrees(-35, 105),
+			LatLngFromDegrees(-90, -120),
+			s1.Degree * 180,
+		},
+		{
+			"Spain to Japan",
+			LatLngFromDegrees(40.4379332, -3.749576),
+			LatLngFromDegrees(35.6733227, 139.6403486),
+			s1.Degree * 29.2,
+		},
+		{
+			"Japan to Spain",
+			LatLngFromDegrees(35.6733227, 139.6403486),
+			LatLngFromDegrees(40.4379332, -3.749576),
+			s1.Degree * -27.2,
+		},
+	} {
+		t.Run(tc.name, func(t *testing.T) {
+			got := EarthInitialBearingFromLatLngs(tc.a, tc.b)
+			if diff := (got - tc.want).Abs(); diff > 0.01 {
+				t.Errorf(
+					"EarthInitialBearingFromLatLngs(%s, %s): got %s, want %s, diff %s",
+					tc.a,
+					tc.b,
+					got,
+					tc.want,
+					diff,
+				)
+			}
+		})
+	}
+}
+
+func TestEarthInitialBearingFromLatLngsUndefinedResultDoesNotCrash(t *testing.T) {
+	// EarthInitialBearingFromLatLngs says if a == b, a == -b, or a is one of
+	// Earth's poles, the return value is undefined. Make sure it returns a real
+	// value (but don't assert what it is) rather than panicking or NaN.
+	for _, tc := range []struct {
+		name string
+		a, b LatLng
+	}{
+		{
+			"North pole prime meridian to Null Island",
+			LatLngFromDegrees(90, 0),
+			LatLngFromDegrees(0, 0),
+		},
+		{
+			"North pole facing east to Guatemala",
+			LatLngFromDegrees(90, 90),
+			LatLngFromDegrees(15, -90),
+		},
+		{
+			"South pole facing west to McMurdo",
+			LatLngFromDegrees(-90, -90),
+			LatLngFromDegrees(-78, 166),
+		},
+		{
+			"South pole anti-prime meridian to Null Island",
+			LatLngFromDegrees(-90, -180),
+			LatLngFromDegrees(0, 0),
+		},
+		{
+			"Jakarta and antipode",
+			LatLngFromDegrees(-6.109, 106.668),
+			LatLngFromDegrees(6.109, -180+106.668),
+		},
+		{
+			"Alert and antipode",
+			LatLngFromDegrees(82.499, -62.350),
+			LatLngFromDegrees(-82.499, 180-62.350),
+		},
+	} {
+		t.Run(tc.name, func(t *testing.T) {
+			got := EarthInitialBearingFromLatLngs(tc.a, tc.b)
+			if math.IsNaN(got.Radians()) || math.IsInf(got.Radians(), 0) {
+				t.Errorf(
+					"EarthInitialBearingFromLatLngs(%s, %s): got %s, want a real value",
+					tc.a,
+					tc.b,
+					got,
+				)
+			}
+		})
+	}
+}
diff --git a/s2/edge_distances_test.go b/s2/edge_distances_test.go
index a943c40e..d47774f9 100644
--- a/s2/edge_distances_test.go
+++ b/s2/edge_distances_test.go
@@ -20,6 +20,7 @@ import (
 
 	"github.com/golang/geo/r3"
 	"github.com/golang/geo/s1"
+	"github.com/google/go-units/unit"
 )
 
 func TestEdgeDistancesCheckDistance(t *testing.T) {
@@ -867,7 +868,7 @@ func TestEdgeDistancesEdgePairMaxDistance(t *testing.T) {
 func TestEdgeDistancesPointToLeft(t *testing.T) {
 	a := PointFromLatLng(LatLngFromDegrees(0, 0))
 	b := PointFromLatLng(LatLngFromDegrees(0, 5)) // east
-	dist := kmToAngle(10 / 1000.0)
+	dist := s1.EarthAngleFromLength(10 * unit.Meter)
 
 	c := PointToLeft(a, b, dist)
 	if got := a.Distance(c).Radians(); !float64Near(got, dist.Radians(), epsilon) {
@@ -882,7 +883,7 @@ func TestEdgeDistancesPointToLeft(t *testing.T) {
 func TestEdgeDistancesPointToRight(t *testing.T) {
 	a := PointFromLatLng(LatLngFromDegrees(0, 0))
 	b := PointFromLatLng(LatLngFromDegrees(0, 5)) // east
-	dist := kmToAngle(10 / 1000.0)
+	dist := s1.EarthAngleFromLength(10 * unit.Meter)
 
 	c := PointToRight(a, b, dist)
 	if got := a.Distance(c).Radians(); !float64Near(got, dist.Radians(), epsilon) {
diff --git a/s2/edge_query_closest_test.go b/s2/edge_query_closest_test.go
index a376290d..0b57ab97 100644
--- a/s2/edge_query_closest_test.go
+++ b/s2/edge_query_closest_test.go
@@ -20,6 +20,7 @@ import (
 
 	"github.com/golang/geo/r3"
 	"github.com/golang/geo/s1"
+	"github.com/google/go-units/unit"
 )
 
 func TestClosestEdgeQueryNoEdges(t *testing.T) {
@@ -238,7 +239,7 @@ func BenchmarkEdgeQueryFindEdges(b *testing.B) {
 				targetType:               queryTypePoint,
 				numTargetEdges:           0,
 				chooseTargetFromIndex:    false,
-				radiusKm:                 1000,
+				radius:                   1000 * unit.Kilometer,
 				maxDistanceFraction:      -1,
 				maxErrorFraction:         -1,
 				targetRadiusFraction:     0.0,
@@ -254,7 +255,7 @@ func BenchmarkEdgeQueryFindEdges(b *testing.B) {
 				targetType:               queryTypePoint,
 				numTargetEdges:           0,
 				chooseTargetFromIndex:    false,
-				radiusKm:                 1000,
+				radius:                   1000 * unit.Kilometer,
 				maxDistanceFraction:      -1,
 				maxErrorFraction:         -1,
 				targetRadiusFraction:     0.0,
@@ -271,7 +272,7 @@ func BenchmarkEdgeQueryFindEdges(b *testing.B) {
 				targetType:               queryTypePoint,
 				numTargetEdges:           0,
 				chooseTargetFromIndex:    false,
-				radiusKm:                 1000,
+				radius:                   1000 * unit.Kilometer,
 				maxDistanceFraction:      -1,
 				maxErrorFraction:         0.1,
 				targetRadiusFraction:     0.0,
@@ -288,7 +289,7 @@ func BenchmarkEdgeQueryFindEdges(b *testing.B) {
 				targetType:               queryTypePoint,
 				numTargetEdges:           0,
 				chooseTargetFromIndex:    false,
-				radiusKm:                 1000,
+				radius:                   1000 * unit.Kilometer,
 				maxDistanceFraction:      -1,
 				maxErrorFraction:         0.01,
 				targetRadiusFraction:     0.0,
@@ -302,7 +303,7 @@ func BenchmarkEdgeQueryFindEdges(b *testing.B) {
 				targetType:               queryTypeEdge,
 				numTargetEdges:           0,
 				chooseTargetFromIndex:    false,
-				radiusKm:                 1000,
+				radius:                   1000 * unit.Kilometer,
 				maxDistanceFraction:      -1,
 				maxErrorFraction:         -1,
 				targetRadiusFraction:     -1.0,
@@ -316,7 +317,7 @@ func BenchmarkEdgeQueryFindEdges(b *testing.B) {
 				targetType:               queryTypeCell,
 				numTargetEdges:           0,
 				chooseTargetFromIndex:    false,
-				radiusKm:                 1000,
+				radius:                   1000 * unit.Kilometer,
 				maxDistanceFraction:      -1,
 				maxErrorFraction:         -1,
 				targetRadiusFraction:     -1.0,
diff --git a/s2/edge_query_test.go b/s2/edge_query_test.go
index 66f82a55..02faa194 100644
--- a/s2/edge_query_test.go
+++ b/s2/edge_query_test.go
@@ -16,13 +16,13 @@ package s2
 
 import (
 	"flag"
-
 	"fmt"
 	"math/rand"
 	"reflect"
 	"testing"
 
 	"github.com/golang/geo/s1"
+	"github.com/google/go-units/unit"
 )
 
 var (
@@ -223,7 +223,7 @@ const edgeQueryTestNumEdges = 100
 const edgeQueryTestNumQueries = 200
 
 // The approximate radius of Cap from which query edges are chosen.
-var testCapRadius = kmToAngle(10)
+var testCapRadius = s1.EarthAngleFromLength(10 * unit.Kilometer)
 
 /*
 // testEdgeQueryWithGenerator is used to perform high volume random testing on EdqeQuery
@@ -340,12 +340,12 @@ func benchmarkEdgeQueryFindClosest(b *testing.B, bmOpts *edgeQueryBenchmarkOptio
 	index := NewShapeIndex()
 	opts := NewClosestEdgeQueryOptions().MaxResults(1).IncludeInteriors(bmOpts.includeInteriors)
 
-	radius := kmToAngle(bmOpts.radiusKm.Radians())
+	angle := s1.EarthAngleFromLength(bmOpts.radius)
 	if bmOpts.maxDistanceFraction > 0 {
-		opts.DistanceLimit(s1.ChordAngleFromAngle(s1.Angle(bmOpts.maxDistanceFraction) * radius))
+		opts.DistanceLimit(s1.ChordAngleFromAngle(s1.Angle(bmOpts.maxDistanceFraction) * angle))
 	}
 	if bmOpts.maxErrorFraction > 0 {
-		opts.MaxError(s1.ChordAngleFromAngle(s1.Angle(bmOpts.maxErrorFraction) * radius))
+		opts.MaxError(s1.ChordAngleFromAngle(s1.Angle(bmOpts.maxErrorFraction) * angle))
 	}
 
 	opts.UseBruteForce(*benchmarkBruteForce)
@@ -389,7 +389,7 @@ type edgeQueryBenchmarkOptions struct {
 	targetType               queryTargetType
 	numTargetEdges           int
 	chooseTargetFromIndex    bool
-	radiusKm                 s1.Angle
+	radius                   unit.Length
 	maxDistanceFraction      float64
 	maxErrorFraction         float64
 	targetRadiusFraction     float64
@@ -401,8 +401,8 @@ type edgeQueryBenchmarkOptions struct {
 // use in an edge query.
 //
 // Approximately numIndexEdges will be generated by the requested generator and
-// inserted. The geometry is generated within a Cap of the radius specified
-// by radiusKm (the index radius). Parameters with fraction in their
+// inserted. The geometry is generated within a Cap of the specified radius
+// (the index radius). Parameters with fraction in their
 // names are expressed as a fraction of this radius.
 //
 // Also generates a set of target geometries for the query, based on the
@@ -436,7 +436,10 @@ func generateEdgeQueryWithTargets(opts *edgeQueryBenchmarkOptions, query *EdgeQu
 	// Replace with r := rand.New(rand.NewSource(opts.randomSeed)) and pass through.
 	r := rand.New(rand.NewSource(opts.randomSeed))
 	opts.randomSeed++
-	indexCap := CapFromCenterAngle(randomPoint(r), opts.radiusKm)
+	indexCap := CapFromCenterAngle(
+		randomPoint(r),
+		s1.EarthAngleFromLength(opts.radius),
+	)
 
 	query.Reset()
 	queryIndex.Reset()
@@ -452,10 +455,11 @@ func generateEdgeQueryWithTargets(opts *edgeQueryBenchmarkOptions, query *EdgeQu
 	}
 
 	for i := 0; i < numTargets; i++ {
-		targetDist := fractionToRadius(opts.centerSeparationFraction, opts.radiusKm.Radians())
+		targetDist := fractionToAngle(opts.centerSeparationFraction, opts.radius)
 		targetCap := CapFromCenterAngle(
 			sampleBoundaryFromCap(CapFromCenterAngle(indexCap.Center(), targetDist)),
-			fractionToRadius(opts.targetRadiusFraction, opts.radiusKm.Radians()))
+			fractionToAngle(opts.targetRadiusFraction, opts.radius),
+		)
 
 		switch opts.targetType {
 		case queryTypePoint:
@@ -535,9 +539,9 @@ func sampleCellFromIndex(index *ShapeIndex) CellID {
 	return iter.CellID()
 }
 
-func fractionToRadius(fraction, radiusKm float64) s1.Angle {
+func fractionToAngle(fraction float64, radius unit.Length) s1.Angle {
 	if fraction < 0 {
 		fraction = -randomFloat64() * fraction
 	}
-	return s1.Angle(fraction) * kmToAngle(radiusKm)
+	return s1.EarthAngleFromLength(unit.Length(fraction) * radius)
 }
diff --git a/s2/loop_test.go b/s2/loop_test.go
index 7abfa21f..b5374d50 100644
--- a/s2/loop_test.go
+++ b/s2/loop_test.go
@@ -22,6 +22,7 @@ import (
 	"github.com/golang/geo/r1"
 	"github.com/golang/geo/r3"
 	"github.com/golang/geo/s1"
+	"github.com/google/go-units/unit"
 )
 
 var (
@@ -1864,7 +1865,8 @@ func BenchmarkLoopContainsPoint(b *testing.B) {
 				b.StopTimer()
 				loops := make([]*Loop, numLoopSamples)
 				for i := range numLoopSamples {
-					loops[i] = RegularLoop(randomPoint(), kmToAngle(10.0), vertices)
+					radius := s1.EarthAngleFromLength(10 * unit.Kilometer)
+					loops[i] = RegularLoop(randomPoint(), radius, vertices)
 				}
 
 				queries := make([][]Point, numLoopSamples)
diff --git a/s2/point_test.go b/s2/point_test.go
index 5da788a9..a768c68b 100644
--- a/s2/point_test.go
+++ b/s2/point_test.go
@@ -18,6 +18,7 @@ import (
 	"math"
 	"testing"
 
+	"github.com/golang/geo/earth"
 	"github.com/golang/geo/r3"
 	"github.com/golang/geo/s1"
 )
@@ -45,7 +46,7 @@ func TestOriginPoint(t *testing.T) {
 	}
 
 	// Check that the origin is not too close to either pole.
-	if dist := math.Acos(OriginPoint().Z) * earthRadiusKm; dist <= 50 {
+	if dist := math.Acos(OriginPoint().Z) * earth.Radius.Kilometers(); dist <= 50 {
 		t.Errorf("Origin point is to close to the North Pole. Got %v, want >= 50km", dist)
 	}
 }
diff --git a/s2/pointcompression_test.go b/s2/pointcompression_test.go
index 3a4707bd..5e6fad57 100644
--- a/s2/pointcompression_test.go
+++ b/s2/pointcompression_test.go
@@ -18,6 +18,9 @@ import (
 	"bytes"
 	"reflect"
 	"testing"
+
+	"github.com/golang/geo/s1"
+	"github.com/google/go-units/unit"
 )
 
 func TestFacesIterator(t *testing.T) {
@@ -33,9 +36,9 @@ func TestFacesIterator(t *testing.T) {
 }
 
 func makeSnappedPoints(nvertices int, level int) []Point {
-	const radiusKM = 0.1
 	center := PointFromCoords(1, 1, 1)
-	pts := regularPoints(center, kmToAngle(radiusKM), nvertices)
+	radius := s1.EarthAngleFromLength(0.1 * unit.Kilometer)
+	pts := regularPoints(center, radius, nvertices)
 	for i, pt := range pts {
 		id := CellFromPoint(pt).ID()
 		if level < id.Level() {
diff --git a/s2/polygon_test.go b/s2/polygon_test.go
index 52511a87..bbd3bc0d 100644
--- a/s2/polygon_test.go
+++ b/s2/polygon_test.go
@@ -21,6 +21,7 @@ import (
 	"testing"
 
 	"github.com/golang/geo/s1"
+	"github.com/google/go-units/unit"
 )
 
 const (
@@ -1144,7 +1145,7 @@ func TestPolygonInvert(t *testing.T) {
 	origin := PointFromLatLng(LatLngFromDegrees(0, 0))
 	pt := PointFromLatLng(LatLngFromDegrees(30, 30))
 	p := PolygonFromLoops([]*Loop{
-		RegularLoop(origin, 1000/earthRadiusKm, 100),
+		RegularLoop(origin, s1.EarthAngleFromLength(1000*unit.Kilometer), 100),
 	})
 
 	if p.ContainsPoint(pt) {
diff --git a/s2/predicates_test.go b/s2/predicates_test.go
index 31212427..c3a0758d 100644
--- a/s2/predicates_test.go
+++ b/s2/predicates_test.go
@@ -22,6 +22,7 @@ import (
 
 	"github.com/golang/geo/r3"
 	"github.com/golang/geo/s1"
+	"github.com/google/go-units/unit"
 )
 
 func TestPredicatesEpsilonForDigits(t *testing.T) {
@@ -305,7 +306,7 @@ func TestPredicatesStableSignFailureRate(t *testing.T) {
 	// that are as collinear as possible and spaced the given distance apart
 	// by counting up the times it returns Indeterminate.
 	failureCount := 0
-	m := math.Tan(spacing / earthRadiusKm)
+	m := math.Tan(s1.EarthAngleFromLength(unit.Length(spacing) * unit.Kilometer).Radians())
 	for range iters {
 		f := randomFrame()
 		a := f.col(0)
diff --git a/s2/regioncoverer_test.go b/s2/regioncoverer_test.go
index a7c86e56..55da8044 100644
--- a/s2/regioncoverer_test.go
+++ b/s2/regioncoverer_test.go
@@ -20,6 +20,9 @@ import (
 	"math/rand"
 	"reflect"
 	"testing"
+
+	"github.com/golang/geo/s1"
+	"github.com/google/go-units/unit"
 )
 
 func TestCovererRandomCells(t *testing.T) {
@@ -337,7 +340,8 @@ func BenchmarkRegionCovererCoveringLoop(b *testing.B) {
 			size := int(math.Pow(2.0, float64(n)))
 			regions := make([]Region, numCoveringBMRegions)
 			for i := range numCoveringBMRegions {
-				regions[i] = RegularLoop(randomPoint(), kmToAngle(10.0), size)
+				radius := s1.EarthAngleFromLength(10 * unit.Kilometer)
+				regions[i] = RegularLoop(randomPoint(), radius, size)
 			}
 			return regions
 		})
diff --git a/s2/s2_test.go b/s2/s2_test.go
index 6445f541..58b9da0d 100644
--- a/s2/s2_test.go
+++ b/s2/s2_test.go
@@ -74,14 +74,6 @@ func float64Near(x, y, ε float64) bool {
 	return math.Abs(x-y) <= ε
 }
 
-// The Earth's mean radius in kilometers (according to NASA).
-const earthRadiusKm = 6371.01
-
-// kmToAngle converts a distance on the Earth's surface to an angle.
-func kmToAngle(km float64) s1.Angle {
-	return s1.Angle(km / earthRadiusKm)
-}
-
 // randomBits returns a 64-bit random unsigned integer whose lowest "num" are random, and
 // whose other bits are zero.
 //
diff --git a/s2/s2_test_test.go b/s2/s2_test_test.go
index 4598fd98..492bb629 100644
--- a/s2/s2_test_test.go
+++ b/s2/s2_test_test.go
@@ -22,25 +22,6 @@ import (
 	"github.com/golang/geo/s1"
 )
 
-func TestKmToAngle(t *testing.T) {
-	tests := []struct {
-		have float64
-		want s1.Angle
-	}{
-		{0.0, 0.0},
-		{1.0, 0.00015696098420815537 * s1.Radian},
-		{earthRadiusKm, 1.0 * s1.Radian},
-		{-1.0, -0.00015696098420815537 * s1.Radian},
-		{-10000.0, -1.5696098420815536300 * s1.Radian},
-		{1e9, 156960.984208155363007 * s1.Radian},
-	}
-	for _, test := range tests {
-		if got := kmToAngle(test.have); !float64Eq(float64(got), float64(test.want)) {
-			t.Errorf("kmToAngle(%f) = %0.20f, want %0.20f", test.have, got, test.want)
-		}
-	}
-}
-
 func numVerticesAtLevel(level int) int {
 	// Sanity / overflow check
 	if level < 0 || level > 14 {