I’ve spent plenty of time recently with Tyler Morgan Wall‘s rayshader R package, learning a lot about visualising geographic data along the way.
I fly paragliders and early on in my exploration, put together a few helper functions to add gps traces from glider flights to rayshader scenes and to make it easier to prepare elevation data. Rayshader itself is fabulous but battling with government portals to find open elevation data and then working out how to stitch it together is less fun, so those early functions have gradually morphed into my first R CRAN package – geoviz.
What’s geoviz? You might already know ggmap, that makes it dead easy to draw maps in ggplot2. Well I wanted to make a ggmap for rayshader; a set of functions that could take you in a couple of lines of code from lat-long coords or a gps trace file, to a rayshader 3D scene, without messing about downloading elevation files and working out how to process them.
A couple of acknowledgements up front… geoviz only exists because rayshader exists. As I say in the git repo, “Rayshader is an awesome bit of kit! I’m just doing some colouring in.” I’m also indebted to Miles McBain’s slippymath, which geoviz leans on to download map tiles from Mapbox, Mapzen and Stamen, and to Michael Sumner for bits and pieces of advice and for geoviz’s first external contributor pull request.
tl:dr, skip to the bottom for the complete script to draw this scene.
Right then, let’s draw… wherever you like! Get yourself some lat-long coordinates for a favourite location. Rayshader is all about 3D viz, so it helps if you pick somewhere mountainous. You can draw somewhere pancake-flat if you want, it just won’t look as interesting. It’s easy to get your coordinates by manually searching and clicking a spot on Google Maps, or you can use ggmap::geocode() but Google’s API restrictions now mean you’ll need to set up a Google API key, which is a bit laborious if you haven’t done it before.
Picked somewhere? Great. I’m going to use Ullswater in the English Lake District for this example. It’s a beautiful spot – a lake in North West England surrounded by mountains.
You’re going to need geoviz version 0.2.1 (or later), which has just made its way to CRAN.
install.packages("geoviz")
The first thing you need is a Digital Elevation Model (DEM). This is a file of terrain height data, that rayshader will use as the base for your scene. You could go and Google for files that cover the area you want to visualise, download them and stitch them together, or you could do this.
library(rayshader)
library(geoviz)
#Set up the area you want to visualise (ggmap::geocode can help you do this straight from R but you'll need a Google API key)
lat <- 54.5750474
long <- -2.8936547
square_km <- 10
#Increase this to ~60 for a higher resolution (but slower) image
max_tiles <- 10
#Get elevation data from Mapzen
dem <- mapzen_dem(lat, long, square_km, max_tiles = max_tiles)
Note that we've set max_tiles to 10 here, which will make the code run faster and put less load on the remote map server, but gives a lower quality image. Set max_tiles = 60 for a high quality scene.
You can use rayshader to visualise the DEM as it is if you like. For more on how rayshader's options work, see its own site here.
#Draw the rayshader scene
elmat = matrix(
raster::extract(dem, raster::extent(dem), method = 'bilinear'),
nrow = ncol(dem),
ncol = nrow(dem)
)
sunangle <- 270
scene <- elmat %>%
sphere_shade(sunangle = sunangle, texture = "bw")
rayshader::plot_3d(
scene,
elmat,
zscale = raster_zscale(dem),
solid = FALSE,
shadow = TRUE,
shadowdepth = -150
)
Now we want to add an overlay image. First, we can try a Stamen Watercolor map, which gives the scene a lovely cartoony look. The overlay is semi-transparent, so that shading on the terrain base shows through.
overlay_image <-
slippy_overlay(
dem,
image_source = "stamen",
image_type = "watercolor",
png_opacity = 0.3
)
#Draw the rayshader scene
scene <- elmat %>%
sphere_shade(sunangle = sunangle, texture = "bw") %>%
add_overlay(overlay_image)
rayshader::plot_3d(
scene,
elmat,
zscale = raster_zscale(dem),
solid = FALSE,
shadow = TRUE,
shadowdepth = -150
)
What about a satellite image? We can do that with geoviz too, but you'll need a mapbox API key. Register for one and then the process works just like it did for the Stamen overlay.
# If you haven't got a mapbox key, go here first: https://docs.mapbox.com/help/glossary/access-token
mapbox_key <- "YOUR MAPBOX KEY"
overlay_image <-
slippy_overlay(
dem,
image_source = "mapbox",
image_type = "satellite",
png_opacity = 0.6,
api_key = mapbox_key
)
#Draw the rayshader scene
scene <- elmat %>%
sphere_shade(sunangle = sunangle, texture = "bw") %>%
add_overlay(overlay_image)
rayshader::plot_3d(
scene,
elmat,
zscale = raster_zscale(dem),
solid = FALSE,
shadow = TRUE,
shadowdepth = -150
)
That's it! The next post will show you how to add a GPS trace to your scene and visualise runs, cycle rides, flights, animal tracking, or whatever else you can think of to strap a GPS to.
You might have noticed that your images didn't look quite as nice as the screenshots on this post. That's because - as mentioned earlier - you need to set max_tiles = 60 for a high quality scene and also because I didn't put rayshader's full range of shadows into the snippets above, since they take a while to calculate. This is also why your scenes might be rectangular at the moment instead of square - more tiles will fix it.
Here's the full script to draw the Stamen Watercolor and Mapbox Satellite scenes, in glorious high resolution with deep shadows. Enjoy!
devtools::install_github("neilcharles/geoviz")
library(rayshader)
library(geoviz)
#Set up the area you want to visualise (ggmap::geocode can help you do this straight from R but you'll need a Google API key)
lat <- 54.5750474
long <- -2.8936547
square_km <- 10
#Increase this to ~60 for a higher resolution (but slower) image
max_tiles <- 60
#Get elevation data from Mapzen
dem <- mapzen_dem(lat, long, square_km, max_tiles = max_tiles)
overlay_image <-
slippy_overlay(
dem,
image_source = "stamen",
image_type = "watercolor",
png_opacity = 0.3,
max_tiles = max_tiles
)
# Or for a satellite image, use this block
# If you haven't got a mapbox key, go here first: https://docs.mapbox.com/help/glossary/access-token
mapbox_key <- "YOUR MAPBOX KEY"
#overlay_image <-
# slippy_overlay(
# dem,
# image_source = "mapbox",
# image_type = "satellite",
# png_opacity = 0.6,
# api_key = mapbox_key
# )
sunangle <- 270
#Draw the rayshader scene
elmat = matrix(
raster::extract(dem, raster::extent(dem), method = 'bilinear'),
nrow = ncol(dem),
ncol = nrow(dem)
)
scene <- elmat %>%
sphere_shade(sunangle = sunangle, texture = "bw") %>%
add_overlay(overlay_image) %>%
#The next two lines create deep shadows but are slow to run at high quality
add_shadow(
ray_shade(
elmat,
anglebreaks = seq(30, 60),
sunangle = sunangle,
multicore = TRUE,
lambert = FALSE,
remove_edges = FALSE
)
) %>%
add_shadow(ambient_shade(elmat, multicore = TRUE, remove_edges = FALSE))
rayshader::plot_3d(
scene,
elmat,
zscale = raster_zscale(dem),
solid = FALSE,
shadow = TRUE,
shadowdepth = -150
)
rgl::view3d(theta =290, phi = 18, zoom = 0.5, fov = 5)
rayshader::render_depth(
focus = 0.5,
fstop = 18,
filename = "scene.png"
)