GDAL and image tiles, the {ceramic} package

news
code
Author

Michael D. Sumner

Published

2023-04-22

A new version of {ceramic} is now on CRAN, version 0.8.0.

The package exists for two purpose

  1. to download tiles from image servers
  2. to load raster data from image servers

NOTE: we need a Mapbox key for the Mapbox servers see, please see ceramic::get_api_key().

The original versions of ceramic didn’t really separate these tasks but the new version does.

Get raster data from online

We can read satellite imagery or elevation with a central point and a buffer (in metres).

library(ceramic)
Loading required package: terra
terra 1.7.23
pt <- cbind(147.3257, -42.8826)

(im <- cc_location(pt, buffer = c(15000, 25000)))
class       : SpatRaster 
dimensions  : 960, 576, 3  (nrow, ncol, nlyr)
resolution  : 52.08333, 52.08333  (x, y)
extent      : 16385222, 16415222, -5319119, -5269119  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857) 
source(s)   : memory
colors RGB  : 1, 2, 3 
names       : red, green, blue 
min values  :   0,     0,    0 
max values  : 255,   255,  252 
(el <- cc_elevation(pt, buffer = c(15000, 25000)))
class       : SpatRaster 
dimensions  : 960, 576, 1  (nrow, ncol, nlyr)
resolution  : 52.08333, 52.08333  (x, y)
extent      : 16385222, 16415222, -5319119, -5269119  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857) 
source(s)   : memory
name        :  lyr.1 
min value   :   -1.9 
max value   : 1264.9 
op <- par(mfcol = c(1, 2))
plot(el, legend = F); plot(im, add = TRUE); plot(el, legend = F, 
                                                 col = hcl.colors(64)[tail(seq_len(64), 45)])

par(op)

These raster objects are in terra ‘SpatRaster’ format (older ceramic used raster package).

These use the Mapbox ‘mapbox.satellite’ and ‘mapbox-terrain-rgb’ tile servers.

Get raster tiles from online

If we want the actual tiles, we can use the original get_tiles() function. (In older versions cc_location() and cc_elevation() would invoke get_tiles, but no longer).

pt <- cbind(147.3257, -42.8826)

imtiles <- get_tiles(pt, buffer = c(15000, 25000))
Preparing to download: 24 tiles at zoom = 12 from 
https://api.mapbox.com/v4/mapbox.satellite/
str(imtiles)
List of 3
 $ files : chr [1:24] "/perm_storage/home/mdsumner/.cache/.ceramic/api.mapbox.com/v4/mapbox.satellite/12/3722/2586.jpg" "/perm_storage/home/mdsumner/.cache/.ceramic/api.mapbox.com/v4/mapbox.satellite/12/3723/2586.jpg" "/perm_storage/home/mdsumner/.cache/.ceramic/api.mapbox.com/v4/mapbox.satellite/12/3724/2586.jpg" "/perm_storage/home/mdsumner/.cache/.ceramic/api.mapbox.com/v4/mapbox.satellite/12/3725/2586.jpg" ...
 $ tiles :List of 2
  ..$ tiles:'data.frame':   24 obs. of  2 variables:
  .. ..$ x: int [1:24] 3722 3723 3724 3725 3722 3723 3724 3725 3722 3723 ...
  .. ..$ y: int [1:24] 2586 2586 2586 2586 2587 2587 2587 2587 2588 2588 ...
  .. ..- attr(*, "out.attrs")=List of 2
  .. .. ..$ dim     : Named int [1:2] 4 6
  .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
  .. .. ..$ dimnames:List of 2
  .. .. .. ..$ x: chr [1:4] "x=3722" "x=3723" "x=3724" "x=3725"
  .. .. .. ..$ y: chr [1:6] "y=2586" "y=2587" "y=2588" "y=2589" ...
  ..$ zoom : int 12
  ..- attr(*, "class")= chr "tile_grid"
 $ extent: num [1:4] 16385222 16415222 -5319119 -5269119
eltiles <- get_tiles(pt, buffer = c(15000, 25000))
Preparing to download: 24 tiles at zoom = 12 from 
https://api.mapbox.com/v4/mapbox.satellite/
head(gsub(ceramic::ceramic_cache(), "",  eltiles$files))
[1] "/api.mapbox.com/v4/mapbox.satellite/12/3722/2586.jpg"
[2] "/api.mapbox.com/v4/mapbox.satellite/12/3723/2586.jpg"
[3] "/api.mapbox.com/v4/mapbox.satellite/12/3724/2586.jpg"
[4] "/api.mapbox.com/v4/mapbox.satellite/12/3725/2586.jpg"
[5] "/api.mapbox.com/v4/mapbox.satellite/12/3722/2587.jpg"
[6] "/api.mapbox.com/v4/mapbox.satellite/12/3723/2587.jpg"

And to see what tiles we have we can materialize their footprint in wk::rct() form, this is way more efficient than having to use polygons.

zoomtiles <- ceramic::ceramic_tiles(imtiles$tiles$zoom)
## sub out the ones we just triggered
zoomtiles <- dplyr::filter(zoomtiles, fullname %in% imtiles$files)
rc <- ceramic::tiles_to_polygon(zoomtiles)
plot(rc)
plot(ext(im), add = TRUE)  ## see the image from above is not the tiles, but the buffer around a point

plot(read_tiles(pt, buffer = c(15000, 25000)))  ## but we can read thos tiles exactly
Preparing to download: 24 tiles at zoom = 12 from 
https://api.mapbox.com/v4/mapbox.satellite/
plot(rc, add = TRUE, border = "white")

Diverse query inputs

Finally, we can use a point and buffer to get imagery, or we can use a spatial object, currently supported are objects from {geos}, {wk}, {terra}, {sf}, {sp}, {raster}, and {stars}.

Here’s an example.

wkarea <- wk::rct(xmin = 5e6, ymin = -2e6, xmax = 2e6, ymax = 3e6, crs = "+proj=laea +lon_0=25 +lat_0=31")
im2 <- cc_location(wkarea)
el2 <- cc_elevation(wkarea)
op <- par(mfcol = c(1, 2))
plot(el2, legend = F); plot(im2, add = TRUE); plot(el2, col = grey.colors(128), legend = FALSE)

par(op)

(Please note that the result is still in Mercator, the query can be in any projection but we’re not matching that here, only its extent - ceramic version 0.8.0 is merely a stepping stone to some of the things we can do better with GDAL).

Previous versions of cc_location() and cc_elevation() were stuck with only ‘zoom’ and ‘max_tiles’ arguments, these make sense when you restrict exactly to the available tiles but when you just want a given area and a resolution, ‘dimension’ is more appropriate.

By default, the dimension is chosen relative to the graphics device. But, we can also specify exactly what we want. (Use zero for one of the dimensions to let the system figure out an appropriate size for the query you have).

cc_location(wkarea)
class       : SpatRaster 
dimensions  : 960, 776, 3  (nrow, ncol, nlyr)
resolution  : 7094.906, 7094.906  (x, y)
extent      : 4203326, 9708973, 105375.6, 6916485  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857) 
source(s)   : memory
colors RGB  : 1, 2, 3 
names       : red, green, blue 
min values  :   0,     0,    0 
max values  : 255,   255,  255 
cc_location(wkarea, dimension = c(400, 700))
class       : SpatRaster 
dimensions  : 700, 400, 3  (nrow, ncol, nlyr)
resolution  : 13766.3, 9730.157  (x, y)
extent      : 4203326, 9709847, 105375.6, 6916485  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857) 
source(s)   : memory
colors RGB  : 1, 2, 3 
names       : red, green, blue 
min values  :   0,     2,    0 
max values  : 255,   255,  255 
cc_location(wkarea, dimension = c(1024, 0))
class       : SpatRaster 
dimensions  : 1267, 1024, 3  (nrow, ncol, nlyr)
resolution  : 5377.463, 5377.463  (x, y)
extent      : 4203326, 9709847, 103240.4, 6916485  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857) 
source(s)   : memory
colors RGB  : 1, 2, 3 
names       : red, green, blue 
min values  :   0,     1,    0 
max values  : 255,   255,  255 

Some other improvements

  • it’s faster, for loading reads data from image servers directly with the GDAL warper API
  • we have separation of tile downloading and interaction from raster loading via GDAL
  • we don’t materialize tiles as polygons, we have the tile description, or the compact rct representation
  • new functions read_tiles(), unpack_rgb(),

See the full list of changes here: https://hypertidy.github.io/ceramic/news/index.html.

ceramic is my final CRAN package that had direct dependencies on rgdal, rgeos, or maptools - in some ways it was the most challenging for me and had loomed as a problem for some time. But, I’ve progressed my own tool kit well in the time and I learnt a lot with this update.

A future version will probably make the separation (are we tiles, or are we loading raster?) more complete, at any rate there are better ways to organize these things now. There are data sources, data readers and writers, data structures, and algorithms. This work is part of a family of tools that aims to make that orthogonality real and accessible.