Last active
February 3, 2018 20:42
-
-
Save eliocamp/74d32690bab4532f7123ef1a18e4d625 to your computer and use it in GitHub Desktop.
Geom for relief shading.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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)) | |
# Esto es lo único que es nuevo. El resto es básicamente copy-paste | |
# de geom_raster y 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)] | |
# Desde geom_raster y 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) | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment