Generates a 3d animation of changing shaded areas using lidar-derived DSM

shade_animation.r

library(raster)
library(tidyverse)
library(lubridate)
library(insol)
library(av)

#input data is a digital surface model (DSM) generated with airborne lidar
dsm1 <- raster("dsm.tif")

#convert to martix (required by rayshader)
dsm_mat = raster_to_matrix(dsm1)

#generate a sequence of sun positions for a selected day and time interval. 
timeseq <- seq(ISOdate(2020,6,22,0),ISOdate(2020,6,22,23),by="20 min")  #June 22, 2020. Change to any day here
jd = JD(timeseq)

## sun position
sp = sunpos(sunvector(jd,49.2827,-123.1207,-7)) #define location and time zone here
sp <- sp %>% as.data.frame() %>% mutate(zenith2 = 90 - zenith) %>% mutate(timeseq=timeseq) %>% filter(zenith2 >0) 

n_frames <- nrow(sp) #number of frames in the animation is equal to the number of sun positions. Shorten the time interval in timeseq to increase the number of frames.

# CAMERA POSITIONS - define how the camera position changes during the animation. One position per frame.
# thanks to @researchremora for this chunk of code

phivechalf <- 30 + 60 * 1/(1 + exp(seq(-7, 20, length.out = n_frames/2)/2))
phivecfull <- c(phivechalf, rev(phivechalf)) 

thetavec <- 0 + 45 * sin(seq(0, 719, length.out = n_frames) * pi/360) * -1

zoomvec <- 0.45 + 0.2 * 1/(1 + exp(seq(-5, 20, length.out = n_frames/2)))
zoomvecfull <- c(zoomvec, rev(zoomvec)) 


#animation frames are generate here:

for (i in 1:nrow(sp)) {
  print(i)
  
  #textured shadow map based on sun position 
  dsm_viz <- dsm_mat %>%
    sphere_shade(texture = "imhof1", sunangle = sp$azimuth[i]) %>%
    add_shadow(ray_shade(dsm_mat, zscale = 1, sunaltitude = sp$zenith2[i], sunangle = sp$azimuth[i]), 0.3) %>%
    add_shadow(ambient_shade(dsm_mat), 0) 

  
  dsm_viz %>%  plot_3d(dsm_mat, zscale = 1, fov = 0, theta = thetavec[i], zoom =zoomvecfull[i], phi = phivecfull[i], windowsize = c(1000, 800))
  
  render_snapshot(filename = sprintf("%i.png", i), 
                  title_text = paste0(hour(sp$timeseq[i]),":",minute(sp$timeseq[i])), title_size=20  )
  rgl::rgl.close()
}

flist <- paste0(1:n_frames,".png")

#output can be a gif animation or an mp4 movie
#animation
gifski::gifski(png_files = flist, gif_file = "animation.gif",delay=0.1, width = 500, height = 400)

#movie
av::av_encode_video(flist, framerate = 30, output = "animation.mp4")