geom_relief.R

## from https://eliocamp.github.io/codigo-r/2018/02/how-to-make-shaded-relief-in-r/

geom_relief <- function(mapping = NULL, data = NULL,
                        stat = "identity", position = "identity",
                        ...,
                        raster = TRUE,
                        interpolate = TRUE,
                        na.rm = FALSE,
                        show.legend = NA,
                        inherit.aes = TRUE) {
  ggplot2::layer(
    data = data,
    mapping = mapping,
    stat = stat,
    geom = GeomRelief,
    position = position,
    show.legend = show.legend,
    inherit.aes = inherit.aes,
    params = list(
      raster = raster,
      interpolate = interpolate,
      na.rm = na.rm,
      ...
    )
  )
}

GeomRelief <- ggplot2::ggproto("GeomRelief", GeomTile,
                               required_aes = c("x", "y", "z"),
                               default_aes = ggplot2::aes(color = NA, fill = "grey35", size = 0.5, linetype = 1,
                                                          alpha = NA, light = "white", dark = "gray20", sun.angle = 60),
                               draw_panel = function(data, panel_scales, coord, raster, interpolate) {
                                 if (!coord$is_linear()) {
                                   stop("non lineal coordinates are not implemented in GeomRelief", call. = FALSE)
                                 } else {
                                   coords <- as.data.table(coord$transform(data, panel_scales))
                                   
                                   # This is the only part that's actually new. The rest is essentially 
                                   # copy-pasted from geom_raster and geom_tile
                                   coords[, sun.angle := (sun.angle + 90)*pi/180]
                                   coords[, dx := .derv(z, x), by = y]
                                   coords[, dy := .derv(z, y), by = x]
                                   coords[, shade := (cos(atan2(-dy, -dx) - sun.angle) + 1)/2]
                                   coords[is.na(shade), shade := 0]
                                   coords[, fill := .rgb2hex(colorRamp(c(dark, light), space = "Lab")(shade)),
                                          by = .(dark, light)]
                                   
                                   # From geom_raster and geom_tile
                                   if (raster == TRUE){
                                     if (!inherits(coord, "CoordCartesian")) {
                                       stop("geom_raster only works with Cartesian coordinates", call. = FALSE)
                                     }
                                     # Convert vector of data to raster
                                     x_pos <- as.integer((coords$x - min(coords$x)) / resolution(coords$x, FALSE))
                                     y_pos <- as.integer((coords$y - min(coords$y)) / resolution(coords$y, FALSE))
                                     
                                     nrow <- max(y_pos) + 1
                                     ncol <- max(x_pos) + 1
                                     
                                     raster <- matrix(NA_character_, nrow = nrow, ncol = ncol)
                                     raster[cbind(nrow - y_pos, x_pos + 1)] <- alpha(coords$fill, coords$alpha)
                                     
                                     # Figure out dimensions of raster on plot
                                     x_rng <- c(min(coords$xmin, na.rm = TRUE), max(coords$xmax, na.rm = TRUE))
                                     y_rng <- c(min(coords$ymin, na.rm = TRUE), max(coords$ymax, na.rm = TRUE))
                                     
                                     grid::rasterGrob(raster,
                                                      x = mean(x_rng), y = mean(y_rng),
                                                      width = diff(x_rng), height = diff(y_rng),
                                                      default.units = "native", interpolate = interpolate
                                     )
                                     
                                   } else {
                                     ggplot2:::ggname("geom_rect", grid::rectGrob(
                                       coords$xmin, coords$ymax,
                                       width = coords$xmax - coords$xmin,
                                       height = coords$ymax - coords$ymin,
                                       default.units = "native",
                                       just = c("left", "top"),
                                       gp = grid::gpar(
                                         col = coords$fill,
                                         fill = alpha(coords$fill, coords$alpha),
                                         lwd = coords$size * .pt,
                                         lty = coords$linetype,
                                         lineend = "butt"
                                       )
                                     ))
                                     
                                   }
                                 }
                               }
)

rect_to_poly <- function(xmin, xmax, ymin, ymax) {
  data.frame(
    y = c(ymax, ymax, ymin, ymin, ymax),
    x = c(xmin, xmax, xmax, xmin, xmin)
  )
}

.rgb2hex <- function(array) {
  rgb(array[, 1], array[, 2], array[, 3], maxColorValue = 255)
}


.derv <- function(x, y, order = 1, cyclical = FALSE, fill = FALSE) {
  N <- length(x)
  d <- y[2] - y[1]
  if (order >= 3) {
    dxdy <- .derv(.derv(x, y, order = 2, cyclical = cyclical, fill = fill),
                  y, order = order - 2, cyclical = cyclical, fill = fill)
  } else {
    if (order == 1) {
      dxdy <- (x[c(2:N, 1)] - x[c(N, 1:(N-1))])/(2*d)
    } else if (order == 2) {
      dxdy <- (x[c(2:N, 1)] + x[c(N, 1:(N-1))] - 2*x)/d^2
    }
    if (!cyclical) {
      if (!fill) {
        dxdy[c(1, N)] <- NA
      }
      if (fill) {
        dxdy[1] <- (-11/6*x[1] + 3*x[2] - 3/2*x[3] + 1/3*x[4])/d
        dxdy[N] <- (11/6*x[N] - 3*x[N-1] + 3/2*x[N-2] - 1/3*x[N-3])/d
      }
    }
    
  }
  return(dxdy)
}