Add Equal Earth projection, beginning basic value checking of project…

…ions (#7)
owlpinetech · Jan 20, 2024 · 0423275 · 0423275
1 parent 8543683
commit 0423275
Show file tree

Hide file tree

Showing 5 changed files with 137 additions and 6 deletions.
diff --git a/bounded_test.go b/bounded_test.go
@@ -94,11 +94,15 @@ func FuzzOrthographicProjectBounded(f *testing.F) {
 }
 
 func FuzzRobinsonProjectBounded(f *testing.F) {
-	projectionBoundedFuzz(f, NewRobinsonProjection())
+	projectionBoundedFuzz(f, NewRobinson())
 }
 
 func FuzzNaturalEarthProjectBounded(f *testing.F) {
-	projectionBoundedFuzz(f, NewNaturalEarthProjection())
+	projectionBoundedFuzz(f, NewNaturalEarth())
+}
+
+func FuzzEqualEarthProjectBounded(f *testing.F) {
+	projectionBoundedFuzz(f, NewEqualEarth())
 }
 
 func projectionBoundedFuzz(f *testing.F, proj Projection) {

diff --git a/invertability_test.go b/invertability_test.go
@@ -86,13 +86,17 @@ func FuzzSinusoidalProjectInverse(f *testing.F) {
 //}
 
 func FuzzRobinsonProjectInverse(f *testing.F) {
-	projectInverseFuzz(f, NewRobinsonProjection())
+	projectInverseFuzz(f, NewRobinson())
 }
 
 func FuzzNaturalEarthProjectInverse(f *testing.F) {
-	projectInverseFuzz(f, NewNaturalEarthProjection())
+	projectInverseFuzz(f, NewNaturalEarth())
 }
 
+//func FuzzEqualEarthProjectInverse(f *testing.F) {
+//	projectInverseFuzz(f, NewEqualEarth())
+//}
+
 func withinTolerance(n1, n2, tolerance float64) bool {
 	if n1 == n2 {
 		return true

diff --git a/project_test.go b/project_test.go
@@ -0,0 +1,88 @@
+package flatsphere
+
+import (
+	"fmt"
+	"math"
+	"testing"
+)
+
+type projectTestCase struct {
+	lat     float64
+	lon     float64
+	expectX float64
+	expectY float64
+}
+
+type inverseTestCase struct {
+	x         float64
+	y         float64
+	expectLat float64
+	expectLon float64
+}
+
+func checkProject(t *testing.T, name string, proj Projection, cases []projectTestCase) {
+	for ind, tc := range cases {
+		t.Run(fmt.Sprintf("%s%d", name, ind), func(t *testing.T) {
+			x, y := proj.Project(tc.lat, tc.lon)
+			if !withinTolerance(x, tc.expectX, 0.000001) || !withinTolerance(y, tc.expectY, 0.000001) {
+				t.Errorf("projection at %f,%f expected %e,%e, but got %e,%e", tc.lat, tc.lon, tc.expectX, tc.expectY, x, y)
+			}
+		})
+	}
+}
+
+func checkInverse(t *testing.T, name string, proj Projection, cases []inverseTestCase) {
+	for ind, tc := range cases {
+		t.Run(fmt.Sprintf("%s%d", name, ind), func(t *testing.T) {
+			lat, lon := proj.Inverse(tc.x, tc.y)
+			if !withinTolerance(lat, tc.expectLat, 0.000001) || !withinTolerance(lon, tc.expectLon, 0.000001) {
+				t.Errorf("inverse at %f,%f expected %e,%e, but got %e,%e", tc.x, tc.y, tc.expectLat, tc.expectLon, lat, lon)
+			}
+		})
+	}
+}
+
+func TestMercatorProjectSanity(t *testing.T) {
+	checkProject(t, "mercator", NewMercator(), []projectTestCase{
+		{0, 0, 0, 0},
+		{math.Pi / 4, math.Pi / 2, math.Pi / 2, math.Log(math.Tan(math.Pi/4) + 1/math.Cos(math.Pi/4))},
+		{-math.Pi / 4, math.Pi / 2, math.Pi / 2, math.Log(math.Tan(-math.Pi/4) + 1/math.Cos(-math.Pi/4))},
+		{math.Pi / 4, -math.Pi / 2, -math.Pi / 2, math.Log(math.Tan(math.Pi/4) + 1/math.Cos(math.Pi/4))},
+	})
+}
+
+func TestMercatorInverseSanity(t *testing.T) {
+	checkInverse(t, "invMercator", NewMercator(), []inverseTestCase{
+		{0, 0, 0, 0},
+		{math.Pi / 2, math.Log(math.Tan(math.Pi/4) + 1/math.Cos(math.Pi/4)), math.Pi / 4, math.Pi / 2},
+		{math.Pi / 2, math.Log(math.Tan(-math.Pi/4) + 1/math.Cos(-math.Pi/4)), -math.Pi / 4, math.Pi / 2},
+		{-math.Pi / 2, math.Log(math.Tan(math.Pi/4) + 1/math.Cos(math.Pi/4)), math.Pi / 4, -math.Pi / 2},
+	})
+}
+
+/*func TestTransverseMercatorProjectSanity(t *testing.T) {
+	xFrom := func(lat, lon float64) float64 {
+		return math.Log((1+math.Sin(lon)*math.Cos(lat))/(1-math.Sin(lon)*math.Cos(lat))) / 2
+	}
+	yFrom := func(lat, lon float64) float64 {
+		return math.Atan(math.Tan(lat) / math.Cos(lon))
+	}
+	checkProject(t, "transverseMercator", NewObliqueProjection(NewMercator(), 0, math.Pi/2, -math.Pi/2), []projectTestCase{
+		{0, 0, 0, 0},
+		{math.Pi / 2, 0, 0, math.Pi / 2},
+		{math.Pi / 4, 0, 0, math.Pi / 4},
+		{-math.Pi / 2, 0, 0, -math.Pi / 2},
+		{math.Pi / 4, math.Pi / 2, xFrom(math.Pi/4, math.Pi/2), yFrom(math.Pi/4, math.Pi/2)},
+		{-math.Pi / 4, math.Pi / 2, math.Pi / 2, -math.Log(math.Tan(-math.Pi/4) + 1/math.Cos(-math.Pi/4))},
+		{math.Pi / 4, -math.Pi / 2, -math.Pi / 2, -math.Log(math.Tan(math.Pi/4) + 1/math.Cos(math.Pi/4))},
+	})
+}
+
+func TestTransverseMercatorInverseSanity(t *testing.T) {
+	checkInverse(t, "invTransverseMercator", NewObliqueProjection(NewMercator(), 0, math.Pi/2, -math.Pi/2), []inverseTestCase{
+		{0, 0, 0, 0},
+		{-math.Pi / 2, math.Log(math.Tan(math.Pi/4) + 1/math.Cos(math.Pi/4)), math.Pi / 4, math.Pi / 2},
+		{math.Pi / 2, -math.Log(math.Tan(-math.Pi/4) + 1/math.Cos(-math.Pi/4)), -math.Pi / 4, math.Pi / 2},
+		{-math.Pi / 2, -math.Log(math.Tan(math.Pi/4) + 1/math.Cos(math.Pi/4)), math.Pi / 4, -math.Pi / 2},
+	})
+}*/
diff --git a/pseudocylindrical.go b/pseudocylindrical.go
@@ -129,3 +129,38 @@ func (e EckertIV) Inverse(x float64, y float64) (float64, float64) {
 func (e EckertIV) PlanarBounds() Bounds {
 	return NewRectangleBounds(4, 2)
 }
+
+// An equal-area pseudocylindrical projection.
+// https://en.wikipedia.org/wiki/Equal_Earth_projection
+type EqualEarth struct{}
+
+var (
+	eeB      float64 = math.Sqrt(3) / 2
+	eeYscale float64 = equalEarthPoly(math.Pi/3) / (math.Pi / 3)
+)
+
+func NewEqualEarth() EqualEarth {
+	return EqualEarth{}
+}
+
+func equalEarthPoly(x float64) float64 {
+	return 0.003796*math.Pow(x, 9) + 0.000893*math.Pow(x, 7) - 0.081106*math.Pow(x, 3) + 1.340264*x
+}
+
+func equalEarthDeriv(x float64) float64 {
+	return 9*0.003796*math.Pow(x, 8) + 7*0.000893*math.Pow(x, 6) - 3*0.081106*math.Pow(x, 2) + 1.340264
+}
+
+func (e EqualEarth) Project(latitude float64, longitude float64) (float64, float64) {
+	theta := math.Asin(math.Sqrt(3) / 2 * math.Sin(latitude))
+	return math.Cos(theta) / eeB / equalEarthDeriv(theta) * longitude, equalEarthPoly(theta)
+}
+
+func (e EqualEarth) Inverse(x float64, y float64) (float64, float64) {
+	theta := newtonsMethod(y/eeYscale, equalEarthPoly, equalEarthDeriv, 1e-6, 1e-15, 125)
+	return math.Asin(math.Sin(theta) / eeB), x * eeB / math.Cos(theta) * equalEarthDeriv(theta)
+}
+
+func (e EqualEarth) PlanarBounds() Bounds {
+	return NewRectangleBounds(2*math.Pi, math.Pi)
+}
diff --git a/tabular.go b/tabular.go
@@ -42,13 +42,13 @@ var (
 
 // Create a new Robinson projection, a well-known instance of a pseudocylindrical projection defined by a table of values.
 // https://en.wikipedia.org/wiki/Robinson_projection
-func NewRobinsonProjection() TabularProjection {
+func NewRobinson() TabularProjection {
 	return NewTabularProjection(robinsonNaturalEarthLatitudes, robinsonLengthRatios, robinsonDistRatios, 4, 0.5072)
 }
 
 // Create a new Natural Earth projection, a well-known instance of a pseudocylindrical projection defined by a table of values.
 // https://en.wikipedia.org/wiki/Natural_Earth_projection
-func NewNaturalEarthProjection() TabularProjection {
+func NewNaturalEarth() TabularProjection {
 	return NewTabularProjection(robinsonNaturalEarthLatitudes, naturalEarthLengthRatios, naturalEarthDistRatios, 4, 0.520)
 }