An R script showing how to create a 3D night-time scene of the 'City of London' borough with {EAlidaR } and {Rayshader}

CityOfLondonRayshade.R

# An R script showing how to create a 3D night-time scene of the 'City of London' borough with {EAlidaR } and {Rayshader}

# If you don't already have the {EAlidaR package} Run:
# devtools::install_github('h-a-graham/EAlidaR')

library(EAlidaR)
library(raster)
library(sf)
library(tidyverse)
library(here)
library(rayshader)

# ------- Create Zipped Shapefile to upload and download Imagery ------------------

# Load city_of_london_sf from {EAlidaR} package (this is just an sf object of the London Borough)
# Then let's save it as a .shp file and zip it so we can upload it to the DEFRA portal

data_folder <- '/data' # declare folder to save data...
sf::st_write(city_of_london_sf, file.path(here(), data_folder,'COL.shp')) # write the sf object as an ESRI shapefile
filelist <- Sys.glob(file.path(here(), data_folder, 'COL.*'))  # Gather all files related to the .shp
zipped_shape <- zip::zipr(file.path(here(), data_folder,'COL.zip'), files = filelist)  # zip the files 
purrr::map(filelist, file.remove) # remove the old shapefile

# Now upload this file here: https://environment.data.gov.uk/DefraDataDownload/?Mode=survey
# Click 'Get Tiles' and select the RGB or Nighttime imagery (FYI the coverage of these aerial images is generally 
# limited but good in London...) and Download!


# ================== Full disclosure: for this next bit I cheated here and used QGIS ==========================
# I couldn't find an easy way to load .ecw data into R - my configuration of rgdal doesn't seem to like it. 
# So I loaded all the downloaded tiles into QGIS and merged all the required tiles. If you know how to load 
# .ecw data in R please let me know! So from here on we'll be dealing with a single file named 'CoL_NighTime.tif'


# ----------------------- Load the Night Time Orthomosaic  -----------------------------------

CoL_Night <- raster::brick('data/CoL_orthos/Night_Merge/CoL_NighTime.tif')

# plotRGB(CoL_Night) # take a look if you like

# get bounds for RGB raster which we use to extract the LiDAR
CoL_bounds <- st_as_sf(st_as_sfc(st_bbox(CoL_Night))) %>%
  st_set_crs(., 27700 )

# get 0.5m DSM data using the EAlidaR package.
CoL_dsm <- EAlidaR::get_area(poly_area = CoL_bounds , resolution=0.5, model_type = 'DSM', merge_tiles = TRUE, crop=TRUE)


# --------------- Set up Matrices for Rayshader --------------

# Function to convert RGB raster brick to a matrix readable with rayshader:
# code adapted from:https://www.tylermw.com/a-step-by-step-guide-to-making-3d-maps-with-satellite-imagery-in-r/
brick_to_matrix <- function(raster_brick){
  names(raster_brick)[1] = "r" 
  names(raster_brick)[2] = "g"
  names(raster_brick)[3] = "b"
  
  r_matrix = rayshader::raster_to_matrix(raster_brick$r)
  g_matrix = rayshader::raster_to_matrix(raster_brick$g)
  b_matrix = rayshader::raster_to_matrix(raster_brick$b)

  rgb_array = array(0,dim=c(nrow(r_matrix),ncol(r_matrix),3))
  
  rgb_array[,,1] = r_matrix/255 #Red layer
  rgb_array[,,2] = g_matrix/255 #Blue layer
  rgb_array[,,3] = b_matrix/255 #Green layer
  
  return(aperm(rgb_array, c(2,1,3)))
}

# generate matrix of RGB imagery
Col_rgb_array = brick_to_matrix(CoL_Night)

# generate matrix of Elevation data
CoL_el_matrix = raster_to_matrix(CoL_dsm)

# -------------- Use {rayshader} to render the 3D scene ----------------
# plot the 3D scene with rayshader...
plot_3d(Col_rgb_array, CoL_el_matrix, windowsize = c(1200,900), shadow=FALSE, solid = FALSE, 
        zscale = 0.6, zoom=0.4, phi=40,theta=-135,fov=70, background = "#000000")

# create a cool depth of field image (Add filename = '...png' to save to disk)
Sys.sleep(0.2)
render_depth(focus = 0.7, focallength = 70, clear = FALSE)

# create the spinning scene using default orbit settings...
render_movie(filename='CoL_NightFlight.mp4', fps = 24)