3D plot of ranked daily streamflow using {rayshader}

rayshader-streamflow.R

library(tidyverse)
library(dataRetrieval)
library(rayshader)

# install {av} to render movie
# install.packages("av")

# fetch daily streamflow from USGS NWIS
df_nwis <- readNWISdv(siteNumbers = "01059000", parameterCd = "00060") %>% 
  renameNWISColumns() %>% 
  addWaterYear()

# compute ranks by water year
df_rank <- df_nwis %>% 
  filter(waterYear > min(waterYear), waterYear < max(waterYear)) %>%  # exclude first/last water years, which are incomplete
  group_by(waterYear) %>% 
  mutate(flowRank = row_number(Flow) - 1) %>%  # daily flow ranks range from 0 to 364
  filter(flowRank <= 365 - 1) # exclude rank 366 (leap years)

# plot
p <- df_rank  %>% 
  ggplot(aes(waterYear, flowRank)) +
  geom_raster(aes(fill = Flow)) +
  scale_fill_viridis_c("Flow (cfs)", trans = "log10", breaks = c(1e3, 1e4, 1e5), labels = scales::comma) +
  scale_x_continuous("Water Year (Oct-Sep)", expand = c(0, 0), breaks = scales::pretty_breaks(n = 10)) +
  scale_y_continuous("Daily Flow Rank\n(sorted within each water year)", expand = c(0, 0), breaks = scales::pretty_breaks(n = 10)) +
  labs(
    title = "Changes in Ranked Daily Flows by Water Year",
    subtitle = "USGS 01059000 Androscoggin River near Auburn, Maine",
    caption = "Data Source: USGS NWIS | Created by: @walkerenvres"
  )
ggsave("2d-plot.png", plot = p, width = 6, height = 5) # save ggplot

# generate rayshader plot
plot_gg(p, multicore = TRUE, width = 6, height = 5, scale = 200)
rgl::par3d(windowRect = c(0,0,1500,1500)) # increase resolution

# save static 3d
render_snapshot(filename = "3d-plot") 

# movie: tilt down
thetas_1 <- rep(0, times = 90)
phis_1 <- seq(90, 30, length.out = 90)

# movie: rotate around
thetas_2 <- seq(0, 360, length.out = 360)
phis_2 <- rep(30, times = 360)

# movie: tilt up
thetas_3 <- thetas_1
phis_3 <- rev(phis_1)

# movie: combine
thetas <- c(thetas_1, thetas_2, thetas_3)
phis <- c(phis_1, phis_2, phis_3)

# render
render_movie("movie.mp4", frames = length(thetas), fps = 30, zoom = 0.8, fov = 10, type = "custom", theta = thetas, phi = phis)