hillShade.R

library(raster)
library(rasterVis)
library(maptools) ## for readShapeLines
library(colorspace) ## for terrain_hcl

## Get data from GADM and DIVA-GIS
setwd(tempdir())
download.file('http://www.gadm.org/data/shp/ETH_adm.zip', 'ETH_adm.zip')
unzip('ETH_adm.zip', exdir='.')
download.file('http://www.diva-gis.org/data/msk_alt/ETH_msk_alt.zip', 'ETH_msk_alt.zip')
unzip('ETH_msk_alt.zip', exdir='.')

## Read DEM and boundaries
proj <- CRS(' +proj=longlat +ellps=WGS84')
ethAdm <- readShapeLines('ETH_adm1.shp', proj4string=proj)
ethAlt <- raster('ETH_msk_alt')

## Compute shaded relief
slope <- terrain(ethAlt, 'slope')
aspect <- terrain(ethAlt, 'aspect')
ethHS <- hillShade(slope=slope, aspect=aspect,
                       angle=20, direction=30)


## DEM color theme: terrain colors and blue for the sea
terrainTheme <- modifyList(rasterTheme(region=terrain_hcl(n=15)),
                           list(panel.background=list(col='skyblue3')))
                           
## hill shade theme: gray colors with semitransparency
hsTheme <- modifyList(GrTheme(),                   
                      list(regions=list(alpha=0.3)))

## DEM alone
levelplot(ethAlt, par.settings=terrainTheme, margin=FALSE, colorkey=FALSE)
## DEM with shaded relief
levelplot(ethAlt, par.settings=terrainTheme, margin=FALSE, colorkey=FALSE) +
  levelplot(ethHS, par.settings=hsTheme)
## DEM with shaded relief and boundaries
levelplot(ethAlt, par.settings=terrainTheme, margin=FALSE, colorkey=FALSE) +
  levelplot(ethHS, par.settings=hsTheme) +
    layer(sp.lines(ethAdm, col='black', lwd=0.3))