From 549fbd81de023b3bb389a4198f87737780b8f50b Mon Sep 17 00:00:00 2001 From: Robert Kleffner Date: Mon, 22 Jan 2024 15:42:04 -0700 Subject: [PATCH] Cassini projection, some great new tests (#8) --- README.md | 5 ++-- bounded_test.go | 12 ++++++--- cylindrical.go | 28 ++++++++++++++++++- healpix.go | 2 +- invertability_test.go | 29 ++++++++++++-------- oblique.go | 6 +++-- oblique_test.go | 63 +++++++++++++++++++++++++++++++++++++++++++ project_test.go | 27 +++++++++++++++++++ 8 files changed, 151 insertions(+), 21 deletions(-) create mode 100644 oblique_test.go diff --git a/README.md b/README.md index a76c967..9519994 100644 --- a/README.md +++ b/README.md @@ -64,7 +64,7 @@ Determine how representative of reality a projection is at a point on the sphere A list of the predefined projections supported by the package. -**Invertability**: all projections in this library have a well defined inverse function. However, only some of them reasonably satisfy the invertability property `x = f_1(f(x))` due to floating-point error in computations, edge-cases resulting in infinities or NaNs, or ambiguity at some type of critical point (commonly the poles or prime meridian). The table below checks of *invertible* for all projections which satisfy the `x = f_1(f(x))` property for all valid spherical locations to a reasonably fine degree of precision, referred to here as **everywhere floating-point invertible**. Oblique transforms of floating-point invertible standard projections do not necessarily share that property. +**Invertability**: all projections in this library have a well defined inverse function. However, only some of them reasonably satisfy the invertability property `x = f_1(f(x))` due to floating-point error in computations, edge-cases resulting in infinities or NaNs, or ambiguity at some type of critical point (commonly the poles or prime meridian). The table below checks of *invertible* for all projections which satisfy the `x = f_1(f(x))` property for all valid spherical locations **except the poles** to a reasonably fine degree of precision, referred to here as **everywhere floating-point invertible**. Oblique transforms of floating-point invertible standard projections do not necessarily share that property. Efforts are ongoing to improve the coverage of this property to more projections where possible. @@ -81,7 +81,7 @@ Efforts are ongoing to improve the coverage of this property to more projections |Miller|:white_check_mark:| |Central|:white_check_mark:| |Sinusoidal|:white_check_mark:| -|HEALPix| | +|HEALPix|:white_check_mark:| |Mollweide| | |Homolosine| | |Eckert IV| | @@ -92,6 +92,7 @@ Efforts are ongoing to improve the coverage of this property to more projections |Orthographic| | |Robinson|:white_check_mark:| |Natural Earth|:white_check_mark:| +|Cassini|:white_check_mark:| ## Credits diff --git a/bounded_test.go b/bounded_test.go index 924cd27..fb5d1b7 100644 --- a/bounded_test.go +++ b/bounded_test.go @@ -81,10 +81,6 @@ 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()) //} @@ -105,6 +101,14 @@ func FuzzEqualEarthProjectBounded(f *testing.F) { projectionBoundedFuzz(f, NewEqualEarth()) } +func FuzzCassiniProjectBounded(f *testing.F) { + projectionBoundedFuzz(f, NewCassini()) +} + +func FuzzTransversePlateCarreeProjectBounded(f *testing.F) { + projectionBoundedFuzz(f, NewObliqueProjection(NewPlateCarree(), 0, math.Pi/2, -math.Pi/2)) +} + func projectionBoundedFuzz(f *testing.F, proj Projection) { f.Add(0.0, 0.0) f.Add(0.0, math.Pi) diff --git a/cylindrical.go b/cylindrical.go index d3625f4..52c761a 100644 --- a/cylindrical.go +++ b/cylindrical.go @@ -1,6 +1,8 @@ package flatsphere -import "math" +import ( + "math" +) // An early standard cylindrical projection useful for navigation. // https://en.wikipedia.org/wiki/Mercator_projection @@ -201,3 +203,27 @@ func (c Central) Inverse(x float64, y float64) (lat float64, lon float64) { func (c Central) PlanarBounds() Bounds { return NewRectangleBounds(2*math.Pi, 2*math.Pi) } + +// A transverse version of the Plate–Carée projection, implemented directly for efficiency. +// https://en.wikipedia.org/wiki/Cassini_projection +type Cassini struct{} + +func NewCassini() Cassini { + return Cassini{} +} + +func (c Cassini) Project(lat float64, lon float64) (float64, float64) { + x := math.Asin(math.Cos(lat) * math.Sin(lon)) + y := math.Atan2(math.Sin(lat), math.Cos(lat)*math.Cos(lon)) + return x, y +} + +func (c Cassini) Inverse(x float64, y float64) (float64, float64) { + lat := math.Asin(math.Cos(x) * math.Sin(y)) + lon := math.Atan2(math.Sin(x), math.Cos(x)*math.Cos(y)) + return lat, lon +} + +func (c Cassini) PlanarBounds() Bounds { + return NewRectangleBounds(math.Pi, 2*math.Pi) +} diff --git a/healpix.go b/healpix.go index 73f6b0e..0e4391f 100644 --- a/healpix.go +++ b/healpix.go @@ -48,7 +48,7 @@ func (h HEALPixStandard) Inverse(x float64, y float64) (float64, float64) { lon := facetX + (x-facetX)/sigma return lat, lon } else { - return math.Copysign(math.Pi/2, y), -math.Pi + return math.Copysign(math.Pi/2, y), x } } diff --git a/invertability_test.go b/invertability_test.go index 12421a9..651d818 100644 --- a/invertability_test.go +++ b/invertability_test.go @@ -53,9 +53,9 @@ func FuzzSinusoidalProjectInverse(f *testing.F) { projectInverseFuzz(f, NewSinusoidal()) } -//func FuzzHEALPixStandardProjectInverse(f *testing.F) { -// projectInverseFuzz(f, NewHEALPixStandard()) -//} +func FuzzHEALPixStandardProjectInverse(f *testing.F) { + projectInverseFuzz(f, NewHEALPixStandard()) +} //func FuzzMollweideProjectInverse(f *testing.F) { // projectInverseFuzz(f, NewMollweide()) @@ -97,6 +97,10 @@ func FuzzNaturalEarthProjectInverse(f *testing.F) { // projectInverseFuzz(f, NewEqualEarth()) //} +func FuzzCassiniProjectInverse(f *testing.F) { + projectInverseFuzz(f, NewCassini()) +} + func withinTolerance(n1, n2, tolerance float64) bool { if n1 == n2 { return true @@ -114,16 +118,19 @@ func projectInverseFuzz(f *testing.F, proj Projection) { f.Add(math.Pi/2, math.Pi) f.Add(66.0, 0.0) f.Fuzz(func(t *testing.T, lat float64, lon float64) { - if math.Abs(lat) > math.Pi/2 { - lat = math.Mod(lat, math.Pi/2) - } - if math.Abs(lon) > math.Pi { - lon = math.Mod(lon, math.Pi) - } + lat = math.Mod(lat, math.Pi/2) + lon = math.Mod(lon, math.Pi) x, y := proj.Project(lat, lon) rlat, rlon := proj.Inverse(x, y) - if !withinTolerance(lat, rlat, 0.00001) || !withinTolerance(lon, rlon, 0.00001) { - t.Errorf("expected %e,%e, got %e,%e", lat, lon, rlat, rlon) + + if withinTolerance(lat, math.Pi/2, 0.0000001) || withinTolerance(lat, -math.Pi/2, 0.0000001) { + if !withinTolerance(lat, rlat, 0.000001) { + t.Errorf("expected lat %e, but got %e", lat, rlat) + } + } else { + if !withinTolerance(lat, rlat, 0.00001) || !withinTolerance(lon, rlon, 0.00001) { + t.Errorf("expected %e,%e, got %e,%e", lat, lon, rlat, rlon) + } } }) } diff --git a/oblique.go b/oblique.go index f2c1118..c626a24 100644 --- a/oblique.go +++ b/oblique.go @@ -114,8 +114,10 @@ func (o ObliqueProjection) TransformToOblique(latitude float64, longitude float6 newLat := math.Asin(preAsin) var newLon float64 inner := math.Sin(latitude)/o.cosPoleLat/math.Cos(newLat) - math.Tan(o.poleLat)*math.Tan(newLat) - if o.poleLat == -math.Pi/2 { - newLon = o.poleLon + rotateLon + if o.poleLat == math.Pi/2 { + newLon = rotateLon + o.poleLon + } else if o.poleLat == -math.Pi/2 { + newLon = -rotateLon + o.poleLon + math.Pi } else if math.Abs(inner) > 1 { if (rotateLon == 0 && latitude < -o.poleLat) || (rotateLon != 0 && latitude < o.poleLat) { newLon = o.poleLon + math.Pi diff --git a/oblique_test.go b/oblique_test.go new file mode 100644 index 0000000..1d04fb2 --- /dev/null +++ b/oblique_test.go @@ -0,0 +1,63 @@ +package flatsphere + +import ( + "math" + "testing" +) + +func FuzzObliqueNoOp(f *testing.F) { + obliqEq := NewObliqueProjection(NewPlateCarree(), math.Pi/2, 0, 0) + cassini := NewPlateCarree() + + f.Add(0.0, 0.0) + f.Add(math.Pi/4, math.Pi/4) + f.Add(math.Pi/2, 0.0) + f.Add(0.0, math.Pi/2) + f.Fuzz(func(t *testing.T, lat float64, lon float64) { + if math.Abs(lat) > math.Pi/2 { + lat = math.Mod(lat, math.Pi/2) + } + if math.Abs(lon) > math.Pi { + lon = math.Mod(lon, math.Pi) + } + + xo, yo := obliqEq.Project(lat, lon) + xc, yc := cassini.Project(lat, lon) + + if !withinTolerance(xo, xc, 0.000001) || !withinTolerance(yo, yc, 0.000001) { + t.Errorf("expected %e,%e, but got %e,%e", xc, yc, xo, yo) + } + }) +} + +func rotatePoint(x, y float64, rad float64) (float64, float64) { + rotX := x*math.Cos(rad) - y*math.Sin(rad) + rotY := x*math.Sin(rad) + y*math.Cos(rad) + return rotX, rotY +} + +func FuzzObliquePlateCarreeVsCassiniProject(f *testing.F) { + obliqEq := NewObliqueProjection(NewPlateCarree(), 0, math.Pi/2, -math.Pi/2) + cassini := NewCassini() + + f.Add(0.0, 0.0) + f.Add(math.Pi/4, math.Pi/4) + f.Add(math.Pi/2, 0.0) + f.Add(0.0, math.Pi/2) + f.Fuzz(func(t *testing.T, lat float64, lon float64) { + if math.Abs(lat) > math.Pi/2 { + lat = math.Mod(lat, math.Pi/2) + } + if math.Abs(lon) > math.Pi { + lon = math.Mod(lon, math.Pi) + } + + xo, yo := obliqEq.Project(lat, lon) + xc, yc := cassini.Project(lat, lon) + xr, yr := rotatePoint(xc, yc, math.Pi/2) + + if !withinTolerance(xo, xr, 0.000001) || !withinTolerance(yo, yr, 0.000001) { + t.Errorf("expected %e,%e, but got %e,%e", xr, yr, xo, yo) + } + }) +} diff --git a/project_test.go b/project_test.go index 88f19f4..c4f3768 100644 --- a/project_test.go +++ b/project_test.go @@ -60,6 +60,33 @@ func TestMercatorInverseSanity(t *testing.T) { }) } +func TestCassiniProjectSanity(t *testing.T) { + checkProject(t, "cassini", NewCassini(), []projectTestCase{ + {0, 0, 0, 0}, + {0, math.Pi, 0, math.Pi}, + {0, math.Pi / 2, math.Pi / 2, 0}, + {math.Pi / 4, 0, 0, math.Pi / 4}, + {-math.Pi / 4, 0, 0, -math.Pi / 4}, + {math.Pi / 2, 0, 0, math.Pi / 2}, + {-math.Pi / 2, 0, 0, -math.Pi / 2}, + {math.Pi / 2, -math.Pi, 0, math.Pi / 2}, + {-math.Pi / 2, math.Pi, 0, -math.Pi / 2}, + }) +} + +func TestCassiniInverseSanity(t *testing.T) { + checkInverse(t, "invCassini", NewCassini(), []inverseTestCase{ + {0, 0, 0, 0}, + {0, math.Pi, 0, math.Pi}, + {math.Pi / 2, 0, 0, math.Pi / 2}, + {0, math.Pi / 4, math.Pi / 4, 0}, + {0, -math.Pi / 4, -math.Pi / 4, 0}, + {0, math.Pi / 2, math.Pi / 2, 0}, + {0, -math.Pi / 2, -math.Pi / 2, 0}, + {math.Pi / 2, math.Pi, 0, 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