Code to make pretty bike ride plots with rayshader and rStrava

3d_ride.R

source("helpers.R")

plot_3d_route_area <- function(stream, zscale = 15){
  
  # create bounding box for the ride
  ui("1 of 5: Locating Route...")
  bbox <- list(
    p1 = list(long = min(stream$lng), lat = min(stream$lat)),
    p2 = list(long = max(stream$lng), lat = max(stream$lat))
  )
  
  # get elevation data for the box
  ui("2 of 5: Getting Route Elevation Data...")
  image_size <- define_image_size(bbox, major_dim = 600)
  file <- get_usgs_elevation_data(bbox, size = image_size$size)
  
  elev_img <- raster::raster(file)
  elev_matrix <- matrix(
    raster::extract(elev_img, raster::extent(elev_img), buffer = 1000), 
    nrow = ncol(elev_img), ncol = nrow(elev_img)
  )
  
  # calculate rayshader layers
  ui("3 of 5: Calculating Rayshader Layers...")
  ambmat <- ambient_shade(elev_matrix, zscale = 30)
  raymat <- ray_shade(elev_matrix, zscale = 30, lambert = TRUE)
  watermap <- detect_water(elev_matrix)
  
  # get the overlay map
  ui("4 of 5: Overlaying Street Data...")
  file <- get_arcgis_map_image(bbox, 
                       map_type = "World_Street_Map",
                       width = image_size$width, 
                       height = image_size$height)
  
  overlay_img <- png::readPNG(file)
  
  # create the plot
  ui("5 of 5: Opening RGL Device...")
  rgl::clear3d()
  elev_matrix %>% 
    sphere_shade(texture = "imhof4") %>% 
    add_water(watermap, color = "imhof4") %>%
    add_overlay(overlay_img, alphalayer = 0.5) %>%
    add_shadow(raymat, max_darken = 0.5) %>%
    add_shadow(ambmat, max_darken = 0.5) %>%
    plot_3d(elev_matrix, zscale = zscale, windowsize = c(1200, 1000),
            water = TRUE, wateralpha = 0,
            theta = 25, phi = 30, zoom = 0.65, fov = 60)
  
  ui("You should now see a map. Call render_route with the results of this function to view your ride!")
  list(
    elev_matrix = elev_matrix,
    stream = stream,
    bbox = bbox,
    image_size = image_size,
    zscale = zscale
  )
}

render_route <- function(route,
                         sample = 10){
  
  # render the route as labels
  
  # calculate nudge amount
  znudge <- 0.005*(max(route$elev_matrix) - min(route$elev_matrix))
  
  pb <- progress::progress_bar$new(total = nrow(route$stream)/sample)
  for (i in seq(1, nrow(route$stream), by = sample)) {
    # translate lat/lon into x,y within this box
    # TODO: This should be more efficient!
    pb$tick()
    pos <- find_image_coordinates(
        long = route$stream$lng[i],
        lat = route$stream$lat[i],
        bbox = route$bbox,
        image_width = route$image_size$width,
        image_height = route$image_size$height
    )
    
    
    # add a dot tracking progress
    render_label(
      route$elev_matrix,
      x = pos$x,
      y = pos$y,
      z = znudge,
      zscale = route$zscale,
      relativez = TRUE,
      offset = 0,
      alpha = 0,
      textsize = 2,
      text = "."
    )
  }

}
  
ui_check <- function(...) {
  phrase <- do.call(paste, list(...))
  cat(paste(crayon::green(clisymbols::symbol$tick), phrase, sep = " "))
  cat("\n")
}

ui <- function(phrase) {
  cat(paste(crayon::red(clisymbols::symbol$circle)), phrase, sep = " ")
  cat("\n")
} 
  

example.R

# get ready
options(rgl.useNULL = FALSE)
library(rStrava)
library(rayshader)
library(rgdal)
library(leaflet)
source('3d_ride.R')

# setup strava client
client <- Sys.getenv("CLIENT")
secret <- Sys.getenv("SECRET")
token <- strava_oauth("Rayshader", client,secret, app_scope = "view_private")
stoken <- httr::config(token = token)

# get most recent strava activity
acts <- get_activity_list(stoken)
stream <- get_activity_streams(acts, stoken, acts = 1)

# plot it
get_activity_streams(acts, stoken, acts = 5) %>% 
 # this generates a 3d map of the activity area 
 plot_3d_route_area(zscale = 10) %>% 
 # this adds dots representing your movement
 render_route()
 
> ◯ 1 of 5: Locating Route...
> ◯ 2 of 5: Getting Route Elevation Data...
> ◯ 3 of 5: Calculating Rayshader Layers...
> ◯ 4 of 5: Overlaying Street Data...                                  
> ◯ 5 of 5: Opening RGL Device...
> ◯ You should now see a map. Call render_route with the results of this function to view your ride!

helpers.R

# Everything here here stolen from https://github.com/wcmbishop/rayshader-demo/blob/master/R/rayshader-gif.R

#' Download a map image from the ArcGIS REST API
#'
#' @param bbox bounding box coordinates (list of 2 points with long/lat values)
#' @param map_type map type to download - options are World_Street_Map, World_Imagery, World_Topo_Map
#' @param file file path to save to. Default is NULL, which will create a temp file.
#' @param width image width (pixels)
#' @param height image height (pixels)
#' @param sr_bbox Spatial Reference code for bounding box
#' 
#' @details This function uses the ArcGIS REST API, specifically the 
#' "Execute Web Map Task" task. You can find links below to a web UI for this
#' rest endpoint and API documentation.
#' 
#' Web UI: https://utility.arcgisonline.com/arcgis/rest/services/Utilities/PrintingTools/GPServer/Export%20Web%20Map%20Task/execute
#' API docs: https://developers.arcgis.com/rest/services-reference/export-web-map-task.htm
#'
#' @return file path for the downloaded .png map image
#'
#' @examples
#' bbox <- list(
#'   p1 = list(long = -122.522, lat = 37.707),
#'   p2 = list(long = -122.354, lat = 37.84)
#' )
#' image_size <- define_image_size(bbox, 600)
#' overlay_file <- get_arcgis_map_image(bbox, width = image_size$width,
#'                                      height = image_size$height)
#' 
get_arcgis_map_image <- function(bbox, map_type = "World_Street_Map", file = NULL, 
                                 width = 400, height = 400, sr_bbox = 4326) {
  require(httr)
  require(glue) 
  require(jsonlite)
  
  url <- parse_url("https://utility.arcgisonline.com/arcgis/rest/services/Utilities/PrintingTools/GPServer/Export%20Web%20Map%20Task/execute")
  
  # define JSON query parameter
  web_map_param <- list(
    baseMap = list(
      baseMapLayers = list(
        list(url = jsonlite::unbox(glue("https://services.arcgisonline.com/ArcGIS/rest/services/{map_type}/MapServer",
                                        map_type = map_type)))
      )
    ),
    exportOptions = list(
      outputSize = c(width, height)
    ),
    mapOptions = list(
      extent = list(
        spatialReference = list(wkid = jsonlite::unbox(sr_bbox)),
        xmax = jsonlite::unbox(max(bbox$p1$long, bbox$p2$long)),
        xmin = jsonlite::unbox(min(bbox$p1$long, bbox$p2$long)),
        ymax = jsonlite::unbox(max(bbox$p1$lat, bbox$p2$lat)),
        ymin = jsonlite::unbox(min(bbox$p1$lat, bbox$p2$lat))
      )
    )
  )
  
  res <- GET(
    url, 
    query = list(
      f = "json",
      Format = "PNG32",
      Layout_Template = "MAP_ONLY",
      Web_Map_as_JSON = jsonlite::toJSON(web_map_param))
  )
  
  if (status_code(res) == 200) {
    body <- content(res, type = "application/json")
    # message(jsonlite::toJSON(body, auto_unbox = TRUE, pretty = TRUE))
    if (is.null(file)) 
      file <- tempfile("overlay_img", fileext = ".png")
    
    img_res <- GET(body$results[[1]]$value$url)
    img_bin <- content(img_res, "raw")
    writeBin(img_bin, file)
    # message(paste("image saved to file:", file))
  } else {
    message(res)
  }
  invisible(file)
}

#' Translate the given long/lat coordinates into an image position (x, y).
#'
#' @param long longitude value
#' @param lat latitude value
#' @param bbox bounding box coordinates (list of 2 points with long/lat values)
#' @param image_width image width, in pixels
#' @param image_height image height, in pixels
#'
#' @return named list with elements "x" and "y" defining an image position
#'
find_image_coordinates <- function(long, lat, bbox, image_width, image_height) {
  x_img <- round(image_width * (long - min(bbox$p1$long, bbox$p2$long)) / abs(bbox$p1$long - bbox$p2$long))
  y_img <- round(image_height * (lat - min(bbox$p1$lat, bbox$p2$lat)) / abs(bbox$p1$lat - bbox$p2$lat))
  list(x = x_img, y = y_img)
}

#' Define image size variables from the given bounding box coordinates.
#'
#' @param bbox bounding box coordinates (list of 2 points with long/lat values)
#' @param major_dim major image dimension, in pixels. 
#'                  Default is 400 (meaning larger dimension will be 400 pixels)
#'
#' @return list with items "width", "height", and "size" (string of format "<width>,<height>")
#'
#' @examples
#' bbox <- list(
#'   p1 = list(long = -122.522, lat = 37.707),
#'   p2 = list(long = -122.354, lat = 37.84)
#' )
#' image_size <- define_image_size(bbox, 600)
#' 
define_image_size <- function(bbox, major_dim = 400) {
  # calculate aspect ration (width/height) from lat/long bounding box
  aspect_ratio <- abs((bbox$p1$long - bbox$p2$long) / (bbox$p1$lat - bbox$p2$lat))
  # define dimensions
  img_width <- ifelse(aspect_ratio > 1, major_dim, major_dim*aspect_ratio) %>% round()
  img_height <- ifelse(aspect_ratio < 1, major_dim, major_dim/aspect_ratio) %>% round()
  size_str <- paste(img_width, img_height, sep = ",")
  list(height = img_height, width = img_width, size = size_str)
}


#' Download USGS elevation data from the ArcGIS REST API.
#'
#' @param bbox bounding box coordinates (list of 2 points with long/lat values)
#' @param size image size as a string with format "<width>,<height>"
#' @param file file path to save to. Default is NULL, which will create a temp file.
#' @param sr_bbox Spatial Reference code for bounding box
#' @param sr_image Spatial Reference code for elevation image
#' 
#' @details This function uses the ArcGIS REST API, specifically the 
#' exportImage task. You can find links below to a web UI for this
#' rest endpoint and API documentation.
#' 
#' Web UI: https://elevation.nationalmap.gov/arcgis/rest/services/3DEPElevation/ImageServer/exportImage
#' API docs: https://developers.arcgis.com/rest/services-reference/export-image.htm
#'
#' @return file path for downloaded elevation .tif file. This can be read with
#' \code{read_elevation_file()}.
#'
#' @examples
#' bbox <- list(
#'   p1 = list(long = -122.522, lat = 37.707),
#'   p2 = list(long = -122.354, lat = 37.84)
#' )
#' image_size <- define_image_size(bbox, 600)
#' elev_file <- get_usgs_elevation_data(bbox, size = image_size$size)
#' 
get_usgs_elevation_data <- function(bbox, size = "400,400", file = NULL, 
                                    sr_bbox = 4326, sr_image = 4326) {
  require(httr)
  
  # TODO - validate inputs
  
  url <- parse_url("https://elevation.nationalmap.gov/arcgis/rest/services/3DEPElevation/ImageServer/exportImage")
  res <- GET(
    url, 
    query = list(
      bbox = paste(bbox$p1$long, bbox$p1$lat, bbox$p2$long, bbox$p2$lat,
                   sep = ","),
      bboxSR = sr_bbox,
      imageSR = sr_image,
      size = size,
      format = "tiff",
      pixelType = "F32",
      noDataInterpretation = "esriNoDataMatchAny",
      interpolation = "+RSP_BilinearInterpolation",
      f = "json"
    )
  )
  
  if (status_code(res) == 200) {
    body <- content(res, type = "application/json")
    # TODO - check that bbox values are correct
    # message(jsonlite::toJSON(body, auto_unbox = TRUE, pretty = TRUE))
    
    img_res <- GET(body$href)
    img_bin <- content(img_res, "raw")
    if (is.null(file)) 
      file <- tempfile("elev_matrix", fileext = ".tif")
    writeBin(img_bin, file)
    # message(paste("image saved to file:", file))
  } else {
    warning(res)
  }
  invisible(file)
}