Tabular projections (Robinson & Natural Earth) (#6)

* Tabular projections (Robinson & Natural Earth)

* Upgrade automated build Go version

.github/workflows/go.yml

@@ -19,7 +19,7 @@ jobs:
     - name: Set up Go
       uses: actions/setup-go@v4
       with:
-        go-version: '1.19'
+        go-version: '1.21'
 
     - name: Build
       run: go build -v ./...

README.md

@@ -7,7 +7,7 @@ A [Go](https://go.dev/) library for converting between spherical and planar coor
 
 ## Prerequisites
 
-- **[Go](https://go.dev/)**: any one of the **three latest major** [releases](https://go.dev/doc/devel/release).
+- **[Go](https://go.dev/)**: requires [Go 1.21](https://go.dev/doc/devel/release#go1.21.0) or higher.
 
 ## Install
 
@@ -83,11 +83,15 @@ Efforts are ongoing to improve the coverage of this property to more projections
 |Sinusoidal|:white_check_mark:|
 |HEALPix| |
 |Mollweide| |
+|Homolosine| |
+|Eckert IV| |
 |Stereographic| |
 |Polar| |
 |Lambert azimuthal| |
 |Gnomonic| |
 |Orthographic| |
+|Robinson|:white_check_mark:|
+|Natural Earth|:white_check_mark:|
 
 ## Credits
 

analysis.go

@@ -30,3 +30,30 @@ func newtonsMethod(
 	}
 	return math.NaN()
 }
+
+// Given a table of x,y values representing points on a function, Aitken interpolation approximates
+// the y value of the described function at a given x value, estimating using the points in the table
+// between xStart and xEnd.
+func aitkenInterpolation(xs []float64, ys []float64, xStart int, xEnd int, atX float64) float64 {
+	samples := xEnd - xStart
+	approximations := make([][]float64, samples)
+
+	approximations[0] = ys[xStart:xEnd] // fill in the zero row
+
+	// iterate recursively, forming a triangular matrix
+	for i := 1; i < samples; i++ {
+		approximations[i] = make([]float64, samples)
+		for j := i; j < samples; j++ {
+			determ := determinant(approximations[i-1][i-1], approximations[i-1][j], xs[xStart+i-1]-atX, xs[xStart+j]-atX)
+			denom := xs[xStart+j] - xs[xStart+i-1]
+			approximations[i][j] = determ / denom
+		}
+	}
+
+	// the last value is the final interpolated value
+	return approximations[samples-1][samples-1]
+}
+
+func determinant(a, b, c, d float64) float64 {
+	return a*d - b*c
+}

analysis_test.go

@@ -5,6 +5,8 @@ import (
 	"testing"
 )
 
+// Verify that our approximation method applied to the square root function compares
+// to the built-in square root function nicely.
 func FuzzNewtonsMethodSqrt(f *testing.F) {
 	f.Add(612.0)
 	f.Add(4.0)
@@ -24,3 +26,21 @@ func FuzzNewtonsMethodSqrt(f *testing.F) {
 		}
 	})
 }
+
+func FuzzAitkenInterpolation(f *testing.F) {
+	xs := []float64{0.0, 1.0, 2.0, 3.0}
+	ys := []float64{5.0, 2.0, -5.0, -10.0}
+	poly := func(x float64) float64 { return 5.0 + x - 5.0*x*x + x*x*x }
+	f.Add(0.0)
+	f.Add(1.0)
+	f.Add(2.0)
+	f.Add(-1.0)
+	f.Fuzz(func(t *testing.T, x float64) {
+		x = math.Mod(x, xs[len(xs)-1])
+		approx := aitkenInterpolation(xs, ys, 0, len(xs), x)
+		polyRes := poly(x)
+		if !withinTolerance(approx, polyRes, 0.0001) {
+			t.Errorf("expected %f, got %f", polyRes, approx)
+		}
+	})
+}

bounded_test.go

@@ -93,6 +93,14 @@ func FuzzOrthographicProjectBounded(f *testing.F) {
 	projectionBoundedFuzz(f, NewOrthographic())
 }
 
+func FuzzRobinsonProjectBounded(f *testing.F) {
+	projectionBoundedFuzz(f, NewRobinsonProjection())
+}
+
+func FuzzNaturalEarthProjectBounded(f *testing.F) {
+	projectionBoundedFuzz(f, NewNaturalEarthProjection())
+}
+
 func projectionBoundedFuzz(f *testing.F, proj Projection) {
 	f.Add(0.0, 0.0)
 	f.Add(0.0, math.Pi)

go.mod

@@ -1,3 +1,3 @@
 module github.com/owlpinetech/flatsphere
 
-go 1.19
+go 1.21

invertability_test.go

@@ -61,6 +61,14 @@ func FuzzSinusoidalProjectInverse(f *testing.F) {
 //	projectInverseFuzz(f, NewMollweide())
 //}
 
+//func FuzzHomolosineProjectInverse(f *testing.F) {
+//	projectInverseFuzz(f, NewHomolosine())
+//}
+
+//func FuzzEckertIVProjectInverse(f *testing.F) {
+//	projectInverseFuzz(f, NewEckertIV())
+//}
+
 //func FuzzStereographicProjectInverse(f *testing.F) {
 //	projectInverseFuzz(f, NewStereographic())
 //}
@@ -77,6 +85,14 @@ func FuzzSinusoidalProjectInverse(f *testing.F) {
 //	projectInverseFuzz(f, NewObliqueProjection(NewMercator(), 0, math.Pi/2, -math.Pi/2))
 //}
 
+func FuzzRobinsonProjectInverse(f *testing.F) {
+	projectInverseFuzz(f, NewRobinsonProjection())
+}
+
+func FuzzNaturalEarthProjectInverse(f *testing.F) {
+	projectInverseFuzz(f, NewNaturalEarthProjection())
+}
+
 func withinTolerance(n1, n2, tolerance float64) bool {
 	if n1 == n2 {
 		return true

tabular.go

@@ -0,0 +1,76 @@
+package flatsphere
+
+import (
+	"math"
+	"slices"
+)
+
+// A class of pseudocylindrical projections defined by a table of ratio values at particular latitutdes.
+// Ratio values between the latitudes with defined entries are computed via interpolation during projection/inverse.
+// The polynomial order parameter is use to control the interpolation expensiveness/accuracy, and should be a positive
+// even number. The y scale parameter is used to control the scaling of latitude projection relative to width, allowing
+// the parallel distance ratios to be input in normalized (-1, 1) range (i.e. as an actual ratio).
+type TabularProjection struct {
+	halfPolynomialOrder int
+	yScale              float64
+	latitudes           []float64
+	parallelLengthRatio []float64
+	parallelDistRatio   []float64
+}
+
+// Create a new pseudocylindrical projection from a table of values. The tables should be sorted by the latitude
+// entry, such that the least latitude, and its corresponding ratio entries, are the first row in the table.
+// Polynomial order should be a positive even integer, bigger = more accurate but more compute. yScale should be
+// a positive float in the range (0,1).
+func NewTabularProjection(
+	latitudes []float64,
+	parallelLengthRatios []float64,
+	parallelDistanceRatios []float64,
+	polynomialOrder int,
+	yScale float64,
+) TabularProjection {
+	return TabularProjection{polynomialOrder / 2, yScale, latitudes, parallelLengthRatios, parallelDistanceRatios}
+}
+
+var (
+	robinsonNaturalEarthLatitudes []float64 = []float64{-90, -85, -80, -75, -70, -65, -60, -55, -50, -45, -40, -35, -30, -25, -20, -15, -10, -5, 0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90}
+	robinsonLengthRatios          []float64 = []float64{0.5322, 0.5722, 0.6213, 0.6732, 0.7186, 0.7597, 0.7986, 0.8350, 0.8679, 0.8962, 0.9216, 0.9427, 0.9600, 0.9730, 0.9822, 0.9900, 0.9954, 0.9986, 1.0000, 0.9986, 0.9954, 0.9900, 0.9822, 0.9730, 0.9600, 0.9427, 0.9216, 0.8962, 0.8679, 0.8350, 0.7986, 0.7597, 0.7186, 0.6732, 0.6213, 0.5722, 0.5322}
+	robinsonDistRatios            []float64 = []float64{-1.0000, -0.9761, -0.9394, -0.8936, -0.8435, -0.7903, -0.7346, -0.6769, -0.6176, -0.5571, -0.4958, -0.4340, -0.3720, -0.3100, -0.2480, -0.1860, -0.1240, -0.0620, 0.0000, 0.0620, 0.1240, 0.1860, 0.2480, 0.3100, 0.3720, 0.4340, 0.4958, 0.5571, 0.6176, 0.6769, 0.7346, 0.7903, 0.8435, 0.8936, 0.9394, 0.9761, 1.0000}
+	naturalEarthLengthRatios      []float64 = []float64{0.5630, 0.6270, 0.6754, 0.7160, 0.7525, 0.7874, 0.8196, 0.8492, 0.8763, 0.9006, 0.9222, 0.9409, 0.9570, 0.9703, 0.9811, 0.9894, 0.9953, 0.9988, 1.0000, 0.9988, 0.9953, 0.9894, 0.9811, 0.9703, 0.9570, 0.9409, 0.9222, 0.9006, 0.8763, 0.8492, 0.8196, 0.7874, 0.7525, 0.7160, 0.6754, 0.6270, 0.5630}
+	naturalEarthDistRatios        []float64 = []float64{-1.0000, -0.9761, -0.9394, -0.8936, -0.8435, -0.7903, -0.7346, -0.6769, -0.6176, -0.5571, -0.4958, -0.4340, -0.3720, -0.3100, -0.2480, -0.1860, -0.1240, -0.0620, 0.0000, 0.0620, 0.1240, 0.1860, 0.2480, 0.3100, 0.3720, 0.4340, 0.4958, 0.5571, 0.6176, 0.6769, 0.7346, 0.7903, 0.8435, 0.8936, 0.9394, 0.9761, 1.0000}
+)
+
+// 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 {
+	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 {
+	return NewTabularProjection(robinsonNaturalEarthLatitudes, naturalEarthLengthRatios, naturalEarthDistRatios, 4, 0.520)
+}
+
+func (t TabularProjection) Project(lat float64, lon float64) (float64, float64) {
+	latDegree := lat * math.Pi / 180
+	xInterp := t.interpolate(latDegree, t.latitudes, t.parallelLengthRatio)
+	yInterp := t.interpolate(latDegree, t.latitudes, t.parallelDistRatio)
+	return lon / math.Pi * xInterp, t.yScale * yInterp
+}
+
+func (t TabularProjection) Inverse(x float64, y float64) (float64, float64) {
+	yNorm := y / t.yScale
+	latInterp := t.interpolate(yNorm, t.parallelDistRatio, t.latitudes)
+	lonInterp := t.interpolate(yNorm, t.parallelDistRatio, t.parallelLengthRatio)
+	return latInterp * 180 / math.Pi, math.Pi * x / lonInterp
+}
+
+func (t TabularProjection) interpolate(at float64, xs []float64, ys []float64) float64 {
+	ind, _ := slices.BinarySearch(t.latitudes, at)
+	return aitkenInterpolation(xs, ys, max(ind-t.halfPolynomialOrder, 0), min(ind+t.halfPolynomialOrder, len(xs)), at)
+}
+
+func (t TabularProjection) PlanarBounds() Bounds {
+	return NewRectangleBounds(2, t.yScale*2)
+}