Simulate some geological layers, plot them in a rgl 3D scene and save as .gif

simulated_geology_3dgif.R

library(tidyverse)
library(raster)
library(magick)
library(NLMR)
library(rgl)

# Let's start by simulating a digital elevation model
sim_DEM <- NLMR::nlm_fbm(
  ncol = 100, nrow = 100, resolution = 10,
  fract_dim = 1.6,
  user_seed = 42, rescale = TRUE
) %>% 
  magrittr::add(-0.5) %>% 
  magrittr::multiply_by(10)

# The DEM is (for obvious reasons) the limit of the geological layers
# If a lithology is to cut through the surface, it is eroded
# In this example every lithology is to be eroded by the layer above.
erode <- function(upper_raster, lower_raster){
  
  # Get the elevation of both raster objects
  up <- raster::values(upper_raster)
  lo <- raster::values(lower_raster)
  
  eroded <- purrr::map2_dbl(
    up, lo, function(x, y){
      if(is.na(x)) { # If there is no upper layer, the lower one persists
        y
      } else if(is.na(y)) { # If there is no lower layer, none will be produced
        NA
      } else if(x < y){ # Upper layer below lower layer! Erosion!
        x # Values of the upper lithology are kept (though it is technically lower)
      } else {
        y # If nothing happens nothing happens
      }
    }
  )
  
  # replace with modified elevation
  raster::setValues(lower_raster, eroded)
}

# the uppermost lithology shall be a silty channel, representing 
# some kind of former river bed
sim_silt <- NLMR::nlm_edgegradient(
  ncol = 100, nrow = 100, resolution = 10, direction = 90
) %>% 
  magrittr::multiply_by(-1) %>% # reverse form
  magrittr::multiply_by(5) %>%  # shift the plane into position 
  erode(sim_DEM, .)             # must not be higher that the topographical surface

# the underlying sand gets completely cut through by the river
sim_sand <- NLMR::nlm_fbm(
  ncol = 100, nrow = 100, resolution = 10,
  fract_dim = 1.6,
  user_seed = 321, rescale = TRUE
) %>% 
  magrittr::add(-0.5) %>% 
  magrittr::multiply_by(5) %>% 
  magrittr::add(-4.5) %>% 
  erode(sim_silt, .)

sim_marl <- NLMR::nlm_distancegradient(
  ncol = 100, nrow = 100, resolution = 10, origin = c(20, 30, 10, 15)
) %>% 
  magrittr::multiply_by(10) %>% # Making the height differences a bit more dramatic
  magrittr::add(-12) %>% 
  erode(sim_sand, .)

# Finally a custom plotting function to make visualization a bit more smooth
gm_surface_3d <- function(raster, exag, color){
  z <- as.matrix(raster)
  x <- 1:nrow(z) 
  y <- 1:nrow(z)
  
  rgl::surface3d(x, y, z*exag, color)
}

# set window size for the 3d plot
r3dDefaults$windowRect <- c(100, 100, 500, 500)

# add all the layers to the scene
walk2(
  .x = list(sim_DEM, sim_silt, sim_sand, sim_marl),
  .y = list("green", "brown", "yellow", "blue"),
  .f =gm_surface_3d,
  exag = 5
)

# make a gif and save to file
movie3d(
  spin3d(), duration = 12, fps = 20,
  dir = paste0(getwd()),
  movie = "sim_geomod"
  )