---
title: "Fetching elevation data from Mapbox"
author: "Miles McBain"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Fetching Elevation data from Mapbox}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
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](https://docs.mapbox.com/help/troubleshooting/access-elevation-data/#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?

```{r}
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.

```{r}
tibrogargan_grid <- bbox_to_tile_grid(tibrogargan,zoom = 15)
```

Now we'll fetch the tiles from Mapbox. 

```{r, eval=FALSE}
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:

```{r, eval=FALSE}
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](https://docs.mapbox.com/help/troubleshooting/access-elevation-data/#decode-data) to decode the RGB values to elevation. This was implemented in the function below by Michael Sumner.

```{r, eval=FALSE}

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:

```{r, eval=FALSE}
tibrogargan_elevation <- decode_elevation(tibrogargan_raster)
raster::plot(tibrogargan_elevation)

```

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

## Rendering DEM image

### Rayshader

```{r, eval=FALSE}
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

```{r, eval=FALSE}
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](https://www.github.com/hypertidy/ceramic) - A wrapper for `quadmesh` + `slippymath` for making 3D surfaces textured with satellite tiles.