Mapa cool

This is a document to show how to make a cool map from a Digital Elevation Model with both raster and rasterVis. First we load the libraries

require(rgdal)
require(raster)
require(rasterVis)

Now let’s go for a DEM

TO<-readGDAL("/Volumes/eltera/Users/adiez/Documents/MDT200-TOLEDO-h30.asc")

/Volumes/eltera/Users/adiez/Documents/MDT200-TOLEDO-h30.asc has GDAL driver AAIGrid 
and has 618 rows and 1124 columns

TO<-raster(TO)

Next step, build a color palette

#El azul es 8,129,242 (#0881F2) y habría que añadirlo a los valores negativos
rojo<-c(17,113,117,149,178,202,222,242,238,242,234,218,194,214,226,222,238,255,250,255,255,255,255)
verde<-c(48,153,170,190,214,226,238,238,222,206,182,157,141,157,174,186,198,206,218,222,230,242,255)
azul<-c(7,89,101,113,117,149,161,161,153,133,129,121,125,145,165,182,210,226,234,230,242,255,255)
feli<-rgb(rojo,verde,azul,max=255)

Now we made some derived products (slope, aspect, hill) from the Digital Elevation Model for further use, but first let’s to add the projection values to the raster

utm = "+init=epsg:25830"
geo = "+init=epsg:4326"
proj4string(TO)<-utm
sl<-terrain(TO)#slope
as<-terrain(TO,opt="aspect")
hi<-hillShade(sl,as)

Let’s check the four maps

par(mfrow=c(2,2))
plot(TO,col=feli,interpolate=TRUE,legend=FALSE)
plot(sl, col=grey(100:0/100), legend=FALSE, main='TOLEDO')

plot(as,col=grey(100:0/100), legend=FALSE)
plot(hi*-1,col=grey(100:0/100), legend=FALSE)

And now the composite map

plot(sl, col=grey(100:0/100), legend=FALSE, main='TOLEDO')
plot(hi, col=grey(100:0/100), legend=FALSE, alpha=0.3,add=TRUE)

plot(TO,col=feli, legend=TRUE, alpha=0.35,add=TRUE)

Now same thing with rasterVis, first we set the Theme

feliTh<-function (region = feli, ...) 
{
    theme <- rasterTheme(region = region, ...)
    theme
}

In ggplot we can save the graphics for further use

m3<-levelplot(TO,par.settings = feliTh,alpha.regions = 0.35,col.regions = feli)
m2<-levelplot(hi,par.settings=GrTheme,alpha.regions = 0.3,col.regions = gray(0:100/100),panel=panel.levelplot.raster)
m1<-levelplot(sl,par.settings=GrTheme,alpha.regions = 0.3,col.regions = gray(100:0/100))
m1<-levelplot(sl,par.settings = GrTheme,at = seq(-1,0, 0.01),interpolate=TRUE,panel=panel.levelplot.raster,col.regions = gray(100:0/100))

And the final map

par(mfrow=c(1,1))
m0<-m1+m2
m3+as.layer(m0, under = TRUE)