Get slippy tiles to work.

bbox.go

@@ -9,6 +9,13 @@ type Extenter interface {
 	Extent() (extent [4]float64)
 }
 
+type PtMinMaxer interface {
+	// Min returns the minimum x and y values
+	Min() Point
+	// Max returns the maximum x and y values
+	Max() Point
+}
+
 // MinMaxer is a wrapper for an Extent that gets min/max of the extent
 type MinMaxer interface {
 	MinX() float64
@@ -51,31 +58,31 @@ func (e *Extent) Edges(cwfn ClockwiseFunc) [][2][2]float64 {
 	}
 }
 
-// MaxX is the larger of the x values.
+// MaxX is the largest of the x values.
 func (e *Extent) MaxX() float64 {
 	if e == nil {
 		return math.MaxFloat64
 	}
 	return e[2]
 }
 
-// MinX  is the smaller of the x values.
+// MinX  is the smallest of the x values.
 func (e *Extent) MinX() float64 {
 	if e == nil {
 		return -math.MaxFloat64
 	}
 	return e[0]
 }
 
-// MaxY is the larger of the y values.
+// MaxY is the largest of the y values.
 func (e *Extent) MaxY() float64 {
 	if e == nil {
 		return math.MaxFloat64
 	}
 	return e[3]
 }
 
-// MinY is the smaller of the y values.
+// MinY is the smallest of the y values.
 func (e *Extent) MinY() float64 {
 	if e == nil {
 		return -math.MaxFloat64
@@ -84,13 +91,13 @@ func (e *Extent) MinY() float64 {
 }
 
 // Min returns the (MinX, MinY) values
-func (e *Extent) Min() [2]float64 {
-	return [2]float64{e[0], e[1]}
+func (e *Extent) Min() Point {
+	return Point{e[0], e[1]}
 }
 
 // Max returns the (MaxX, MaxY) values
-func (e *Extent) Max() [2]float64 {
-	return [2]float64{e[2], e[3]}
+func (e *Extent) Max() Point {
+	return Point{e[2], e[3]}
 }
 
 // XSpan is the distance of the Extent in X or inf
@@ -242,7 +249,7 @@ func NewExtentFromGeometry(g Geometry) (*Extent, error) {
 	return &e, nil
 }
 
-// Contains will return whether the given  extent is inside of the  extent.
+// Contains will return whether the given  extent is inside the  extent.
 func (e *Extent) Contains(ne MinMaxer) bool {
 	// Nil extent contains the world.
 	if e == nil {
@@ -333,6 +340,7 @@ func (e *Extent) Clone() *Extent {
 //	+--------------+----------+      |
 //	               |             B   |
 //	               +-----------------+
+//
 // For example the for the above Box A intersects Box B at the area surround by C.
 //
 // If the Boxes don't intersect does will be false, otherwise ibb will be the intersect.

cmd/main.go

@@ -55,7 +55,7 @@ func readInputWKT(filename string) (geom.Geometry, error) {
 }
 
 type outfile struct {
-	tile   *slippy.Tile
+	tile   slippy.Tile
 	format string
 }
 type outfilefile struct {
@@ -83,7 +83,7 @@ func (off *outfilefile) WriteWKTGeom(geos ...geom.Geometry) *outfilefile {
 	return off
 }
 
-func newOutFile(tile *slippy.Tile, tag string) outfile {
+func newOutFile(tile slippy.Tile, tag string) outfile {
 	path := fmt.Sprintf("%v/%v/%v", tile.Z, tile.X, tile.Y)
 	if tag != "" {
 		path = fmt.Sprintf("%v/%v", path, tag)
@@ -96,6 +96,17 @@ func newOutFile(tile *slippy.Tile, tag string) outfile {
 	}
 }
 
+// MvtTileDim is the number of pixels in a tile
+const MvtTileDim = 4096.0
+
+func PixelToNative(g slippy.TileGridder, z slippy.Zoom) (float64, error) {
+	ext, err := slippy.Extent(g, slippy.Tile{Z: z})
+	if err != nil {
+		return 0, err
+	}
+	return ext.XSpan() / MvtTileDim, nil
+}
+
 func main() {
 	flag.Parse()
 	if len(flag.Args()) < 2 || *help {
@@ -126,37 +137,40 @@ func main() {
 		fmt.Fprintf(os.Stderr, "Unabled to parse y: %v", err)
 		usage()
 	}
-	tile := slippy.NewTile(uint(z), uint(x), uint(y))
+	tile := slippy.Tile{
+		Z: slippy.Zoom(z),
+		X: uint(x),
+		Y: uint(y),
+	}
 	fileTemplate := newOutFile(tile, *tag)
 	geo, err := readInputWKT(flag.Args()[1])
 	if err != nil {
 		fmt.Fprintf(os.Stderr, "Unabled to parse/open `%v` : %v", os.Args, err)
 		usage()
 	}
 	ctx := context.Background()
-	/*
-		plywkt, err := wkt.EncodeString(geo)
-		if err != nil {
-			panic(err)
-		}
-		fmt.Printf("Polygon:\n%v\n", plywkt)
-	*/
 	order := winding.Order{}
-	grid3857, _ := slippy.NewGrid(3857)
+	grid3857 := slippy.NewGrid(3857, 0)
 
 	var clipRegion *geom.Extent
 	{
-		webs := slippy.PixelsToNative(grid3857, tile.Z, uint(*buffer))
+		webs, err := PixelToNative(grid3857, tile.Z)
+		if err != nil {
+			panic(err)
+		}
 		ext, _ := slippy.Extent(grid3857, tile)
 		clipRegion = ext.ExpandBy(webs)
 	}
 
 	if *simplifyGeo {
+		tol, err := PixelToNative(grid3857, tile.Z)
+		if err != nil {
+			panic(err)
+		}
 		simp := simplify.DouglasPeucker{
-			Tolerance: slippy.PixelsToNative(grid3857, tile.Z, 10.0),
+			Tolerance: tol,
 		}
 
-		var err error
 		geo, err = planar.Simplify(ctx, simp, geo)
 		if err != nil {
 			fmt.Fprintf(os.Stderr, "Unabled to simplify geo : %v", err)

go.mod

@@ -1,6 +1,6 @@
 module github.com/go-spatial/geom
 
-go 1.13
+go 1.22
 
 require (
 	github.com/arolek/p v0.0.0-20191103215535-df3c295ed582
@@ -10,6 +10,10 @@ require (
 	github.com/mattn/go-sqlite3 v1.12.0
 	github.com/mattn/goveralls v0.0.3-0.20180319021929-1c14a4061c1c
 	github.com/pborman/uuid v1.2.0
-	github.com/stretchr/testify v1.4.0 // indirect
 	golang.org/x/tools v0.0.0-20191114222411-4191b8cbba09
 )
+
+require (
+	github.com/google/uuid v1.0.0 // indirect
+	github.com/stretchr/testify v1.4.0 // indirect
+)

go.sum

@@ -29,7 +29,6 @@ golang.org/x/text v0.3.0/go.mod h1:NqM8EUOU14njkJ3fqMW+pc6Ldnwhi/IjpwHt7yyuwOQ=
 golang.org/x/tools v0.0.0-20191114222411-4191b8cbba09 h1:f3ZhcxGnJxVhcsQgKPvfAg2pjdeWBBgDuO5XfmGnoU0=
 golang.org/x/tools v0.0.0-20191114222411-4191b8cbba09/go.mod h1:b+2E5dAYhXwXZwtnZ6UAqBI28+e2cm9otk0dWdXHAEo=
 golang.org/x/xerrors v0.0.0-20190717185122-a985d3407aa7/go.mod h1:I/5z698sn9Ka8TeJc9MKroUUfqBBauWjQqLJ2OPfmY0=
-gopkg.in/check.v1 v0.0.0-20161208181325-20d25e280405 h1:yhCVgyC4o1eVCa2tZl7eS0r+SDo693bJlVdllGtEeKM=
 gopkg.in/check.v1 v0.0.0-20161208181325-20d25e280405/go.mod h1:Co6ibVJAznAaIkqp8huTwlJQCZ016jof/cbN4VW5Yz0=
 gopkg.in/yaml.v2 v2.2.2 h1:ZCJp+EgiOT7lHqUV2J862kp8Qj64Jo6az82+3Td9dZw=
 gopkg.in/yaml.v2 v2.2.2/go.mod h1:hI93XBmqTisBFMUTm0b8Fm+jr3Dg1NNxqwp+5A1VGuI=

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,153 @@
+package slippy
+
+import (
+	"math"
+
+	"github.com/go-spatial/geom"
+)
+
+/*
+ *
+ * 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
+
+	// floatVariance is used to compare floating point numbers, and to deal with float drift
+	//
+	//  This is mainly used to nudge the calculated tile values into place.
+	//  If the calculated number is extremely close to the border — basically on it — bumping it so that it falls
+	//    into the right most or bottom most tiles, where it should be. Since the variance from the floating point
+	//    calculation could be either positive or negative. If it's negative it would end up in the tile that is to the
+	//    left (leading possibility to a negative tile number) or top tile.
+	//  Take for example a calculated value of 6.99999999 v.s.7.000001 as a float these are practically the same number,
+	//    but we need them to be 7+, as we will truncate the float in the next step
+	//    (the fractional part is the percentage into the tile where the pixel is), thus it is better to bump toward 7
+	//    by the way of a very small step. This is the floatVariance value we use.
+	floatVariance = 0.000001
+)
+
+// 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
+}
+
+// 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) {
+	if tileSize == 0 {
+		tileSize = DefaultTileSize
+	}
+	// bound it because we have a top of the world problem
+	if lat < -Lat4326Max {
+		return int(z.N() - 1)
+	}
+
+	if lat > Lat4326Max {
+		return 0
+	}
+	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
+	}
+	worldSize := float64(tileSize) * z.N()
+
+	// Convert the Degree to radians as most of the math functions work in radians
+	radLat := Degree2Radians(45 + lat/2)
+	// normalize lat
+	latTan := math.Tan(radLat)
+	latNormalized := math.Log(latTan)
+
+	// compute the pixel value for y:
+	yPxRaw := (180 - Radians2Degree(latNormalized)) / 360
+	yPx = yPxRaw * worldSize
+	// instead of getting 7.0 we can end up with 6.9999999999, etc... use floatVariance to correct for such cases
+	return yPx + floatVariance
+}
+
+// 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
+	}
+
+	if lon <= -Lon4326Max {
+		return 0
+	}
+
+	if lon >= Lon4326Max {
+		return int(z.N() - 1)
+	}
+
+	tileX := lon2Px(tileSize, z, lon)
+	tileX = tileX / float64(tileSize)
+	// Truncate to get the tile
+	return int(tileX)
+
+}
+
+// lonPx will return the pixel coordinate for the lon. 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 lon2Px(tileSize uint32, z Zoom, lon float64) (xPx float64) {
+	if tileSize == 0 {
+		tileSize = DefaultTileSize
+	}
+	worldSize := float64(tileSize) * z.N()
+	lonNormalized := 180 + lon
+	// compute the pixel value for x:
+	xPxRaw := lonNormalized / 360
+	xPx = xPxRaw * worldSize
+	// instead of getting 7.0 we can end up with 6.9999999999, etc... use floatVariance to correct for such cases
+	return xPx + floatVariance
+}
+
+func PtFromLatLon(lat, lon float64) geom.Point {
+	return geom.Point{lon, 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
+}