Relief Shade in R

relief-shade.R

library(data.table)
library(ggplot2)
data(volcano)
volcano <- as.data.table(melt(volcano, varnames = c("x", "y"),
                              value.name = "h"))


volcano[, c("dx", "dy") := metR::Derivate(h ~ x + y)]  # from my package. It calculates directional derivatives. 
volcano[, angle := atan2(-dy, -dx)]

sun.angle <- pi/3
ggplot(volcano, aes(x, y)) +
   geom_raster(aes(fill = cos(angle + sun.angle)), alpha = 1, interpolate = TRUE) +
   scale_fill_gradient2(low = "white", high = "white", mid = "gray20", 
                        midpoint = sun.angle, guide = "none") +
   coord_fixed() +
   theme_void()