Skip to content

Commit

Permalink
A few elliptical/circular projections (#9)
Browse files Browse the repository at this point in the history
* Aitoff projection

* Add an example for Go pkg documentation

* Simple distortion calculation tests added

* Add Hammer projection and tests

* A couple fixes

* A couple more fixes for now

* Another readme update

* Lagrange projection, fixing a few planar bounds that diverge
  • Loading branch information
robertkleffner authored Feb 10, 2024
1 parent 549fbd8 commit 586f5c0
Show file tree
Hide file tree
Showing 10 changed files with 243 additions and 31 deletions.
3 changes: 3 additions & 0 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -93,6 +93,9 @@ Efforts are ongoing to improve the coverage of this property to more projections
|Robinson|:white_check_mark:|
|Natural Earth|:white_check_mark:|
|Cassini|:white_check_mark:|
|Aitoff| |
|Hammer| |
|Lagrange|:white_check_mark:|

## Credits

Expand Down
6 changes: 1 addition & 5 deletions azimuthal.go
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
package flatsphere

import (
"fmt"
"math"
)

Expand All @@ -24,7 +23,7 @@ func (s Stereographic) Inverse(x float64, y float64) (float64, float64) {
}

func (s Stereographic) PlanarBounds() Bounds {
return NewRectangleBounds(4, 4)
return NewBounds(math.Inf(-1), math.Inf(-1), math.Inf(1), math.Inf(1))
}

// An ancient equidistant azimuthal projection.
Expand Down Expand Up @@ -64,9 +63,6 @@ func (p LambertAzimuthal) Project(latitude float64, longitude float64) (float64,
func (p LambertAzimuthal) Inverse(x float64, y float64) (float64, float64) {
r := math.Hypot(x, y)
rLat, rLon := math.Asin(1-2*r*r), math.Atan2(x, -y)
if math.IsNaN(rLat) {
fmt.Printf("nan lat: r = %e, x = %e, y = %e, presin = %e, valid = %v\n", r, x, y, 1-2*r*r, 1-2*r*r < -1)
}
return rLat, rLon
}

Expand Down
20 changes: 17 additions & 3 deletions bounded_test.go
Original file line number Diff line number Diff line change
Expand Up @@ -65,9 +65,9 @@ func FuzzHomolosineProjectBounded(f *testing.F) {
projectionBoundedFuzz(f, NewHomolosine())
}

func FuzzEckertIVProjectBounded(f *testing.F) {
projectionBoundedFuzz(f, NewEckertIV())
}
//func FuzzEckertIVProjectBounded(f *testing.F) {
// projectionBoundedFuzz(f, NewEckertIV())
//}

func FuzzStereographicProjectBounded(f *testing.F) {
projectionBoundedFuzz(f, NewStereographic())
Expand Down Expand Up @@ -109,7 +109,21 @@ func FuzzTransversePlateCarreeProjectBounded(f *testing.F) {
projectionBoundedFuzz(f, NewObliqueProjection(NewPlateCarree(), 0, math.Pi/2, -math.Pi/2))
}

func FuzzAitoffProjectBounded(f *testing.F) {
projectionBoundedFuzz(f, NewAitoff())
}

func FuzzHammerProjectBounded(f *testing.F) {
projectionBoundedFuzz(f, NewHammer())
}

func FuzzLagrangeProjectBounded(f *testing.F) {
projectionBoundedFuzz(f, NewLagrange())
}

func projectionBoundedFuzz(f *testing.F, proj Projection) {
f.Add(109.95574287564276, 17.0)
f.Add(-15.707963267948964, -0.09817477042468103)
f.Add(0.0, 0.0)
f.Add(0.0, math.Pi)
f.Add(math.Pi/2, math.Pi/4)
Expand Down
11 changes: 11 additions & 0 deletions bounds.go
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,17 @@ func (b Bounds) Within(x float64, y float64) bool {
y <= b.YMax
}

// Construct a new bounding area from the minimum and maximum along each of the two axes.
// The center point will be a position halfway between the minimum and maximum.
func NewBounds(xmin, ymin, xmax, ymax float64) Bounds {
return Bounds{
XMin: xmin,
YMin: ymin,
XMax: xmax,
YMax: ymax,
}
}

// Construct a bounding area containing the circle described by the given
// radius, centered on the origin.
func NewCircleBounds(radius float64) Bounds {
Expand Down
6 changes: 3 additions & 3 deletions cylindrical.go
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ func (m Mercator) Inverse(x float64, y float64) (lat float64, lon float64) {
}

func (m Mercator) PlanarBounds() Bounds {
return NewRectangleBounds(2*math.Pi, 2*math.Pi)
return NewBounds(-math.Pi, math.Inf(-1), math.Pi, math.Inf(1))
}

// A special case of the equirectangular projection which allows for easy conversion between
Expand Down Expand Up @@ -160,7 +160,7 @@ func (g GallStereographic) Inverse(x float64, y float64) (lat float64, lon float
}

func (g GallStereographic) PlanarBounds() Bounds {
return NewRectangleBounds(2*math.Pi, 1.5*math.Pi)
return NewBounds(-math.Pi, math.Inf(-1), math.Pi, math.Inf(1))
}

// A compromise cylindrical projection intended to resemble Mercator with less distortion at the poles.
Expand Down Expand Up @@ -201,7 +201,7 @@ 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)
return NewBounds(-math.Pi, math.Inf(-1), math.Pi, math.Inf(1))
}

// A transverse version of the Plate–Carée projection, implemented directly for efficiency.
Expand Down
22 changes: 22 additions & 0 deletions distortion_test.go
Original file line number Diff line number Diff line change
@@ -0,0 +1,22 @@
package flatsphere

import "testing"

func Test_DistortionNull(t *testing.T) {
proj := NewMercator()
eqArea, eqAngular := DistortionAt(proj, 0, 0)

if !withinTolerance(eqArea, 0.0, 0.000001) || !withinTolerance(eqAngular, 0.0, 0.000001) {
t.Errorf("expected 0, 0 distortion at mercator equator, got %f, %f", eqArea, eqAngular)
}

eqArea = AreaDistortionAt(proj, 0, 0)
if !withinTolerance(eqArea, 0.0, 0.000001) {
t.Errorf("expected 0 area distortion at mercator equator, got %f", eqArea)
}

eqAngular = AngularDistortionAt(proj, 0, 0)
if !withinTolerance(eqAngular, 0.0, 0.000001) {
t.Errorf("expected 0 angular distortion at mercator equator, got %f", eqAngular)
}
}
28 changes: 28 additions & 0 deletions example_reproject_test.go
Original file line number Diff line number Diff line change
@@ -0,0 +1,28 @@
package flatsphere

import (
"fmt"
"math"
)

// Example_Reproject demonstrates taking a position in one projected plane and converting it into the analogous position
// in a different projection.
func Example_reproject() {
// determine the original projection of the data
original := NewMercator()
// decide on a new desired projection of the data
target := NewCassini()

// take a point within the PlanarBounds() of the original projection
ox, oy := math.Pi, math.Pi
// get a latitude and longitude on the sphere from the original projected point
lat, lon := original.Inverse(ox, oy)
// get a projected point in the desired projection
cx, cy := target.Project(lat, lon)

fmt.Printf("original x, y: %f, %f\n", ox, oy)
fmt.Printf("target x, y: %f, %f\n", cx, cy)

// Output: original x, y: 3.141593, 3.141593
// target x, y: 0.000000, 1.657170
}
14 changes: 13 additions & 1 deletion invertability_test.go
Original file line number Diff line number Diff line change
Expand Up @@ -101,6 +101,18 @@ func FuzzCassiniProjectInverse(f *testing.F) {
projectInverseFuzz(f, NewCassini())
}

//func FuzzAitoffProjectInverse(f *testing.F) {
// projectInverseFuzz(f, NewAitoff())
//}

//func FuzzHammerProjectInverse(f *testing.F) {
// projectInverseFuzz(f, NewHammer())
//}

func FuzzLagrangeProjectInverse(f *testing.F) {
projectInverseFuzz(f, NewLagrange())
}

func withinTolerance(n1, n2, tolerance float64) bool {
if n1 == n2 {
return true
Expand All @@ -125,7 +137,7 @@ func projectInverseFuzz(f *testing.F, proj Projection) {

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)
t.Errorf("expected lat %e, but got %e from %e, %e", lat, rlat, x, y)
}
} else {
if !withinTolerance(lat, rlat, 0.00001) || !withinTolerance(lon, rlon, 0.00001) {
Expand Down
116 changes: 116 additions & 0 deletions lenticular.go
Original file line number Diff line number Diff line change
@@ -0,0 +1,116 @@
package flatsphere

import (
"math"
)

// A compromise azimuthal projection stretched into an elliptical shape.
// https://en.wikipedia.org/wiki/Aitoff_projection
type Aitoff struct{}

func NewAitoff() Aitoff {
return Aitoff{}
}

func (ai Aitoff) Project(lat float64, lon float64) (float64, float64) {
a := math.Acos(math.Cos(lat) * math.Cos(lon/2))
if a == 0 {
return 0, 0
}
x := 2 * math.Cos(lat) * math.Sin(lon/2) * a / math.Sin(a)
y := math.Sin(lat) * a / math.Sin(a)
return x, y
}

func (ai Aitoff) Inverse(x float64, y float64) (float64, float64) {
interLat, interLon := NewPolar().Inverse(x/2, y)
transLat, transLon := NewObliqueAspect(0, 0, 0).TransformToOblique(interLat, interLon)
return transLat, transLon * 2
}

func (a Aitoff) PlanarBounds() Bounds {
return NewEllipseBounds(math.Pi, math.Pi/2)
}

// An elliptical equal-area projection.
// https://en.wikipedia.org/wiki/Hammer_projection
type Hammer struct{}

func NewHammer() Hammer {
return Hammer{}
}

func (h Hammer) Project(lat float64, lon float64) (float64, float64) {
z := math.Sqrt(1 + math.Cos(lat)*math.Cos(lon/2))
x := 2 * math.Cos(lat) * math.Sin(lon/2) / z
y := math.Sin(lat) / z
return x, y
}

func (h Hammer) Inverse(x float64, y float64) (float64, float64) {
z := math.Sqrt(1 - x*x/8 - y*y/2)
shift := 0.0
if math.Hypot(x/2, y) > 1 {
shift = 2 * math.Pi
if math.Signbit(x) {
shift = -shift
}
}
preAsin := z * y * math.Sqrt2
if preAsin > 1 && preAsin < 1+1e-9 {
preAsin = 1
}
if preAsin < -1 && preAsin > -1-1e-9 {
preAsin = -1
}
lat := math.Asin(preAsin)
lon := 2*math.Atan(math.Sqrt(0.5)*z*x/(2*z*z-1)) + shift
return lat, lon
}

func (h Hammer) PlanarBounds() Bounds {
return NewEllipseBounds(2, 1)
}

// A circular conformal projection of the whole earth.
type Lagrange struct{}

func NewLagrange() Lagrange {
return Lagrange{}
}

func (l Lagrange) Project(lat float64, lon float64) (float64, float64) {
p := (1 + math.Sin(lat)) / (1 - math.Sin(lat))
v := math.Pow(p, 0.25)
c := (v+1/v)/2 + math.Cos(lon/2)
if math.IsInf(c, 0) {
if math.Signbit(lat) {
return 0, -1
} else {
return 0, 1
}
}
x := math.Sin(lon/2) / c
y := (v - 1/v) / (2 * c)
return x, y
}

func (l Lagrange) Inverse(x float64, y float64) (float64, float64) {
r2 := x*x + y*y
th := 2 * y / (1 + r2)
if th == 1 {
if math.Signbit(y) {
return -math.Pi / 2, 0
} else {
return math.Pi / 2, 0
}
}
t := math.Pow((1+th)/(1-th), 2)
lat := math.Asin((t - 1) / (t + 1))
lon := 2 * math.Atan2(2*x, 1-r2)
return lat, lon
}

func (l Lagrange) PlanarBounds() Bounds {
return NewCircleBounds(1)
}
48 changes: 29 additions & 19 deletions oblique.go
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,7 @@ import (
"math"
)

type ObliqueProjection struct {
orig Projection
type ObliqueAspect struct {
poleLat float64
poleLon float64
poleTheta float64
Expand All @@ -14,9 +13,8 @@ type ObliqueProjection struct {
cosPoleLat float64 // cached precompute of cos(poleLat)
}

func NewObliqueProjection(original Projection, poleLat float64, poleLon float64, poleTheta float64) ObliqueProjection {
return ObliqueProjection{
orig: original,
func NewObliqueAspect(poleLat float64, poleLon float64, poleTheta float64) ObliqueAspect {
return ObliqueAspect{
poleLat: poleLat,
poleLon: poleLon,
poleTheta: poleTheta,
Expand All @@ -25,22 +23,10 @@ func NewObliqueProjection(original Projection, poleLat float64, poleLon float64,
}
}

func (o ObliqueProjection) Project(latitude float64, longitude float64) (float64, float64) {
return o.orig.Project(o.TransformFromOblique(latitude, longitude))
}

func (o ObliqueProjection) Inverse(x float64, y float64) (float64, float64) {
return o.TransformToOblique(o.orig.Inverse(x, y))
}

func (o ObliqueProjection) PlanarBounds() Bounds {
return o.orig.PlanarBounds()
}

// Applies the pole shift and rotation of the Oblique projection transform to the given
// input latitude and longitude points, so that the returned latitude/longitude are able
// to be used for the non-transformed 'original' projection.
func (o ObliqueProjection) TransformFromOblique(latitude float64, longitude float64) (float64, float64) {
func (o ObliqueAspect) TransformFromOblique(latitude float64, longitude float64) (float64, float64) {
var newLat, newLon float64

poleRelCos := math.Cos(o.poleLon - longitude)
Expand Down Expand Up @@ -101,7 +87,7 @@ func coerceAngle(angle float64) float64 {
// Given a latitude/longitude in the non-transformed 'original' projection space, applies
// the pole shift and rotation of the Oblique projection so that the returned latitude/longitude
// are in the Oblique projection space.
func (o ObliqueProjection) TransformToOblique(latitude float64, longitude float64) (float64, float64) {
func (o ObliqueAspect) TransformToOblique(latitude float64, longitude float64) (float64, float64) {
rotateLon := longitude + o.poleTheta

preAsin := o.sinPoleLat*math.Sin(latitude) - o.cosPoleLat*math.Cos(latitude)*math.Cos(rotateLon)
Expand Down Expand Up @@ -135,3 +121,27 @@ func (o ObliqueProjection) TransformToOblique(latitude float64, longitude float6
}
return newLat, newLon
}

type ObliqueProjection struct {
orig Projection
ObliqueAspect
}

func NewObliqueProjection(original Projection, poleLat float64, poleLon float64, poleTheta float64) ObliqueProjection {
return ObliqueProjection{
original,
NewObliqueAspect(poleLat, poleLon, poleTheta),
}
}

func (o ObliqueProjection) Project(latitude float64, longitude float64) (float64, float64) {
return o.orig.Project(o.TransformFromOblique(latitude, longitude))
}

func (o ObliqueProjection) Inverse(x float64, y float64) (float64, float64) {
return o.TransformToOblique(o.orig.Inverse(x, y))
}

func (o ObliqueProjection) PlanarBounds() Bounds {
return o.orig.PlanarBounds()
}

0 comments on commit 586f5c0

Please sign in to comment.