geom_relief_rayshader.R

geom_relief <- function(mapping = NULL, data = NULL,
                        stat = "identity", position = "identity",
                        ...,
                        anglebreaks = 60,
                        sunangle = 60,
                        shadow = TRUE,
                        maxsearch = 100,
                        lambert = TRUE,
                        zscale = 1,
                        multicore = TRUE,
                        remove_edges = TRUE,
                        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(
        anglebreaks = anglebreaks,
        sunangle = sunangle,
        shadow = shadow,
        maxsearch = maxsearch,
        lambert = lambert,
        zscale = zscale,
        multicore = multicore,
        remove_edges = remove_edges,
        raster = raster,
        interpolate = interpolate,
        na.rm = na.rm,
        ...
      )
    )
}


GeomRelief <- ggplot2::ggproto("GeomRelief", ggplot2::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"),
  draw_panel = function(data, panel_scales, coord, raster, interpolate, shadow,
                        anglebreaks, sunangle, maxsearch, lambert, zscale, 
                        multicore, remove_edges) {
      `%>%` <- dplyr::`%>%`
  
      coords <- coord$transform(data, panel_scales)
      heightmap <- as.matrix(reshape2::dcast(coords, x ~ y, value.var = "z")[, -1])
      
      if (shadow == TRUE) {
      coords$z <- c(rayshader::rayshade(heightmap = heightmap, 
                                        anglebreaks = anglebreaks, 
                                        sunangle = sunangle, 
                                        maxsearch = maxsearch,
                                        lambert = lambert, 
                                        zscale = zscale, 
                                        multicore = multicore,
                                        remove_edges = FALSE))
      } else {
        coords$z <- c(rayshader::lambshade(heightmap = heightmap, 
                                          rayangle = mean(anglebreaks), 
                                          sunangle = sunangle,
                                          zscale = zscale, 
                                          zero_negative = TRUE,
                                          remove_edges = FALSE))
      }
      
      if (remove_edges == TRUE) {
        xs <- unique(coords$x)
        ys <- unique(coords$y)
        coords <- subset(coords, x %in% xs[3:(length(xs)-2)] &
                                 y %in% ys[3:(length(ys)-2)])
      }
      
      coords <- coords %>%
        dplyr::group_by(light, dark) %>%
        dplyr::mutate(fill = scales::colour_ramp(c(dark[1], light[1]))(z)) %>%
        dplyr::ungroup()
      
      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"
          )
        ))
        
      }
  }
)


v <- reshape2::melt(volcano)
ggplot(v, aes(Var1, Var2)) +
  geom_relief(aes(z = value))