In this example, we will learn about 3-D visualisation in R using the package rayshader. The code used in this tutorial was modified from https://www.rayshader.com/. First, we will create a digital elevation model (DEM) with a 1-m spatial resolution. Further, we will add an aerial image at a 20-cm spatial resolution above the DEM. Finally, we will make an animation with 3-D visualisation. The study area is a Rummu quarry (59.227398, 24.193968) - a man-made heap that originated due to the mining of Vasalemma limestone and marble. This place is very popular during the hot season since it is possible to swim and dive in transparent water (link).

Lat’s start with packages installation. We have to install packages from GitHub (not from CRAN! ), otherwise, you will face the problem with plotting. Package remotes is designed to enable downloading the packages from GitHub repositories.

install.packages("remotes")
library(remotes)
remotes::install_github("tylermorganwall/rayshader") 
remotes::install_github("dmurdoch/rgl")

install.packages("magick") # package implements functionality for manipulating 
                            # high-dimensional arrays using efficient vectorised methods

Call the packages we will use:

library(magick)
library(rayshader)
library(raster)
library(rayrender)

Download the DEM data from Landboard Estonia. link

Select 1-m resolution GeoTIFF data for map sheet number 63613. After the file is downloaded - place it to your working directory.

Set your working directory. This should be the folder where you can save the outputs.

setwd("/YOUR/WORKING/DIRECTORY")

Import downloaded raster:

DEM  <-  raster::raster( "63613_dem_1m.tif")   # It works if file "63613_dem_1m.tif" is located in your working directory
crs(DEM) # see the coordinate system
## Coordinate Reference System:
## Deprecated Proj.4 representation:
##  +proj=lcc +lat_0=57.5175539305556 +lon_0=24 +lat_1=59.3333333333333
## +lat_2=58 +x_0=500000 +y_0=6375000 +ellps=GRS80 +units=m +no_defs 
## WKT2 2019 representation:
## PROJCRS["Estonian Coordinate System of 1997 (with axis order normalized for visualization)",
##     BASEGEOGCRS["EST97",
##         DATUM["Estonia 1997",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4180]],
##     CONVERSION["Estonian National Grid",
##         METHOD["Lambert Conic Conformal (2SP)",
##             ID["EPSG",9802]],
##         PARAMETER["Latitude of false origin",57.5175539305556,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8821]],
##         PARAMETER["Longitude of false origin",24,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8822]],
##         PARAMETER["Latitude of 1st standard parallel",59.3333333333333,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8823]],
##         PARAMETER["Latitude of 2nd standard parallel",58,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8824]],
##         PARAMETER["Easting at false origin",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8826]],
##         PARAMETER["Northing at false origin",6375000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8827]]],
##     CS[Cartesian,2],
##         AXIS["easting (Y)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["northing (X)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]]]

Plot the DEM file

plot(DEM,col = terrain.colors(n=50))

DEM covers a big area while we need a small part of it. Rummu quarry is the small area in the bottom left corner of the figure. Create the polygon to clip DEM to the area of interest:

polygon <- as(extent(510800, 511300, 6565300, 6565750), 'SpatialPolygons')

crs(polygon) <- crs(DEM) # assign projection to the polygon

crs(polygon) # polygon and DEM have the same coordinate system
## Coordinate Reference System:
## Deprecated Proj.4 representation:
##  +proj=lcc +lat_0=57.5175539305556 +lon_0=24 +lat_1=59.3333333333333
## +lat_2=58 +x_0=500000 +y_0=6375000 +ellps=GRS80 +units=m +no_defs 
## WKT2 2019 representation:
## PROJCRS["Estonian Coordinate System of 1997 (with axis order normalized for visualization)",
##     BASEGEOGCRS["EST97",
##         DATUM["Estonia 1997",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4180]],
##     CONVERSION["Estonian National Grid",
##         METHOD["Lambert Conic Conformal (2SP)",
##             ID["EPSG",9802]],
##         PARAMETER["Latitude of false origin",57.5175539305556,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8821]],
##         PARAMETER["Longitude of false origin",24,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8822]],
##         PARAMETER["Latitude of 1st standard parallel",59.3333333333333,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8823]],
##         PARAMETER["Latitude of 2nd standard parallel",58,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8824]],
##         PARAMETER["Easting at false origin",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8826]],
##         PARAMETER["Northing at false origin",6375000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8827]]],
##     CS[Cartesian,2],
##         AXIS["easting (Y)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["northing (X)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]]]
DEM_crop <- crop(DEM, polygon) # crop DEM by polygon

plot(DEM_crop,col = terrain.colors(n=50))

To use rayshader package, we need to convert DEM from raster to a matrix

DEM_matrix  <-  raster_to_matrix(DEM_crop)

Let’s try the first image:

DEM_matrix %>%
  sphere_shade(texture = "desert") %>% # `imhof1`,`imhof2`,`imhof3`,`imhof4`,`desert`, `bw`
  plot_map()

You can set more parameters (see more: https://www.rayshader.com/reference/sphere_shade.html)
Try to change the ‘sunangle’ and ‘texture’. For example:

DEM_matrix %>%
  sphere_shade(sunangle = 200, texture = "unicorn") %>%
  plot_map()

Create your own texture and increase the color intensity (default is 1):

DEM_matrix %>%
  sphere_shade(texture=create_texture("#FFA900",
                                      "#FF7600",
                                      "#A2CDCD",
                                      "#CD113B",
                                      "black"),
               colorintensity = 5) %>%
  plot_map()

plot_3d() function in rayshader allows plotting in 3-D:

DEM_matrix %>%
  sphere_shade(texture = "unicorn", 
               sunangle = 10) %>%
  plot_3d(DEM_matrix, 
         baseshape = "hex", # Shape of the base: "rectangle","circle","hex"
          solidcolor = 'black', # Base color
          solidlinecolor = '#191A19', #Base edge line color
          background = '#1F1D36',# Color of the background
          zscale = 0.4, # ratio between the x and y spacing
          fov = 0, # field-of-view angle
          theta = 60, # Rotation around z-axis
          zoom = 0.7, # Zoom factor
          phi = 45, # Azimuth angle
          windowsize = c(500, 500))

This figure is plotted in the new window where you can rotate DEM. Find the most beautiful perspective and run the following script:

Sys.sleep(0.2)
render_snapshot()

After your figure appeared in the Plot window - you can download it (if needed).

You can get a spectacular visualisation with render_highquality() function. It is process-demanding, so wait for 1-3 minutes:

render_highquality(lightdirection = 100, #position of the light angle around the scene
                   lightaltitude = 45, #Angle above the horizon that the light is located.
                   lightintensity = 400, # Intensity of the light.
                   clamp_value = 10, 
                   title_text = "Rummu quarry",
                   title_color = "white",
                   title_bar_color = "black", scene_elements = NULL,
                   camera_location = NULL, camera_lookat = NULL,
                   parallel = TRUE, width = 1000, height = 1000, 
                   samples = 1000) # depending on your computer use from 200 to 2000

rgl::rgl.close() #Close window with image

Rummu heap is located nearby the water body. We can show this water object, assuming that everything that is below 23 m is flooded.

DEM_matrix %>%
  sphere_shade(texture = "imhof1",  sunangle = 100) %>%
  plot_3d(DEM_matrix, 
          baseshape = "circle", # Shape of the base: "rectangle","circle","hex"
          solidcolor = 'grey20', # Base color
          solidlinecolor = 'grey20', #Base edge line color
          zscale = 0.4, # ratio between the x and y spacing
          fov = 0, # field-of-view angle
          theta = 60, # Rotation around z-axis
          zoom = 0.7, # Zoom factor
          phi = 45, # Azimuth angle
          windowsize = c(500, 500),
          water = TRUE, 
          waterdepth = 23, # test 50 and observe some catastrophic scenario
          wateralpha = 0.5, watercolor = "#233aa1",
          waterlinecolor = "white", waterlinealpha = 0.5) 

Sys.sleep(0.2)

render_snapshot()

rgl::rgl.close()

Now the DEM is ready and we can start with an aerial image.
Download the image from Landboard Estonia link.
Insert the same number of the sheet as you used to download DEM. Place the image in your working directory folder:

image  <-  stack("63613.tif") # load RGB image

plotRGB(image)

The area of the image is too big. We can crop the image using the polygon object we created before.

image_crop <- raster::crop(image, polygon)

plotRGB(image_crop)

Convert raster to the matrix. Each band (Red, Green, Blue) should be converted separately.

names(image_crop) <-  c("r","g","b")

image_crop_R <-  rayshader::raster_to_matrix(image_crop$r) # Red band
image_crop_G <-  rayshader::raster_to_matrix(image_crop$g) # Green band
image_crop_B <-  rayshader::raster_to_matrix(image_crop$b) # Blue band

image_crop_RGB <-  array(0,dim=c(nrow(image_crop_R),ncol(image_crop_R),3))

image_crop_RGB[,,1] <-  image_crop_R #Red layer
image_crop_RGB[,,2] <-  image_crop_G #Blue layer
image_crop_RGB[,,3] <-  image_crop_B #Green layer

image_crop_RGB <-  aperm(image_crop_RGB, c(2,1,3))

image_crop_RGB <-  scales::rescale(image_crop_RGB,to=c(0,1)) # rescale image

plot_map(image_crop_RGB)

Now we can plot RGB image above the DEM:

plot_3d(image_crop_RGB, # RGB image
        DEM_matrix, # DEM 
        windowsize = c(700,600),
        zscale = 0.6, 
        shadowdepth = -50,
        zoom=0.7, phi=45,theta=-45,fov=70, 
        background = "#D5D5D5", 
        shadowcolor = "#5e5e5e")

You can add labels to the plot (read more https://www.rayshader.com/reference/render_label.html). For example, add one viewing point and diving place:

render_label(DEM_matrix,x = 180, y = 240,
             altitude=30, zscale=0.6, 
             text = "Viewing point",
             textcolor = "black",
             linecolor="black",
             dashed = TRUE,
             alpha = 0.9,
             textsize = 1.5)

render_label(DEM_matrix,x = 390, y = 370,
             altitude=60, zscale=0.6, 
             text = "Diving place",
             textcolor = "black",
             linecolor="black",
             dashed = TRUE,
             alpha = 0.9,
             textsize = 1.5)

render_snapshot(title_text = "Rummu quarry \nData: Land Board ",
                title_color = "white", 
                title_bar_alpha = 1)

After you select the nice view, we can make a video around the heap using av package.

install.packages('av')
library(av)

First, we will create 360 images for each degree of view (overall 360 degrees):

angles <- seq(0,360,length.out = 361)[-1] # one image for one degree (360 in total)

for(i in 1:360) {
  render_camera(theta=-45+angles[i])
  render_snapshot(filename = sprintf("Rummu_quarry%i.png", i), 
                  title_text = "Rummu quarry \nData: Land Board",
                  title_color = "white", 
                  title_bar_alpha = 1)
}

rgl::rgl.close()

Second, we can create a video with those images:

av::av_encode_video(sprintf("Rummu_quarry%d.png",seq(1,360,by=1)), 
                    framerate = 30, output = "Rummu_quarry.mp4")

Find the video in your working directory (getwd()).

Is it similar?


In the future, you can use rayshader with ggplot2 to make your images in 3-D. Read more information here: link and try to make some 3-D plots in the next lesson.

Author: Iuliia Burdun, Anto Aasa
Supervisors: Anto Aasa & Lika Zhvania
LTOM.02.041
Last update: 2022-11-25 09:05:55