# load necessary packages:
library(maptools)
library(gstat)
library(rgdal)
# download the gridded maps:
download.file("http://spatial-analyst.net/book/system/files/meuse.zip", destfile=paste(getwd(), "meuse.zip", sep="/"))
grid.list <- c("ahn.asc", "dist.asc", "ffreq.asc", "soil.asc")
# unzip the maps in a loop:
for(j in grid.list){
fname <- zip.file.extract(file=j, zipname="meuse.zip")
file.copy(fname, paste("./", j, sep=""), overwrite=TRUE)
}

# load grids to R:
meuse.grid <- readGDAL(grid.list[1])
# fix the layer name:
names(meuse.grid)[1] <- sub(".asc", "", grid.list[1])
for(i in grid.list[-1]) {
   meuse.grid@data[sub(".asc", "", i[1])] <- readGDAL(paste(i))$band1
}
names(meuse.grid)
# reformat maps where needed:
meuse.grid$ffreq <- as.factor(meuse.grid$ffreq)
# pixels per class:
table(meuse.grid$ffreq)
meuse.grid$soil <- as.factor(meuse.grid$soil)
table(meuse.grid$soil)
str(meuse.grid)

# plot the maps together:
ffreq.plt <- spplot(meuse.grid["ffreq"], col.regions=grey(runif(length(levels(meuse.grid$ffreq)))), main="Flooding frequency classes")
dist.plt <- spplot(meuse.grid["dist"], col.regions=grey(rev(seq(0,1,0.025))), main="Distance to river")
soil.plt <- spplot(meuse.grid["soil"], col.regions=grey(runif(length(levels(meuse.grid$ffreq)))), main="Soil type classes")
ahn.plt <- spplot(meuse.grid["ahn"], col.regions=grey(rev(seq(0,1,0.025))), main="Elevation (cm)")
print(ffreq.plt, split=c(1,1,4,1), more=T)
print(dist.plt, split=c(2,1,4,1), more=T)
print(ahn.plt, split=c(3,1,4,1), more=T)
print(soil.plt, split=c(4,1,4,1), more=F)