Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 3 additions & 1 deletion planar/makevalid/makevalid_test.go
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,8 @@ import (
"strings"
"testing"

"github.com/go-spatial/proj"

"github.com/go-spatial/geom/encoding/wkt"
"github.com/go-spatial/geom/slippy"

Expand Down Expand Up @@ -132,7 +134,7 @@ func checkMakeValid(tb testing.TB) {
didClip: true,
},
"issue#70_full": {
ClipBox: slippy.NewTile(13, 8054, 2677).Extent3857().ExpandBy(64.0),
ClipBox: slippy.NewTile(13, 8054, 2677).Extent3857(proj.WebMercator).ExpandBy(64.0),
MultiPolygon: func() *geom.MultiPolygon {
b, err := ioutil.ReadFile(`testdata/issue70.polygon`)
if err != nil {
Expand Down
76 changes: 27 additions & 49 deletions slippy/projections.go
Original file line number Diff line number Diff line change
@@ -1,78 +1,56 @@
package slippy

import "math"
import (
"fmt"
"math"

// ==== lat lon (aka WGS 84) ====
"github.com/go-spatial/geom"
)

// Lat2Tile takes a zoom and a lat to produce the lon
func Lat2Tile(zoom uint, lat float64) (y uint) {
latRad := lat * math.Pi / 180

return uint(math.Exp2(float64(zoom))*
(1.0-math.Log(
math.Tan(latRad)+
(1/math.Cos(latRad)))/math.Pi)) /
2.0
}

// Lon2Tile takes in a zoom and lon to produce the lat
func Lon2Tile(zoom uint, lon float64) (x uint) {
return uint(math.Exp2(float64(zoom)) * (lon + 180.0) / 360.0)
type extents struct {
NativeExtents *geom.Extent
WGS84Extents *geom.Extent
Grid TileGrid
}

// Tile2Lon will return the west most longitude
func Tile2Lon(zoom, x uint) float64 { return float64(x)/math.Exp2(float64(zoom))*360.0 - 180.0 }

// Tile2Lat will return the north most latitude
func Tile2Lat(zoom, y uint) float64 {
var n float64 = math.Pi
if y != 0 {
n = math.Pi - 2.0*math.Pi*float64(y)/math.Exp2(float64(zoom))
}

return 180.0 / math.Pi * math.Atan(0.5*(math.Exp(n)-math.Exp(-n)))
// SupportedProjections contains supported projection native and lat/long extents as well as tile layout ratio
var SupportedProjections = map[uint]extents{
3857: extents{NativeExtents: &geom.Extent{-20026376.39, -20048966.10, 20026376.39, 20048966.10}, WGS84Extents: &geom.Extent{-180.0, -85.0511, 180.0, 85.0511}, Grid: GetGrid(3857)},
4326: extents{NativeExtents: &geom.Extent{-180.0, -90.0, 180.0, 90.0}, WGS84Extents: &geom.Extent{-180.0, -90.0, 180.0, 90.0}, Grid: GetGrid(4326)},
}

// ==== Web Mercator ====

// WebMercatorMax is the max size in meters of a tile
const WebMercatorMax = 20037508.34

// Tile2WebX returns the side of the tile in the -x side
func Tile2WebX(zoom uint, n uint) float64 {
// Tile2WebX returns the side of the tile in the -x side in webmercator
func Tile2WebX(zoom uint, n uint, srid uint) float64 {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The extra arg isn't really necessary here. If it's being used elsewhere to satisfy func (uint, uint, uint) float64 signature, a closure should be made there. Since there there isn't a type requirement in this package it shouldn't be forced here.

res := (WebMercatorMax * 2) / math.Exp2(float64(zoom))

return -WebMercatorMax + float64(n)*res
}

// Tile2WebY returns the side of the tile in the +y side
func Tile2WebY(zoom uint, n uint) float64 {
// Tile2WebY returns the side of the tile in the +y side in webmercator
func Tile2WebY(zoom uint, n uint, srid uint) float64 {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The extra arg isn't really necessary here. If it's being used elsewhere to satisfy func (uint, uint, uint) float64 signature, a closure should be made there. Since there there isn't a type requirement in this package it shouldn't be forced here.

res := (WebMercatorMax * 2) / math.Exp2(float64(zoom))

return WebMercatorMax - float64(n)*res
}

// WebX2Tile returns the column of the tile given the web mercator x value
func WebX2Tile(zoom uint, x float64) uint {
res := (WebMercatorMax * 2) / math.Exp2(float64(zoom))

return uint((x + WebMercatorMax) / res)
}

// WebY2Tile returns the row of the tile given the web mercator y value
func WebY2Tile(zoom uint, y float64) uint {
res := (WebMercatorMax * 2) / math.Exp2(float64(zoom))

return uint(-(y - WebMercatorMax) / res)
}

// ==== pixels ====

// MvtTileDim is the number of pixels in a tile
const MvtTileDim = 4096.0

// Pixels2Webs scalar conversion of pixels into web mercator units
// PixelsToProjectedUnits scalar conversion of pixels into projected units
// TODO (@ear7h): perhaps rethink this
func Pixels2Webs(zoom uint, pixels uint) float64 {
return WebMercatorMax * 2 / math.Exp2(float64(zoom)) * float64(pixels) / MvtTileDim
func PixelsToProjectedUnits(zoom uint, pixels uint, srid uint) float64 {
Copy link
Contributor

@ear7h ear7h Nov 26, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This could be renamed FromPixels Pixels2Native and have a func (TileGrid, zoom, pixels uint) float signature

switch srid {
case 3857:
return WebMercatorMax * 2 / math.Exp2(float64(zoom)) * float64(pixels) / MvtTileDim
case 4326:
return 360.0 / math.Exp2(float64(zoom)) * float64(pixels) / MvtTileDim / 2
default:
panic(fmt.Sprintf("unsupported srid: %v", srid))
}
}
83 changes: 62 additions & 21 deletions slippy/tile.go
Original file line number Diff line number Diff line change
@@ -1,9 +1,11 @@
package slippy

import (
"errors"
"fmt"
"math"

"errors"
"github.com/go-spatial/proj"

"github.com/go-spatial/geom"
)
Expand Down Expand Up @@ -33,15 +35,16 @@ type Tile struct {
// NewTileMinMaxer returns the smallest tile which fits the
// geom.MinMaxer. Note: it assumes the values of ext are
// EPSG:4326 (lng/lat)
func NewTileMinMaxer(ext geom.MinMaxer) *Tile {
upperLeft := NewTileLatLon(MaxZoom, ext.MaxY(), ext.MinX())
//TODO (meilinger): we need this anymore?
func NewTileMinMaxer(ext geom.MinMaxer, tileSRID uint) *Tile {
upperLeft := NewTileLatLon(MaxZoom, ext.MaxY(), ext.MinX(), tileSRID)
point := &geom.Point{ext.MaxX(), ext.MinY()}

var ret *Tile

for z := uint(MaxZoom); int(z) >= 0 && ret == nil; z-- {
upperLeft.RangeFamilyAt(z, func(tile *Tile) error {
if tile.Extent4326().Contains(point) {
upperLeft.RangeFamilyAt(z, tileSRID, func(tile *Tile, srid uint) error {
if tile.Extent4326(tileSRID).Contains(point) {
ret = tile
return errors.New("stop iter")
}
Expand All @@ -54,9 +57,10 @@ func NewTileMinMaxer(ext geom.MinMaxer) *Tile {
}

// NewTileLatLon instantiates a tile containing the coordinate with the specified zoom
func NewTileLatLon(z uint, lat, lon float64) *Tile {
x := Lon2Tile(z, lon)
y := Lat2Tile(z, lat)
func NewTileLatLon(z uint, lat, lon float64, srid uint) *Tile {
grid := GetGrid(srid)
x := grid.Lon2XIndex(z, lon)
y := grid.Lat2YIndex(z, lat)

return &Tile{
Z: z,
Expand All @@ -73,13 +77,16 @@ func minmax(a, b uint) (uint, uint) {
}

// FromBounds returns a list of tiles that make up the bound given. The bounds should be defined as the following lng/lat points [4]float64{west,south,east,north}
func FromBounds(bounds *geom.Extent, z uint) []Tile {
func FromBounds(bounds *geom.Extent, z uint, tileSRID uint) []Tile {
if bounds == nil {
return nil
}

minx, maxx := minmax(Lon2Tile(z, bounds[0]), Lon2Tile(z, bounds[2]))
miny, maxy := minmax(Lat2Tile(z, bounds[1]), Lat2Tile(z, bounds[3]))
grid := GetGrid(tileSRID)

minx, maxx := minmax(grid.Lon2XIndex(z, bounds[0]), grid.Lon2XIndex(z, bounds[2]))
miny, maxy := minmax(grid.Lat2YIndex(z, bounds[1]), grid.Lat2YIndex(z, bounds[3]))

// tiles := make([]Tile, (maxx-minx)*(maxy-miny))
var tiles []Tile
for x := minx; x <= maxx; x++ {
Expand All @@ -94,30 +101,64 @@ func FromBounds(bounds *geom.Extent, z uint) []Tile {
// ZXY returns back the z,x,y of the tile
func (t Tile) ZXY() (uint, uint, uint) { return t.Z, t.X, t.Y }

// Extent gets the extent of the tile in the units of the tileSRID
func (t Tile) NativeExtent(tileSRID uint) *geom.Extent {
if _, ok := SupportedProjections[tileSRID]; !ok {
panic(fmt.Sprintf("unsupported tileSRID %v", tileSRID))
}

grid := GetGrid(tileSRID)
pts := []float64{grid.XIndex2Lon(t.Z, t.X), grid.YIndex2Lat(t.Z, t.Y+1), grid.XIndex2Lon(t.Z, t.X+1), grid.YIndex2Lat(t.Z, t.Y)}

// No need to go further, we've already got the WGS84 extents
if tileSRID == proj.WGS84 {
return geom.NewExtent(
[2]float64{pts[0], pts[1]},
[2]float64{pts[2], pts[3]},
)
}

pts, err := proj.Convert(proj.EPSGCode(tileSRID), pts)
if err != nil {
panic(fmt.Sprintf("error converting %v to %v", pts, tileSRID))
}

return geom.NewExtent(
[2]float64{pts[0], pts[1]},
[2]float64{pts[2], pts[3]},
)
}

// Extent3857 returns the tile's extent in EPSG:3857 (aka Web Mercator) projection
func (t Tile) Extent3857() *geom.Extent {
func (t Tile) Extent3857(tileSRID uint) *geom.Extent {
if tileSRID != 3857 {
// Can't necessarily get webmercator extent for 4326 tile
panic("unable to get 3857 extent on 4326 tile")
}
return geom.NewExtent(
[2]float64{Tile2WebX(t.Z, t.X), Tile2WebY(t.Z, t.Y+1)},
[2]float64{Tile2WebX(t.Z, t.X+1), Tile2WebY(t.Z, t.Y)},
[2]float64{Tile2WebX(t.Z, t.X, tileSRID), Tile2WebY(t.Z, t.Y+1, tileSRID)},
[2]float64{Tile2WebX(t.Z, t.X+1, tileSRID), Tile2WebY(t.Z, t.Y, tileSRID)},
)
}

// Extent4326 returns the tile's extent in EPSG:4326 (aka lat/long)
func (t Tile) Extent4326() *geom.Extent {
// Extent4326 returns the tile's extent in EPSG:4326 (aka lat/long) given the tilespace's SRID
func (t Tile) Extent4326(tileSRID uint) *geom.Extent {
grid := GetGrid(tileSRID)
return geom.NewExtent(
[2]float64{Tile2Lon(t.Z, t.X), Tile2Lat(t.Z, t.Y+1)},
[2]float64{Tile2Lon(t.Z, t.X+1), Tile2Lat(t.Z, t.Y)},
[2]float64{grid.XIndex2Lon(t.Z, t.X), grid.YIndex2Lat(t.Z, t.Y+1)},
[2]float64{grid.XIndex2Lon(t.Z, t.X+1), grid.YIndex2Lat(t.Z, t.Y)},
)
}

// RangeFamilyAt calls f on every tile vertically related to t at the specified zoom
// TODO (ear7h): sibling support
func (t Tile) RangeFamilyAt(zoom uint, f func(*Tile) error) error {
//TODO (meilinger): should this be part of TileGrid?
func (t Tile) RangeFamilyAt(zoom, srid uint, f func(tile *Tile, srid uint) error) error {
// handle ancestors and self
if zoom <= t.Z {
mag := t.Z - zoom
arg := NewTile(zoom, t.X>>mag, t.Y>>mag)
return f(arg)
return f(arg, srid)
}

// handle descendants
Expand All @@ -129,7 +170,7 @@ func (t Tile) RangeFamilyAt(zoom uint, f func(*Tile) error) error {

for x := leastX; x < leastX+delta; x++ {
for y := leastY; y < leastY+delta; y++ {
err := f(NewTile(zoom, x, y))
err := f(NewTile(zoom, x, y), srid)
if err != nil {
return err
}
Expand Down
107 changes: 107 additions & 0 deletions slippy/tile_grid.go
Original file line number Diff line number Diff line change
@@ -0,0 +1,107 @@
package slippy

import (
"fmt"
"math"
)

// TileGrid contains the tile layout, including ability to get WGS84 coordinates for tile extents
type TileGrid interface {
Copy link
Contributor

@ear7h ear7h Nov 26, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

From my understanding the TileGrid is what defines the rules for the tiles and how they're placed. If this is correct, I think the following changes should be made:

  • the name should be changed to TileScheme which alludes to the rules which the tiles are created and laid out. "Grid" makes me thing of a structure which already exists (in memory); like a "dish" vs a "recipe"
  • the implementation should use the Tile type wherever possible. By doing this, the Tile and TileGrid types can be simplified and some functions moved to stand alone with func(TileGrid, Tile .... signatures. I'm thinking the new types and signatures would look like this:
type TileGrid interface {
    // returns the srid of the projection of this tiling scheme
    // the geom types from the other methods will be in this projection.
    SRID() int

    // this should return the (exclusive) maximum tile in this zoom. AKA:
    //    Tile{z, MaxX+1, MaxY+1}.
    // it is exclusive so this function can be easily used for boundary checks
    Size(z uint) Tile 

    // converts from a point (in the tile scheme's projection) and zoom to a tile.
    // ok will be false if the point is not valid in the projection
    FromNative(z uint, pt geom.Point) (t *Tile, ok bool)

    // Retuns the tile's upper left point. ok will be false if the tile is not valid.
    // <strike>one thing to note about implementing this is that you can compare tiles directly:
    //    ok := tileScheme.Max(tile.Z) >= tile</strike>
    // second thing to note is that the 
    ToNative(Tile) (pt geom.Point, ok bool)
}

// no need for the srid because f can be made a closure
func RangeFamilyAt(TileGrid, Tile, f func(*Tile) error) error

// Tile.ExtentXXXX aren't needed because of  ToNative
// NewTileLatLon aren't needed beacuse of FromNative

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I believe that we'd still need a mechanism to support getting a tile's extent in both projected as well as WGS84 units to support downstream usage within Tegola. Part of the reasoning for separating tiling logic out to support equirectangular tile schemes is in support of tile rendering as "4326" and a non-square tile grid. But, again, I very well could be missing some of the vision for these changes. I'd love to chat further so that I'm on the same page!

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I was thinking that the conversion from one coordinate system to another should be decoupled from slippy (it was not before, but it is something we were looking to change). The implementer of Grid works only with its native units and then the calling code can use the proj package to convert to WGS84 or other projection systems.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Gotcha. Makes sense.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ear7h re: you can compare tiles directly, I didn't realize that structs in golang offered an ordering mechanism (or operator overloading) but am having trouble a mechanism to support this construct. Is there something I'm missing or misunderstanding?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@meilinger My mistake, that doesn't work 😬 They can be compared for equality but not ordering.

ContainsIndex(zoom, x, y uint) bool
GridSize(zoom uint) (x, y uint)
MaxXY(zoom uint) (maxx, maxy uint)
Lat2YIndex(zoom uint, lat float64) (gridy uint)
Lon2XIndex(zoom uint, lon float64) (gridx uint)
XIndex2Lon(zoom, x uint) (lon float64)
YIndex2Lat(zoom, y uint) (lat float64)
}

func GetGrid(srid uint) TileGrid {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For convention, this function should be renamed NewGrid. Also, as a clarification, are there external limitations (other programs, a spec, convention etc.) which prevent the usage of square tiles with 4326 projection?

Copy link
Contributor Author

@meilinger meilinger Nov 29, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I believe that, from my experience, when someone mentions "4326" as the projection (which it is not), they are typically referring to equirectangular projection using lat/lon and WGS84 geoid upon a 2x1 tile ratio scheme (using square tiles, but 2x1 tile extent). But that is completely notional, not part of a spec and my experience is only a single data-point in usages.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think I understand what's going on bit better now, thank you for the explanation. I always get confused with srid, projections, and datums. From what I've found after a brief search it seem like some programs expect the behavior you describe. So, I think the functionality here is good.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I hear ya, I still get confused as to difference of them as well.

switch srid {
case 4326:
return &grid{tileExtentRatio: 2, srid: srid}
case 3857:
return &grid{tileExtentRatio: 1, srid: srid}
default:
panic(fmt.Sprintf("unsupported srid: %v", srid))
}
}

type grid struct {
tileExtentRatio float64
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you document what this ratio is? x/y or y/x

srid uint
}

func (g *grid) ContainsIndex(zoom, x, y uint) bool {
xsize, ysize := g.GridSize(zoom)
if x < xsize && y < ysize {
return true
}

return false
}

func (g *grid) GridSize(zoom uint) (x, y uint) {
dim := uint(math.Exp2(float64(zoom)))
return uint(float64(dim) * g.tileExtentRatio), dim
}

func (g *grid) MaxXY(zoom uint) (maxx, maxy uint) {
xsize, ysize := g.GridSize(zoom)

return xsize - 1, ysize - 1
}

func (g *grid) Lat2YIndex(zoom uint, lat float64) (gridy uint) {
switch g.srid {
case 3857:
latRad := lat * math.Pi / 180
return uint(math.Exp2(float64(zoom))*
(1.0-math.Log(
math.Tan(latRad)+
(1/math.Cos(latRad)))/math.Pi)) /
2.0
case 4326:
return uint(math.Exp2(float64(zoom)) * -(lat - 90.0) / (360.0 / g.tileExtentRatio))
default:
panic(fmt.Sprintf("unsupported srid: %v", g.srid))
}
}

func (g *grid) Lon2XIndex(zoom uint, lon float64) (gridx uint) {
switch g.srid {
case 3857:
fallthrough
case 4326:
return uint(math.Exp2(float64(zoom)) * (lon + 180.0) / (360.0 / g.tileExtentRatio))
default:
panic(fmt.Sprintf("unsupported srid: %v", g.srid))
}
}

func (g *grid) XIndex2Lon(zoom, x uint) (lon float64) {
switch g.srid {
case 3857:
fallthrough
case 4326:
return float64(x)/math.Exp2(float64(zoom))*(360.0/g.tileExtentRatio) - 180.0
default:
panic(fmt.Sprintf("unsupported srid: %v", g.srid))
}
}

func (g *grid) YIndex2Lat(zoom, y uint) (lat float64) {
switch g.srid {
case 3857:
var n = math.Pi
if y != 0 {
n = math.Pi - 2.0*math.Pi*float64(y)/math.Exp2(float64(zoom))
}

return 180.0 / math.Pi * math.Atan(0.5*(math.Exp(n)-math.Exp(-n)))
case 4326:
return -(180.0/math.Exp2(float64(zoom))*float64(y) - 90.0)
default:
panic(fmt.Sprintf("unsupported srid: %v", g.srid))
}
}
Loading