cape_town.R

#### All thanks to the wonderful Tyler Morgan-Wall for creating the fantastic Rayshader package ####
devtools::install_github("tylermorganwall/rayshader")
install.packages("needs")
needs(tidyverse, rayshader, magrittr, raster)

# Cape Town DEM downloaded from the City of Cape Town Open Data Portal: http://web1.capetown.gov.za/web1/OpenDataPortal/DatasetDetail?DatasetName=Digital%20elevation%20model
unzip("~/Downloads/10m_Grid_GeoTiff.zip")
CT <- raster("10m_Grid_GeoTiff/10m_BA.tif")
# plot the raster to make sure everything looks okay, and to get an idea of what area it covers
plot(CT)

# print the raster dimensions, to serve as a guide when cropping to our area of interest
CT@extent
# crop to our area of interest: the waterfront & table mountain
crop_CT <- crop(CT, extent(-61180, -47000, -3764020, -3751020))
# aggregate the raster by a factor of two (shrink the dimensions to make calculations and rendering faster, and given the very fine granularity of the DEM we still end up with very detailed data)
crop_CT <- aggregate(crop_CT, fact=2, fun='mean')
# plot the result to check that everything looks good
plot(crop_CT)
CTdem <- crop_CT

# convert our cropped DEM to a raster, ready for Rayshader to work its magic
CTdem = matrix(raster::extract(crop_CT,raster::extent(crop_CT),buffer=1000),
               nrow=ncol(crop_CT),ncol=nrow(crop_CT))

# replace any NA values with -1 (in some cases I have found water gets treated as NA, which can do weird things in the rendering stage. By converting NAs to -1s, we’re making sure anything below sea level is treated as water, rather than a void)
CTdem <- CTdem %>% replace(., is.na(.), -1)

# calculate shadows
raymatCT = ray_shade(CTdem)
ambmatCT = ambient_shade(CTdem)

# plot a top-down view of the resulting shaded relief map, complete with water and shadows
CTdem %>%
  sphere_shade(texture = create_texture("#d2d0b7","#68624a","#86847e","#585141","#b0b77b")) %>%
  add_water(detect_water(CTdem), color="#105b70") %>% 
  add_shadow(raymatCT) %>%
  add_shadow(ambmatCT) %>%
  plot_map()

# print the numbers of rows and columns in our DEM matrix, to serve as a reference for positioning labels in our 3D plot
ncol(CTdem)
nrow(CTdem)

# render a 3D view of our chosen scene, including shadows, a water layer and labels for key points (the render_label() and render_snapshot() commands must be run separately after the image is first rendered). Voila!
CTdem %>%
  sphere_shade(texture = create_texture("#d2d0b7","#68624a","#86847e","#585141","#b0b77b")) %>%
  add_water(detect_water(CTdem), color="#105b70") %>% 
  add_shadow(raymatCT) %>%
  add_shadow(ambmatCT) %>%
  plot_3d(CTdem, zscale=20, fov=0, theta=165, phi=4, windowsize=c(1000,800), zoom=0.35,
          water=TRUE, waterdepth=0, wateralpha=0.3,watercolor="lightblue",
          waterlinecolor="white",waterlinealpha=0.5)
render_label(CTdem,x=235,y=370, z=100,zscale=20,
             text = "Lion’s Head", textsize = 1, linewidth = 2, antialias = T)
render_label(CTdem,x=335,y=200, z=20,zscale=20,
             text = "Table Mountain", textsize = 1, linewidth = 0, antialias = T)
render_label(CTdem,x=473,y=265, z=400,zscale=20,
             text = "Devil’s Peak", textsize = 1, linewidth = 2, antialias = T)
render_snapshot()