From a5a278ddea2ce230395d2d406a3b921558f90aa2 Mon Sep 17 00:00:00 2001 From: Robert Kleffner Date: Wed, 17 Jan 2024 20:50:16 -0700 Subject: [PATCH] Add an important new generalized test category, and fix a couple bugs with it (#2) Verify projected points are within the projection's planar bounds, add Eckert IV projection, fix bounds for HEALPix --- bounded_test.go | 111 ++++++++++++++++++++++++++++++++++++++++++ healpix.go | 4 +- invertability_test.go | 28 +++++------ pseudocylindrical.go | 28 +++++++++++ 4 files changed, 155 insertions(+), 16 deletions(-) create mode 100644 bounded_test.go diff --git a/bounded_test.go b/bounded_test.go new file mode 100644 index 0000000..d8c42e0 --- /dev/null +++ b/bounded_test.go @@ -0,0 +1,111 @@ +package flatsphere + +import ( + "math" + "testing" +) + +func FuzzMercatorProjectBounded(f *testing.F) { + projectionBoundedFuzz(f, NewMercator()) +} + +func FuzzPlateCarreeProjectBounded(f *testing.F) { + projectionBoundedFuzz(f, NewPlateCarree()) +} + +func FuzzEquirectangularPositiveProjectBounded(f *testing.F) { + projectionBoundedFuzz(f, NewEquirectangular(45*math.Pi/180)) +} + +func FuzzEquirectangularNegativeProjectBounded(f *testing.F) { + projectionBoundedFuzz(f, NewEquirectangular(-45*math.Pi/180)) +} + +func FuzzLambertCylindricalProjectBounded(f *testing.F) { + projectionBoundedFuzz(f, NewCylindricalEqualArea(0)) +} + +func FuzzBehrmannProjectBounded(f *testing.F) { + projectionBoundedFuzz(f, NewBehrmann()) +} + +func FuzzGallOrthographicProjectBounded(f *testing.F) { + projectionBoundedFuzz(f, NewGallOrthographic()) +} + +func FuzzHoboDyerProjectBounded(f *testing.F) { + projectionBoundedFuzz(f, NewHoboDyer()) +} + +func FuzzGallStereographicProjectBounded(f *testing.F) { + projectionBoundedFuzz(f, NewGallStereographic()) +} + +func FuzzMillerProjectBounded(f *testing.F) { + projectionBoundedFuzz(f, NewMiller()) +} + +func FuzzCentralProjectBounded(f *testing.F) { + projectionBoundedFuzz(f, NewCentral()) +} + +func FuzzSinusoidalProjectBounded(f *testing.F) { + projectionBoundedFuzz(f, NewSinusoidal()) +} + +func FuzzHEALPixStandardProjectBounded(f *testing.F) { + projectionBoundedFuzz(f, NewHEALPixStandard()) +} + +func FuzzMollweideProjectBounded(f *testing.F) { + projectionBoundedFuzz(f, NewMollweide()) +} + +func FuzzHomolosineProjectBounded(f *testing.F) { + projectionBoundedFuzz(f, NewHomolosine()) +} + +func FuzzEckertIVProjectBounded(f *testing.F) { + projectionBoundedFuzz(f, NewEckertIV()) +} + +func FuzzStereographicProjectBounded(f *testing.F) { + projectionBoundedFuzz(f, NewStereographic()) +} + +func FuzzPolarProjectBounded(f *testing.F) { + projectionBoundedFuzz(f, NewPolar()) +} + +func FuzzLambertAzimuthalProjectBounded(f *testing.F) { + projectionBoundedFuzz(f, NewLambertAzimuthal()) +} + +//func FuzzTransverseMercatorProjectBounded(f *testing.F) { +// projectionBoundedFuzz(f, NewObliqueProjection(NewMercator(), 0, math.Pi/2, -math.Pi/2)) +//} + +//func FuzzGnomonicProjectBounded(f *testing.F) { +// projectionBoundedFuzz(f, NewGnomonic()) +//} + +func FuzzOrthographicProjectBounded(f *testing.F) { + projectionBoundedFuzz(f, NewOrthographic()) +} + +func projectionBoundedFuzz(f *testing.F, proj Projection) { + f.Add(0.0, 0.0) + f.Add(0.0, math.Pi) + f.Add(math.Pi/2, math.Pi/4) + f.Add(math.Pi/2, 0.0) + f.Add(-math.Pi/2, -math.Pi/4) + f.Add(math.Pi/2, math.Pi) + f.Fuzz(func(t *testing.T, lat float64, lon float64) { + lat = math.Mod(lat, math.Pi/2) + lon = math.Mod(lon, math.Pi) + x, y := proj.Project(lat, lon) + if !proj.PlanarBounds().Within(x, y) { + t.Errorf("projected point for %e,%e (equal to %e,%e) was outside of the planar bounds for %v", lat, lon, x, y, proj) + } + }) +} diff --git a/healpix.go b/healpix.go index 8ba6047..73f6b0e 100644 --- a/healpix.go +++ b/healpix.go @@ -54,8 +54,8 @@ func (h HEALPixStandard) Inverse(x float64, y float64) (float64, float64) { func (h HEALPixStandard) PlanarBounds() Bounds { return Bounds{ - XMin: 0, - XMax: 2 * math.Pi, + XMin: -math.Pi, + XMax: math.Pi, YMin: -math.Pi / 2, YMax: math.Pi / 2, } diff --git a/invertability_test.go b/invertability_test.go index 4693f6a..b12e10d 100644 --- a/invertability_test.go +++ b/invertability_test.go @@ -21,59 +21,59 @@ func FuzzEquirectangularNegativeProjectInverse(f *testing.F) { projectInverseFuzz(f, NewEquirectangular(-45*math.Pi/180)) } -func FuzzLambertCylindrical(f *testing.F) { +func FuzzLambertCylindricalProjectInverse(f *testing.F) { projectInverseFuzz(f, NewCylindricalEqualArea(0)) } -func FuzzBehrmann(f *testing.F) { +func FuzzBehrmannProjectInverse(f *testing.F) { projectInverseFuzz(f, NewBehrmann()) } -func FuzzGallOrthographic(f *testing.F) { +func FuzzGallOrthographicProjectInverse(f *testing.F) { projectInverseFuzz(f, NewGallOrthographic()) } -func FuzzHoboDyer(f *testing.F) { +func FuzzHoboDyerProjectInverse(f *testing.F) { projectInverseFuzz(f, NewHoboDyer()) } -func FuzzGallStereographic(f *testing.F) { +func FuzzGallStereographicProjectInverse(f *testing.F) { projectInverseFuzz(f, NewGallStereographic()) } -func FuzzMiller(f *testing.F) { +func FuzzMillerProjectInverse(f *testing.F) { projectInverseFuzz(f, NewMiller()) } -func FuzzCentral(f *testing.F) { +func FuzzCentralProjectInverse(f *testing.F) { projectInverseFuzz(f, NewCentral()) } -func FuzzSinusoidal(f *testing.F) { +func FuzzSinusoidalProjectInverse(f *testing.F) { projectInverseFuzz(f, NewSinusoidal()) } -//func FuzzHEALPixStandard(f *testing.F) { +//func FuzzHEALPixStandardProjectInverse(f *testing.F) { // projectInverseFuzz(f, NewHEALPixStandard()) //} -//func FuzzMollweide(f *testing.F) { +//func FuzzMollweideProjectInverse(f *testing.F) { // projectInverseFuzz(f, NewMollweide()) //} -//func FuzzStereographic(f *testing.F) { +//func FuzzStereographicProjectInverse(f *testing.F) { // projectInverseFuzz(f, NewStereographic()) //} -//func FuzzPolar(f *testing.F) { +//func FuzzPolarProjectInverse(f *testing.F) { // projectInverseFuzz(f, NewPolar()) //} -//func FuzzLambertAzimuthal(f *testing.F) { +//func FuzzLambertAzimuthalProjectInverse(f *testing.F) { // projectInverseFuzz(f, NewLambertAzimuthal()) //} -//func FuzzTransverseMercator(f *testing.F) { +//func FuzzTransverseMercatorProjectInverse(f *testing.F) { // projectInverseFuzz(f, NewObliqueProjection(NewMercator(), 0, math.Pi/2, -math.Pi/2)) //} diff --git a/pseudocylindrical.go b/pseudocylindrical.go index 4008ea1..b0d620d 100644 --- a/pseudocylindrical.go +++ b/pseudocylindrical.go @@ -101,3 +101,31 @@ func (h Homolosine) Inverse(x float64, y float64) (lat float64, lon float64) { func (h Homolosine) PlanarBounds() Bounds { return NewRectangleBounds(2*math.Pi, math.Pi) } + +// An equal-area pseudocylindrical projection, in which the polar lines are half the size of the equator. +// https://en.wikipedia.org/wiki/Eckert_IV_projection +type EckertIV struct{} + +func NewEckertIV() EckertIV { + return EckertIV{} +} + +func (e EckertIV) Project(latitude float64, longitude float64) (float64, float64) { + theta := newtonsMethod(latitude, + func(t float64) float64 { return t + math.Sin(2*t)/2 + 2*math.Sin(t) }, + func(t float64) float64 { return 1 + math.Cos(2*t) + 2*math.Cos(t) }, + 1e-4, 1e-15, 125) + return longitude / math.Pi * (1 + math.Cos(theta)), math.Sin(theta) +} + +func (e EckertIV) Inverse(x float64, y float64) (float64, float64) { + theta := math.Asin(y) + latNumer := theta + math.Sin(2*theta)/2 + 2*math.Sin(theta) + lat := math.Asin(latNumer / (2 + math.Pi/2)) + lon := x / (1 + math.Cos(theta)) * math.Pi + return lat, lon +} + +func (e EckertIV) PlanarBounds() Bounds { + return NewRectangleBounds(4, 2) +}