Skip to content

Commit

Permalink
Updated slippy tiles to work as it use to.
Browse files Browse the repository at this point in the history
But kept the ability to have different SRID's.
  • Loading branch information
gdey committed Jul 25, 2024
1 parent c26cad9 commit 1b94fbb
Show file tree
Hide file tree
Showing 20 changed files with 1,244 additions and 963 deletions.
16 changes: 12 additions & 4 deletions bbox.go
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,13 @@ type Extenter interface {
Extent() (extent [4]float64)
}

type PtMinMaxer interface {
// Min returns the minimum x and y values
Min() [2]float64
// Max returns the maximum x and y values
Max() [2]float64
}

// MinMaxer is a wrapper for an Extent that gets min/max of the extent
type MinMaxer interface {
MinX() float64
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand Down
39 changes: 30 additions & 9 deletions cmd/main.go
Original file line number Diff line number Diff line change
Expand Up @@ -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 {
Expand Down Expand Up @@ -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)
Expand All @@ -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.TileGriddor, 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 {
Expand All @@ -116,17 +127,21 @@ func main() {
fmt.Fprintf(os.Stderr, "Unabled to parse z: %v", err)
usage()
}
x, err := strconv.ParseUint(parts[1], 10, 64)
x, err := strconv.ParseInt(parts[1], 10, 64)
if err != nil {
fmt.Fprintf(os.Stderr, "Unabled to parse x: %v", err)
usage()
}
y, err := strconv.ParseUint(parts[2], 10, 64)
y, err := strconv.ParseInt(parts[2], 10, 64)
if err != nil {
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: int(x),
Y: int(y),
}
fileTemplate := newOutFile(tile, *tag)
geo, err := readInputWKT(flag.Args()[1])
if err != nil {
Expand All @@ -142,21 +157,27 @@ func main() {
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)
Expand Down
6 changes: 6 additions & 0 deletions point.go
Original file line number Diff line number Diff line change
Expand Up @@ -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] }

Expand Down
143 changes: 143 additions & 0 deletions slippy/maths.go
Original file line number Diff line number Diff line change
@@ -0,0 +1,143 @@
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

// floatDrift is used to compare floating point numbers, and to deal with float drift
floatDrift = 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 floatDrift to correct for such cases
return yPx + floatDrift
}

// 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 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 Zoom, lon float64) (xPx float64) {
if tileSize == 0 {
tileSize = DefaultTileSize
}
worldSize := float64(tileSize) * z.N()
lonNormalized := 180 + lon
// compute the pixel value for y:
xPxRaw := lonNormalized / 360
xPx = xPxRaw * worldSize
// instead of getting 7.0 we can end up with 6.9999999999, etc... use floatDrift to correct for such cases
return xPx + floatDrift
}

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
}
82 changes: 82 additions & 0 deletions slippy/maths_openstreet_maps_algo.go.old
Original file line number Diff line number Diff line change
@@ -0,0 +1,82 @@
package slippy

import (
"fmt"
"math"
)

// Algo from
// from ref: https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#ECMAScript_(JavaScript/ActionScript,_etc.)

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
}

func deg2num(z Zoom, lon, lat float64) (x int, y int) {
return degLon2x(z, lon), degLat2y(z, lat)
}

func num2deg(z Zoom, x, y int) (lat float64, long float64) {
lat = y2deg(z, y)
long = x2deg(z, x)
return lat, long
}

func degLat2y(z Zoom, lat float64) (y int) {
return int(deg2y(z, lat))
}

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 deg2x(z Zoom, lon float64) float64 {
n := math.Exp2(float64(z))
return math.Floor((lon + 180.0) / 360.0 * n)
}

func deg2y(z Zoom, lat float64) float64 {
n := z.N()

// bound it because we have a top of the world problem
if lat < -85.0511 {
return n - 1
}

if lat > 85.0511 {
return 0
}

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
}

Loading

0 comments on commit 1b94fbb

Please sign in to comment.