A 3D Lake District in R and Rayshader

If you follow me on Twitter, you’ll have seen that I’ve got quite excited about a new R package called Rayshader. It lets you do similar things to my¬†terrain shading efforts in Google Earth, but with much more control and cooler looking output.

In the UK, we’re lucky to have the OS terrain 50 dataset, which is 50m resolution elevation data for the whole of the mainland, available for free.

If you want to follow along with the code below, download the OS 50 data and then unzip the following “OST50GRID” zip files into a folder. You’ll also see “OST50CONT” files in the archive, but these are shapefiles for drawing contour lines and aren’t what you need for this project.

The following files will give you a map of the English Lake District:
NY10, NY11, NY12, NY20, NY21, NY22, NY30, NY31, NY32, NY40, NY41, NY42, SD19, SD29, SD39, SD49

The code snippets below are reusable for any UK location and you can look up which files you’ll need to map a particular area here.

Let’s draw a map of the Lake District, using @tylermorganwall‘s example code.



library(tidyverse)
library(rayshader)

#------------- Set Up Data ------------------------------------

#Set input directory
os50_height_dir <- "PATH_TO_UNZIPPED_OS50_FILES/"

#Load all terrain files in input directory
raster_layers <- tibble(filename = list.files(os50_height_dir, "*.asc$")) %>% 
  mutate(raster =
           map(filename, .f = ~raster::raster(rgdal::readGDAL(paste0(os50_height_dir, .))))
  ) %>% 
  pull(raster)

#Combine raster layers
raster_layers$fun <- mean
raster_mosaic <- do.call(raster::mosaic, raster_layers)

#------------- Draw Scene --------------------------------------

elmat = matrix( raster::extract(raster_mosaic, raster::extent(raster_mosaic), buffer = 1000), nrow = ncol(raster_mosaic), ncol = nrow(raster_mosaic) )

elmat %>%
  sphere_shade(zscale = 50, sunangle = 270, texture = "imhof4") %>%
  add_water(detect_water(elmat, min_area = 50), color="imhof4") %>%
  add_shadow(ray_shade(elmat, anglebreaks = seq(60, 90), sunangle = 270, multicore = TRUE, lambert = FALSE)) %>%
  add_shadow(ambient_shade(elmat, multicore = TRUE)) %>% 
  plot_3d(elmat, zscale = 50, soliddepth = -75)

 

Very pretty!

The reason for writing this post is that I wanted to do something a little bit different from the example code and overlay an image onto the terrain before ray shading it. You could use a satellite photo, but what if you want to colour the terrain by height and then shade it? What if you want to draw green valleys and pale hilltops?

I tried and I couldn’t get it to work, but then Tyler – the package’s creator – helpfully pointed me in the right direction.

Assuming you’ve already run the code above to load the elevation files and create a raster, this block will use that same raster to create a png image as the base for the rayshaded surface.



#-------------- Create an elevation shaded png -----------------------------

#Create a palette
col_ramp <- colorRampPalette(c("#54843f", "grey", "white"))

#Set up image output
png("elevation_shading.png", width=ncol(raster_mosaic), height=nrow(raster_mosaic), units = "px", pointsize = 1)

par(mar = c(0,0,0,0), xaxs = "i", yaxs = "i") #Parameters to create a borderless image

raster::image(
  raster_mosaic,
  col = col_ramp(12),
  maxpixels = raster::ncell(raster_mosaic),
  axes = FALSE
)

dev.off()

#Load generated png image from disk
terrain_image <- png::readPNG("elevation_shading.png")

#----------------- Draw Scene -----------------------------------------------

terrain_image %>% 
  add_water(detect_water(elmat, min_area = 50, remove_edges = FALSE), color="imhof4") %>%
  add_shadow(ray_shade(elmat, anglebreaks = seq(60, 90), sunangle = 270, multicore = TRUE, lambert = FALSE, remove_edges = FALSE)) %>%
  add_shadow(ambient_shade(elmat, multicore = TRUE, remove_edges = FALSE)) %>% 
  plot_3d(elmat, zscale = 50, soliddepth = -75)

And that’s the effect we were looking for!

I hope this was helpful. You could use it for all kinds of colour schemes and maybe even animate a summer to winter transition with green hilltops replaced by snowy peaks. I’m going to try to work out how to combine multiple images because I’ve got a feeling Rayshader’s sphere_shade() function might look good overlaid onto an elevation colour scheme like a multiplicative Photoshop layer.

There’ll be a new blog post if that works…