R Script to plot hillShade data over and DEM raster using ggplot2 Uses gridExtra package to push ggplot2 objects to particular viewports

Hillshade DEM.R

# Script for question posted on Stack Overflow
#	Load relevant libraries
library(ggplot2)
library(raster)
library(gridExtra)


vplayout <- function(x, y) {
	viewport(layout.pos.row = x, layout.pos.col = y)
	}


#	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 <- paste(tmp,"\\GHA_HS.asc",sep="")
f2 <- paste(tmp,"\\GHA_DEM.asc",sep="")
download.file(url1,f1)	#File is ~ 5,655Kb
download.file(url2,f2)	#File is ~ 2,645Kb


#	Create rasters from downloaded files
hs <- raster("GHA_HS.asc")
dem <- raster("GHA_DEM.asc")


#	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)

#	This won't work because ggplot2 won't work with multiple continuous colour scales (deliberate design convention to avoid users making complicated plots which cloud data meaning)
#	However there are a few use cases where overplotting of data may be useful
#	We can do this using the gridExtra package...

grid.newpage()
pushViewport(viewport(layout = grid.layout(1, 1)))
print(p1 , vp = vplayout( 1 , 1 ) )
pushViewport(viewport(layout = grid.layout(1, 1)))
print(p2 , vp = vplayout( 1 , 1 ) )


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


#	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?

I ran this code but it produces this image:

Adding the following code in p2 to remove background, grids and panel borders doesn't help either :

Theme(
plot.background = element_blank()
,panel.grid.major = element_blank()
,panel.grid.minor = element_blank()
,panel.border = element_blank()
)

Hey @simonohanlon101 and @GlobalConservationist
I also struggled with this and started doing your approach. After all your steps, I tried a different thing, that is stacking both the DEM and hillshade, and then use the hillshade values as alpha for the ggplot.

rast <- stack(DEMs,hill)
rastdf <- rasterToPoints(rast)
 rastdf <- data.frame(rastdf)
colnames(rastdf) <- c("X","Y","DEM", "hillshade")

To costumize the alpha range, I reescaled the hillshade values. As I have some NA after obtaining the hillshade from the DEM. I made a custom function, that assigns a value to the NAs and I can also customize the range of transparency for a better effect:

rescale_hill <- function(x, minx, maxx){  # input is the vector to reescale, the minimum and maximum desired
  d <- numeric(length(x)) #make a new vector
   max_old<- max(x, na.rm = T) #find the min and max values
   min_old <- min(x, na.rm = T)
  for(i in 1:length(x)){
    if(is.na(x[i])== T){ d[i]<- minx} #for the NAs assign the min value
    else{
      d[i] <- ((maxx-minx)/(max_old- min_old ))*(x[i]-max_old) +maxx
      }
  }
  d
}

Then I can reescale the hillshade and get rid of NAs:

rastdf$hill <- rescale_hill(rastdf$hillshade, 1, 0.5)

And add it in ggplot:

ggplot()+
  geom_raster(data=rastdf, aes(X, Y, fill = DEM), alpha=I(rastdf$hill))
...

That was before:

And after:

Not great but at least gives a bit of a hillshade effect