Multidirectional hillshade - example with code

[jniedballa]

Multidirectional hillshades often look better than the default unidirectional hillshades computed by shade(), but are not natively supported. While they can be created with multiple calls to shade() and combing the output, below is a suggested modification to shade() which accepts multiple values in arguments angle and direction and automatically combines the resulting rasters to multidirectional hillshades:

shade <- function(slope, aspect, angle=45, direction=0, normalize=FALSE, filename="", overwrite=FALSE, combine = "mean", ...) {
  
  x <- c(slope[[1]], aspect[[1]])
  
  stopifnot(length(angle) == length(direction))
  
  if(length(angle) == 1) {
    
    direction <- direction[1] * pi/180
    zenith <- (90 - angle[1]) * pi/180
    
    if (normalize) {
      fun <- function(slp, asp) {
        shade <- cos(slp) * cos(zenith) + sin(slp) * sin(zenith) * cos(direction-asp)
        shade[shade < 0] <- 0
        shade * 255
      }
    } else {
      fun <- function(slp, asp) { cos(slp) * cos(zenith) + sin(slp) * sin(zenith) * cos(direction-asp) }
    }
    return(lapp(x, fun=fun, filename=filename, overwrite=overwrite, wopt=list(...)))
  }
  
  if(length(angle) > 1) {
    
    combine <- match.arg(combine, choices = c("mean", "Reduce"))
    
    hillsh_list <- lapply(1:length(angle), FUN = function(i) {
      direction_i <- direction[i] * pi/180
      zenith_i <- (90 - angle[i])*pi/180
      
      
      if (normalize) {
        fun <- function(slp, asp) {
          shade <- cos(slp) * cos(zenith_i) + sin(slp) * sin(zenith_i) * cos(direction_i-asp)
          shade[shade < 0] <- 0
          shade * 255
        }
      } else {
        fun <- function(slp, asp) { cos(slp) * cos(zenith_i) + sin(slp) * sin(zenith_i) * cos(direction_i-asp) }
      }
      lapp(x, fun=fun)
      
    })
    
    # combination 1: mean
    if(combine == "mean") hillsh_multi <- mean(rast(hillsh_list))
    
    # combination 2: Reduce
    if(combine == "Reduce") hillsh_multi <- Reduce(mean, hillsh_list)
    
    
    if(filename != "") {
      writeRaster(hillsh_multi, filename=filename, overwrite=overwrite, wopt=list(...))
    }
    return(hillsh_multi)
  }
}

It offers two options for combining the rasters. I found Reduce to look better when adding a light source almost in zenith (highlights ridgelines nicely, similar to the look of Terrain in Google Maps), and I'm sure the look can be further improved with other angle/direction settings.

Example:

#  sample data don't really do justice due to low resolution
f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)
alt <- disagg(r, 10, method="bilinear")
slope <- terrain(alt, "slope", unit="radians")
aspect <- terrain(alt, "aspect", unit="radians")

hill <- shade(slope, aspect, 40, 270)

hill_multi_mean <- shade(slope, aspect, 
               angle = c(45, 45, 45, 80),
               direction = c(225, 270, 315, 135),
               combine = "mean")

hill_multi_reduce <- shade(slope, aspect, 
                         angle = c(45, 45, 45, 80),
                         direction = c(225, 270, 315, 135),
                         combine = "Reduce")


plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4), main = "Default hillshade")
plot(hill_multi_mean, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4), main = "Multidirectional hillshade (mean)")
plot(hill_multi_reduce, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4), main = "Multidirectional hillshade (Reduce)")

Examples with 30m SRTM, (about 10x10km):

[rhijmans]

Thanks, shade is now vectorized for the "angle" and "direction" arguments. I have not added an argument to summarize the output, as that seems easy enough to do in a next step.

library(terra)
f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)
alt <- disagg(r, 10, method="bilinear")
slope <- terrain(alt, "slope", unit="radians")
aspect <- terrain(alt, "aspect", unit="radians")

h <- shade(slope, aspect, angle = c(45, 45, 45, 80), direction = c(225, 270, 315, 135))
h <- Reduce(mean, h)