Fetching elevation data from Mapbox

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)
#>     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       7      4      7      4     1     1           1    3
#> 3      14      9     14      9     1     1           1    4
#> 4      29     18     29     18     1     1           1    5
#> 5      59     36     59     36     1     1           1    6
#> 6     118     73    118     73     1     1           1    7
#> 7     236    147    236    147     1     1           1    8
#> 8     473    295    473    295     1     1           1    9
#> 9     947    591    947    591     1     1           1   10
#> 10   1894   1183   1894   1183     1     1           1   11
#> 11   3788   2366   3788   2366     1     1           1   12
#> 12   7576   4732   7576   4732     1     1           1   13
#> 13  15152   9464  15153   9465     2     2           4   14
#> 14  30304  18929  30306  18931     3     3           9   15
#> 15  60609  37859  60612  37862     4     4          16   16
#> 16 121219  75719 121225  75724     6     7          42   17
#> 17 242438 151439 242451 151449    11    14         154   18

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)

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)

Rendering DEM image

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)

library(quadmesh)

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

See Also

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