How to Slippymath

Miles McBain

2019-06-29

slippymath provides tools to assist your work with slippy map tile servers and the map tiles they provide. The main things it can help with are: * determining what zoom level to use * identifying the tiles you need to plot your data over * identifying where lonlat coordinates fall in tile space * determining the bounding boxes of tiles * compositing tiles into a single spatially referenced raster.

Downloading tiles does not get a built-in helper. This is to make slippymath as service agnostic as possible. There is an example of downloading from Mapbox in the README. You may also enjoy the ceramic package by Michael Sumner.

Determining zoom level

The main assumption slippymath makes is that you can calculate a bounding box for your spatial data. Most spatial tools have this capability - see sf::st_bbox().

Given a bounding box in longitude and latitude (EPSG: 4326), you can execute a ‘tile query’ to see how many slippy map tiles your data would occupy at various zoom levels:

library(slippymath)

uluru_bbox <-
    c(xmin = 131.02084,
      xmax = 131.0535,
      ymin = -25.35461,
      ymax = -25.33568)

bbox_tile_query(uluru_bbox)
#>     x_min  y_min  x_max  y_max y_dim x_dim total_tiles zoom
#> 1       3      2      3      2     1     1           1    2
#> 2       6      4      6      4     1     1           1    3
#> 3      13      9     13      9     1     1           1    4
#> 4      27     18     27     18     1     1           1    5
#> 5      55     36     55     36     1     1           1    6
#> 6     110     73    110     73     1     1           1    7
#> 7     221    146    221    146     1     1           1    8
#> 8     442    293    442    293     1     1           1    9
#> 9     884    586    884    586     1     1           1   10
#> 10   1769   1173   1769   1173     1     1           1   11
#> 11   3538   2346   3539   2346     1     2           2   12
#> 12   7077   4692   7078   4692     1     2           2   13
#> 13  14154   9384  14156   9385     2     3           6   14
#> 14  28309  18769  28312  18771     3     4          12   15
#> 15  56619  37538  56625  37542     5     7          35   16
#> 16 113239  75076 113251  75084     9    13         117   17
#> 17 226478 150153 226502 150168    16    25         400   18

When choosing your zoom level consider: Tiles are 256 x 256 pixel squares which roughly range in size from 30 - 40kb. 100’s of tiles will be 10’s of megabytes. You may also want to check rate limits imposed by the server you intend to pull the tiles from.

slippymath can also choose a zoom level for your tiles, given a tile budget. See below.

Identifying tiles to plot over

slippymath calls the set of tiles you need to plot all your data over a “Tile Grid” or tg for short. You can get a tile grid given a bounding box and tile zoom level OR a bounding box and tile budget (the zoom will be chosen to fit the budget).

For example using the Uluru bounding box from above and zoom level 14:

bbox_to_tile_grid(uluru_bbox, zoom = 14)
#> $tiles
#>       x    y
#> 1 14154 9384
#> 2 14155 9384
#> 3 14156 9384
#> 4 14154 9385
#> 5 14155 9385
#> 6 14156 9385
#> 
#> $zoom
#> [1] 14
#> 
#> attr(,"class")
#> [1] "tile_grid"

And using the same bounding box with a budget of 15 tiles:

bbox_to_tile_grid(uluru_bbox, max_tiles = 15)
#> $tiles
#>        x     y
#> 1  28309 18769
#> 2  28310 18769
#> 3  28311 18769
#> 4  28312 18769
#> 5  28309 18770
#> 6  28310 18770
#> 7  28311 18770
#> 8  28312 18770
#> 9  28309 18771
#> 10 28310 18771
#> 11 28311 18771
#> 12 28312 18771
#> 
#> $zoom
#> [1] 15
#> 
#> attr(,"class")
#> [1] "tile_grid"

A tile grid is an object that contains a dataframe of tile coordinates ($tiles) and a $zoom level.

Identifying the tile for a single point

Occasionally it may make sense to find the tile at a given zoom for a single point. slippymath has functions to convert between longitude and latitude tile coordinate space:

lonlat_to_tilenum(131.02084, -25.35461, zoom = 14)
#> $x
#> [1] 14154
#> 
#> $y
#> [1] 9385

tilenum_to_lonlat(14154, 9385, zoom = 14)
#> $lon
#> [1] 131.001
#> 
#> $lat
#> [1] -25.34403

In this example, the exact lon lat we started with is not returned in the conversion from tile coordinate space. The lon lat returned is always the top left corner point of the tile.

Determining tile bounding boxes

Once you have tiles, you will likely want to know their spatial bounding boxes as this is essential for compositing them. If you use the built-in compositing function, this is taken care of for you.

Find the bounding boxes of a tile grid like so:

uluru_grid <- bbox_to_tile_grid(uluru_bbox, max_tiles = 15)
tile_grid_bboxes(uluru_grid)
#> [[1]]
#>     xmin     ymin     xmax     ymax 
#> 14584185 -2918060 14585408 -2916837 
#> attr(,"class")
#> [1] "bbox"
#> attr(,"crs")
#> $epsg
#> [1] 3857
#> 
#> $proj4string
#> [1] "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
#> 
#> attr(,"class")
#> [1] "crs"
#> 
#> [[2]]
#>     xmin     ymin     xmax     ymax 
#> 14585408 -2918060 14586631 -2916837 
#> attr(,"class")
#> [1] "bbox"
#> attr(,"crs")
#> $epsg
#> [1] 3857
#> 
#> $proj4string
#> [1] "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
#> 
#> attr(,"class")
#> [1] "crs"
#> 
#> [[3]]
#>     xmin     ymin     xmax     ymax 
#> 14586631 -2918060 14587854 -2916837 
#> attr(,"class")
#> [1] "bbox"
#> attr(,"crs")
#> $epsg
#> [1] 3857
#> 
#> $proj4string
#> [1] "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
#> 
#> attr(,"class")
#> [1] "crs"
#> 
#> [[4]]
#>     xmin     ymin     xmax     ymax 
#> 14587854 -2918060 14589077 -2916837 
#> attr(,"class")
#> [1] "bbox"
#> attr(,"crs")
#> $epsg
#> [1] 3857
#> 
#> $proj4string
#> [1] "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
#> 
#> attr(,"class")
#> [1] "crs"
#> 
#> [[5]]
#>     xmin     ymin     xmax     ymax 
#> 14584185 -2919283 14585408 -2918060 
#> attr(,"class")
#> [1] "bbox"
#> attr(,"crs")
#> $epsg
#> [1] 3857
#> 
#> $proj4string
#> [1] "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
#> 
#> attr(,"class")
#> [1] "crs"
#> 
#> [[6]]
#>     xmin     ymin     xmax     ymax 
#> 14585408 -2919283 14586631 -2918060 
#> attr(,"class")
#> [1] "bbox"
#> attr(,"crs")
#> $epsg
#> [1] 3857
#> 
#> $proj4string
#> [1] "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
#> 
#> attr(,"class")
#> [1] "crs"
#> 
#> [[7]]
#>     xmin     ymin     xmax     ymax 
#> 14586631 -2919283 14587854 -2918060 
#> attr(,"class")
#> [1] "bbox"
#> attr(,"crs")
#> $epsg
#> [1] 3857
#> 
#> $proj4string
#> [1] "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
#> 
#> attr(,"class")
#> [1] "crs"
#> 
#> [[8]]
#>     xmin     ymin     xmax     ymax 
#> 14587854 -2919283 14589077 -2918060 
#> attr(,"class")
#> [1] "bbox"
#> attr(,"crs")
#> $epsg
#> [1] 3857
#> 
#> $proj4string
#> [1] "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
#> 
#> attr(,"class")
#> [1] "crs"
#> 
#> [[9]]
#>     xmin     ymin     xmax     ymax 
#> 14584185 -2920506 14585408 -2919283 
#> attr(,"class")
#> [1] "bbox"
#> attr(,"crs")
#> $epsg
#> [1] 3857
#> 
#> $proj4string
#> [1] "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
#> 
#> attr(,"class")
#> [1] "crs"
#> 
#> [[10]]
#>     xmin     ymin     xmax     ymax 
#> 14585408 -2920506 14586631 -2919283 
#> attr(,"class")
#> [1] "bbox"
#> attr(,"crs")
#> $epsg
#> [1] 3857
#> 
#> $proj4string
#> [1] "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
#> 
#> attr(,"class")
#> [1] "crs"
#> 
#> [[11]]
#>     xmin     ymin     xmax     ymax 
#> 14586631 -2920506 14587854 -2919283 
#> attr(,"class")
#> [1] "bbox"
#> attr(,"crs")
#> $epsg
#> [1] 3857
#> 
#> $proj4string
#> [1] "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
#> 
#> attr(,"class")
#> [1] "crs"
#> 
#> [[12]]
#>     xmin     ymin     xmax     ymax 
#> 14587854 -2920506 14589077 -2919283 
#> attr(,"class")
#> [1] "bbox"
#> attr(,"crs")
#> $epsg
#> [1] 3857
#> 
#> $proj4string
#> [1] "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
#> 
#> attr(,"class")
#> [1] "crs"

## for a single tile use tg_bb()

Slippy maps use a special mercator projection: EPSG: 3857. The tile bounding boxes are expressed in this projection.

Compositing tiles to a spatial raster

A simple function to composite tiles to a single raster with the correct spatial extent (again in the EPSG 3857 projection) is provided. It constructs individual raster and subsequently calls raster::merge which has issues with manual calls to gc which can slow things down with many tiles.

To composite, it is assumed you have a tile grid object, and a list of image files that correspond in order to the tiles in the tile grid. See the README example for a way to download tiles and get such a list.

raster_out <- compose_tile_grid(tile_grid, images)

Finally a convenience function is provided to convert the raster to a .png if needed:

raster_to_png(raster_out, "uluru.png")