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.
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.
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.
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.
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.
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.
Finally a convenience function is provided to convert the raster to a .png if needed: