From 04232755f9f54fd534ec6545379e18a7b9819751 Mon Sep 17 00:00:00 2001 From: Robert Kleffner Date: Sat, 20 Jan 2024 12:29:10 -0700 Subject: [PATCH] Add Equal Earth projection, beginning basic value checking of projections (#7) --- bounded_test.go | 8 +++- invertability_test.go | 8 +++- project_test.go | 88 +++++++++++++++++++++++++++++++++++++++++++ pseudocylindrical.go | 35 +++++++++++++++++ tabular.go | 4 +- 5 files changed, 137 insertions(+), 6 deletions(-) create mode 100644 project_test.go diff --git a/bounded_test.go b/bounded_test.go index 9456767..924cd27 100644 --- 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 index f72db01..12421a9 100644 --- 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 new file mode 100644 index 0000000..88f19f4 --- /dev/null +++ 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 index b0d620d..86a7760 100644 --- 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 index 6bba6ab..9fd4b5c 100644 --- 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) }