wip

point.go

@@ -38,6 +38,12 @@ func (p Point) X() float64 { return p[0] }
 // Y is the y coordinate of a point in the projection
 func (p Point) Y() float64 { return p[1] }
 
+// Lon is the lon coordinate of a point in the projection
+func (p Point) Lon() float64 { return p[0] }
+
+// Lat is the lat coordinate of a point in the projection
+func (p Point) Lat() float64 { return p[1] }
+
 // MaxX is the same as X
 func (p Point) MaxX() float64 { return p[0] }
 

slippy/maths.go

@@ -0,0 +1,239 @@
+package slippy
+
+import (
+	"fmt"
+	"github.com/go-spatial/geom"
+	"github.com/go-spatial/proj"
+	"math"
+)
+
+/*
+ *
+ * This file should contain the basic math function for converting
+ * between coordinates that are internal to the system.
+ *
+ * Much of the math here is derived from two sources:
+ *  ref: https://maplibre.org/maplibre-native/docs/book/design/coordinate-system.html#11
+ *  ref: https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#ECMAScript_(JavaScript/ActionScript,_etc.)
+ */
+
+const (
+	// DefaultTileSize is the tile size used if the given tile size is 0.
+	DefaultTileSize = 256
+	// Lat4326Max is the maximum degree for latitude on an SRID 4326 map
+	Lat4326Max = 85.05112
+	// Lon4326Max is the maximum degree for longitude on an SRID 4326 map
+	Lon4326Max = 180
+)
+
+type Zoom uint
+
+func (z Zoom) N() float64 { return math.Exp2(float64(z)) }
+
+// Tile2 describes a slippy tile.
+type Tile2 struct {
+	// zoom
+	Z Zoom
+	// column
+	X int
+	// row
+	Y int
+}
+
+// Degree2Radians converts degrees to radians
+func Degree2Radians(degree float64) float64 {
+	return degree * math.Pi / 180
+}
+
+// Radians2Degree converts radians to degrees
+func Radians2Degree(radians float64) float64 {
+	return radians * 180 / math.Pi
+}
+
+func convert3857DegTo4326Deg(z Zoom, lon, lat float64) (Zoom, float64, float64) {
+	latRad := Degree2Radians(lat)
+	nLat := math.Atan(math.Sinh(latRad))
+	nLat = nLat * 180 / math.Pi
+
+	fmt.Printf("conversion: 3857(%v %v) -> (%v %v)\n", lon, lat, lon, nLat)
+	return z, lon, nLat
+}
+
+func convert4326DegTo3857Deg(z Zoom, lon, lat float64) (Zoom, float64, float64) {
+
+	latRad := Degree2Radians(lat)
+	invLat := 1.0 / math.Cos(latRad)
+	tanLat := math.Tan(latRad)
+	sumLat := tanLat + invLat
+	nLat := math.Log(sumLat)
+	dLat := nLat * 180 / math.Pi
+
+	fmt.Printf("conversion: 4326(%v %v) -> (%v %v)\n", lon, lat, lon, dLat)
+	return z, lon, dLat
+}
+
+// lat2Num will return the Y coordinate for the tile at the given Z.
+//
+//	Lat is assumed to be in degrees in SRID 3857 coordinates
+//	If tileSize == 0 then we will use a tileSize of DefaultTileSize
+func lat2Num(tileSize uint32, z Zoom, lat float64) (y int) {
+	//fmt.Printf("new lat2Num: tileSize=%v, z=%v, lat=%v\n", tileSize, z, lat)
+	//return degLat2y(z, lat)
+	fmt.Printf("old lat2Num: tileSize=%v, z=%v\n", tileSize, z)
+	if tileSize == 0 {
+		tileSize = DefaultTileSize
+	}
+	tileY := lat2Px(tileSize, z, lat)
+	tileY = tileY / float64(tileSize)
+	// Truncate to get the tile
+	return int(tileY)
+
+}
+
+// lat2Px will return the pixel coordinate for the lat. This can return
+// a pixel that is outside the extents of the map, this just means
+// the drawing is happening in the buffered area usually done for stitching
+// purposes.
+func lat2Px(tileSize uint32, z Zoom, lat float64) (yPx float64) {
+	if tileSize == 0 {
+		tileSize = DefaultTileSize
+	}
+	n := uint32(1 << z)
+	worldSize := tileSize * n
+
+	// Convert the Degree to radians as most of the math functions work in radians
+	radLat := Degree2Radians(45 + lat/2)
+	// normalize lat
+	latNormalized := math.Log(math.Tan(radLat))
+
+	// compute the pixel value for y:
+	yPxRaw := (180 - Radians2Degree(latNormalized)) / 360
+	yPx = yPxRaw * float64(worldSize)
+	return yPx
+}
+
+// lon2Num will return the Y coordinate for the tile at the given Z.
+//
+//	Lat is assumed to be in degrees in SRID 3857 coordinates
+//	If tileSize == 0 then we will use a tileSize of DefaultTileSize
+func lon2Num(tileSize uint32, z Zoom, lon float64) (x int) {
+	/*
+		if tileSize == 0 {
+			tileSize = DefaultTileSize
+		}
+		tileX := lon2Px(tileSize, z, lon)
+		tileX = tileX / float64(tileSize)
+		// Truncate to get the tile
+		return int(tileX)
+	*/
+	return degLon2x(z, lon)
+
+}
+
+// lonPx will return the pixel coordinate for the lon. This can return
+// a pixels that is outside the extents of the map, this just means
+// the drawing is happening in the buffered area usually done for stitching
+// purposes.
+func lon2Px(tileSize uint32, z uint32, lon float64) (xPx float64) {
+	if tileSize == 0 {
+		tileSize = DefaultTileSize
+	}
+	worldSize := tileSize * (1 << z)
+	lonNormalized := 180 + lon
+	// compute the pixel value for y:
+	xPxRaw := lonNormalized / 360
+	xPx = xPxRaw * float64(worldSize)
+	return xPx
+}
+
+// from ref: https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#ECMAScript_(JavaScript/ActionScript,_etc.)
+func deg2y(z Zoom, lat float64) float64 {
+	n := z.N()
+	radLat := Degree2Radians(lat)
+	lat3857 := math.Log(math.Tan(radLat) + 1.0/math.Cos(radLat))
+	y := lat3857 / math.Pi
+	y = 1.0 - y
+	y = y / 2.0
+	y = y * n
+	y = math.Floor(y)
+	return y
+}
+func degLat2y(z Zoom, lat float64) (y int) {
+	return int(deg2y(z, lat))
+}
+func deg2x(z Zoom, lon float64) float64 {
+	n := math.Exp2(float64(z))
+	return math.Floor((lon + 180.0) / 360.0 * n)
+}
+func degLon2x(z Zoom, lon float64) (x int) {
+	n := z.N()
+	x = int(deg2x(z, lon))
+	if float64(x) >= n {
+		x = int(n - 1)
+	}
+	return x
+}
+func deg2num(z Zoom, lon, lat float64) (x int, y int) {
+	return degLon2x(z, lon), degLat2y(z, lat)
+}
+
+func x2deg(z Zoom, x int) float64 {
+	n := z.N()
+	long := float64(x) / n
+	long = long * 360.0
+	long = long - 180.0
+	return long
+}
+
+func y2deg(z Zoom, y int) float64 {
+	n := math.Pi - 2.0*math.Pi*float64(y)/z.N()
+	lat := 180.0 / math.Pi * math.Atan(0.5*(math.Exp(n)-math.Exp(-n)))
+	return lat
+}
+
+func num2deg(z Zoom, x, y int) (lat float64, long float64) {
+	lat = y2deg(z, y)
+	long = x2deg(z, x)
+	return lat, long
+}
+
+type Grid4326 struct {
+	// TileSize if 0 will default to 256
+	TileSize uint32
+}
+
+func (Grid4326) SRID() proj.EPSGCode { return proj.EPSG4326 }
+func (Grid4326) Size(z Zoom) (Tile2, bool) {
+	n := int(math.Exp2(float64(z)))
+	return Tile2{
+		Z: z,
+		X: n,
+		Y: n,
+	}, true
+}
+
+func (g Grid4326) tileSize() uint32 {
+	if g.TileSize == 0 {
+		return uint32(256)
+	}
+	return g.TileSize
+}
+
+// FromNative will convert a pt in 3857 coordinates and a zoom to a Tile coordinate
+func (g Grid4326) FromNative(z Zoom, pt geom.Point) (tile Tile2, err error) {
+	y := lat2Num(g.tileSize(), z, pt.Lat())
+	x := lon2Num(g.tileSize(), z, pt.Lon())
+	return Tile2{
+		Z: z,
+		X: x,
+		Y: y,
+	}, nil
+}
+func PtFromLatLon(lat, lon float64) geom.Point {
+	return geom.Point{lon, lat}
+}
+func (g Grid4326) ToNative(tile Tile2) (pt geom.Point, err error) {
+	lat := y2deg(tile.Z, tile.Y)
+	lon := x2deg(tile.Z, tile.X)
+	return PtFromLatLon(lat, lon), nil
+}