Fetching elevation data from Mapbox

In slippymath: Slippy Map Tile Tools

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

There are some slippy map tile providers that can serve you tiles that represent Digital Elevation Model (DEM) data, rather than map imagery. This is the kind of data you want to make 3D maps with rayshader or quadmesh. Mapbox have an API for DEM tiles called Mapbox Terrain-RGB that we will use in this example.

Get a Mapbox API key

Sign up for Mapbox and generate yourself an API access token. For my testing I have used a 'public' token with styles:tiles. Add that API token to your .Renviron file (in your home directory) as:

MAPBOX_API_KEY=<YOUR REALLY LONG TOKEN HERE>

Fetch the RGB tiles for your bounding box

How many tiles?

library(slippymath)

tibrogargan<- c(xmin = 152.938485, ymin = -26.93345, xmax = 152.956467, 
               ymax = -26.921463)

slippymath::bbox_tile_query(tibrogargan)

It's a small area so we don't need a lot. 9 tiles as zoom 15 looks good. Let's get the tile grid.

tibrogargan_grid <- bbox_to_tile_grid(tibrogargan,zoom = 15)

Now we'll fetch the tiles from Mapbox.

library(glue)
library(purrr)
library(curl)

mapbox_query_string <-
paste0("https://api.mapbox.com/v4/mapbox.terrain-rgb/{zoom}/{x}/{y}.jpg90",
       "?access_token=",
       Sys.getenv("MAPBOX_API_KEY"))

tibro_tiles <-
pmap(.l = tibrogargan_grid$tiles,
     zoom = tibrogargan_grid$zoom,

     .f = function(x, y, zoom){
       outfile <- glue("{x}_{y}.jpg")
       curl_download(url = glue(mapbox_query_string),
            destfile = outfile)
       outfile
     }
     )

Stitch tiles into a raster

We composite the images to raster with compose_tile_grid and view the raster:

tibrogargan_raster <- slippymath::compose_tile_grid(tibrogargan_grid, tibro_tiles)
raster::plot(tibrogargan_raster)

r knitr::include_graphics("tibrogargan_layers.png")

We see this mix of layers including a psychedelic blue-green RGB landscape because in terrain-rgb tiles, the RGB values to provide additional precision to the elevation information. We'll decode this in the next section.

Converting RGB tiles to DEM tiles

Mapbox provide a formula to decode the RGB values to elevation. This was implemented in the function below by Michael Sumner.

decode_elevation <- function(dat) {
  height <-  -10000 + ((dat[[1]] * 256 * 256 + dat[[2]] * 256 + dat[[3]]) * 0.1)
  raster::projection(height) <- "+proj=merc +a=6378137 +b=6378137"
  height
}

When we apply it to tibrogargan_raster we get:

tibrogargan_elevation <- decode_elevation(tibrogargan_raster)
raster::plot(tibrogargan_elevation)

r knitr::include_graphics("tibrogargan_decoded.png")

Rendering DEM image

Rayshader

library(magrittr)
library(rayshader)

elevation_mat <- t(raster::as.matrix(tibrogargan_elevation))

shadow_mat <- ray_shade(elevation_mat)

elevation_mat %>%
  sphere_shade(progbar = FALSE, sunangle = 45) %>%
  add_shadow(shadow_mat) %>%
  plot_3d(elevation_mat,
          zscale = 7,
          phi = 30,
          theta = 135)

r knitr::include_graphics("tibrogargan_rayshade.png")

Quadmesh

library(quadmesh)

elevation_mesh <- quadmesh(tibrogargan_elevation)
rgl::shade3d(elevation_mesh, col = "light green")

r knitr::include_graphics("tibrogargan_quadmesh_rgl.png")

See Also

*ceramic - A wrapper for quadmesh + slippymath for making 3D surfaces textured with satellite tiles.

slippymath documentation built on June 28, 2019, 5:04 p.m.