ggplot2: raster plotting does not work as expected when setting alpha values

Question

First post here, I hope I'm observing website etiquette. I couldn't find and answer on the site and I previously posted this to a ggplot2 specific group, but no solutions as yet.

Basically I am trying to overlay two rasters using ggplot2 and require the top one to be semi-transparent. I have a hillShade raster which is computed from an elevation data raster, and I wish to overlay the elevation raster onto the hillshade raster so the resulting plot doesn't look 'flat'. You can see what I mean in the reproducible R code below.

Using base graphics I can achieve the desired result and I have included an example in code below to make it clear what I mean, but I need to do this in ggplot2.

I can't get it to work in ggplot2. Combining the rasters makes the colours go funny (I can plot each one ok by itself). Can anyone help or point me in the right direction. Self contained, reproducible code example included below. (Sorry for the length, but I thought better to be clear).

#   Load relevant libraries
library(ggplot2)
library(raster)


#   Download sample raster data of Ghana from my Dropbox
oldwd <- getwd()
tmp <- tempdir()
setwd(tmp)
url1 <- "http://dl.dropbox.com/s/xp4xsrjn3vb5mn5/GHA_HS.asc"
url2 <- "http://dl.dropbox.com/s/gh7gzou9711n5q7/GHA_DEM.asc"
f1 <- file.path(tmp,"GHA_HS.asc")
f2 <- file.path(tmp,"GHA_DEM.asc")
download.file(url1,f1)  #File is ~ 5,655Kb
download.file(url2,f2)  #File is ~ 2,645Kb


#   Create rasters from downloaded files
hs <-  raster(f1)
dem <- raster(f2)


#   Plot with base graphics to show desired output
plot(hs,col=grey(1:100/100),legend=F)
plot(dem,col=rainbow(100),alpha=0.4,add=T,legend=F)


#   Convert rasters TO dataframes for plotting with ggplot
hdf <- rasterToPoints(hs); hdf <- data.frame(hdf)
colnames(hdf) <- c("X","Y","Hill")
ddf <- rasterToPoints(dem); ddf <- data.frame(ddf)
colnames(ddf) <- c("X","Y","DEM")


#   Create vectors for colour breaks
b.hs <- seq(min(hdf$Hill),max(hdf$Hill),length.out=100)
b.dem <- seq(min(ddf$DEM),max(ddf$DEM),length.out=100)


#   Plot DEM layer with ggplot()
p1 <- ggplot()+
    layer(geom="raster",data=ddf,mapping=aes(X,Y,fill=DEM))+
    scale_fill_gradientn(name="Altitude",colours = rainbow(100),breaks=b.dem)+
    scale_x_continuous(name=expression(paste("Longitude (",degree,")")),limits=c(-4,2),expand=c(0,0))+
    scale_y_continuous(name=expression(paste("Latitude (",degree,")")),limits=c(4,12),expand=c(0,0))+
    coord_equal()
print(p1)


#   Plot hillShade layer with ggplot()
p2 <- ggplot()+
    layer(geom="raster",data=hdf,mapping=aes(X,Y,fill=Hill))+
    scale_fill_gradientn(colours=grey(1:100/100),breaks=b.hs,guide="none")+
    scale_x_continuous(name=expression(paste("Longitude (",degree,")")),limits=c(-4,2),expand=c(0,0))+
    scale_y_continuous(name=expression(paste("Latitude (",degree,")")),limits=c(4,12),expand=c(0,0))+
    coord_equal()
print(p2)


#   Try to plot both together with transparency on the DEM layer
p3 <- ggplot(hdf)+
    geom_raster(aes(X,Y,fill=Hill))+
    scale_fill_gradientn(colours=grey(1:100/100),breaks=b.hs,guide="none")+
    scale_x_continuous(name=expression(paste("Longitude (",degree,")")),limits=c(-4,2),expand=c(0,0))+
    scale_y_continuous(name=expression(paste("Latitude (",degree,")")),limits=c(4,12),expand=c(0,0))+
    geom_raster(data=ddf,aes(X,Y,fill=DEM),alpha=I(0.4))+
    scale_fill_gradientn(name="Altitude",colours = rainbow(100),breaks=b.dem)+
    coord_equal()
 print(p3)


#   Cleanup downloaded files and return to previous wd
unlink(tmp,recursive=T)
setwd(oldwd)

My questions are as follows:

Q1: How can I make the layers of p3 look like they do when plotted with base graphics in the example above?

Q2: How can I more sensibly specify colour scales so I don't have a ridiculous legend on the RHS?

Solution 1

Q1: You can't have different fill scales on different layers. One workaround is to use the fill aesthetic for the DEM and the alpha aesthetic for hillshade. Unfortunately, geom_raster doesn't seem to use the alpha aesthetic the way I expected. You can get the same effect with geom_tile, it just takes longer:

ggplot(hdf) +
  geom_raster(data=ddf,aes(X,Y,fill=DEM)) +
  scale_fill_gradientn(name="Altitude",colours = rainbow(100),breaks=b.dem) +
  geom_tile(aes(X,Y,alpha=Hill), fill = "grey20") +
  scale_alpha(range = c(0, 0.5)) +
  scale_x_continuous(name=expression(paste("Longitude (",degree,")")),
    limits=c(-4,2),expand=c(0,0)) +
  scale_y_continuous(name=expression(paste("Latitude (",degree,")")),
    limits=c(4,12),expand=c(0,0)) +
  coord_equal() 

Q2: Check out ?guide_colorbar. It doesn't work very nicely with your 100 colour breaks, but with fewer it's pretty good.

ggplot(hdf)+
  geom_raster(data=ddf,aes(X,Y,fill=DEM))+
  scale_fill_gradientn(name="Altitude",colours = rainbow(20))+
  guides(fill = guide_colorbar()) +
  geom_tile(aes(X,Y,alpha=Hill), fill = "grey20") +
  scale_alpha(range = c(0, 0.5)) +
  scale_x_continuous(name=expression(paste("Longitude (",degree,")")),
    limits=c(-4,2),expand=c(0,0)) +
  scale_y_continuous(name=expression(paste("Latitude (",degree,")")),
    limits=c(4,12),expand=c(0,0)) +
  coord_equal() 

Solution 2

The alpha in raster will be supported in the next version, so you can draw by:

ggplot(NULL, aes(X, Y)) + 
  geom_raster(data = ddf, aes(fill = DEM)) + 
  geom_raster(data = hdf, aes(alpha = Hill)) +
  scale_fill_gradientn(name="Altitude",colours = rainbow(20))+
  guides(fill = guide_colorbar()) +
  scale_alpha(range = c(0, 0.5), guide = "none") +
  scale_x_continuous(name=expression(paste("Longitude (",degree,")")), limits=c(-4,2),expand=c(0,0)) +
  scale_y_continuous(name=expression(paste("Latitude (",degree,")")), limits=c(4,12),expand=c(0,0)) +
  coord_equal()

anyway, very beautiful plot.

If you want to use this immediately, try installing from github:

library(devtools)
install_github("ggplot2", "kohske", "fix/geom-raster-alpha")

Note that geom_tile and geom_raster look different in some devices. Perhaps raster is better for your purpose.