beautiful-thematic-maps: alternative code without using ggplot2

code_without_ggplot2.R

## https://timogrossenbacher.ch/2016/12/beautiful-thematic-maps-with-ggplot2-only/
library(raster)
library(sp)
library(rgdal)
library(viridisLite)
library(rasterVis)

## Read age data
data <- read.csv("input/avg_age_15.csv", stringsAsFactors = F)
## Read shapefile
gde_15 <- readOGR("input/geodata/gde-1-1-15.shp", layer = "gde-1-1-15")
## Join data with shapefile
mapData <- merge(gde_15, data, by.x = 'BFS_ID', by.y = 'bfs_id')
## Import the raster file 
relief <- raster("input/geodata/02-relief-georef-clipped-resampled.tif")

## Number of classes for the age data
no_classes <- 6

## Colors using the magma palette
pal <- rev(
    magma(n = no_classes,
          ## Limit the first and last colors
          begin = 0.1,
          end = 0.8))

## Compute breaks with quantiles
quantiles <- quantile(mapData$avg_age_15, 
                      probs = seq(0, 1,
                                  length.out = no_classes + 1))

## Not needed here, but a better alternative to original code (a "for"
## loop)
## labels <- paste(signif(quantiles[-7], 4),
##                 signif(quantiles[-1], 4),
##                 sep = '-')

## Create a trellis object using sp::spplot to display the
## SpatialPolygons object.
psp <- spplot(mapData["avg_age_15"],
              col.regions = pal,
              ## define the points where the colors change
              at = quantiles,
              ## set the border color and width
              col = 'lightgray', lwd = 0.2,
              par.settings = list(axis.line=list(col="transparent")),
              ## adjust the legend
              colorkey =
                  list(space = 'bottom',
                       height = 0.7,
                       labels = list(
                           at = quantiles,
                           labels = signif(quantiles, 3),
                           rot = 45)
                       )
              )

## Create a trellis object to display the relief using
## rasterVis::levelplot
plevel <- levelplot(relief,
                    maxpixels = 2e6,
                    margin = FALSE,
                    par.settings = GrTheme)

## Ok, ready to go: overlay these trellis objects to display the
## result
psp + as.layer(plevel, under = TRUE)

Very very admirable! Finally someone proving that using ggplot2 is not a must at all for high quality graphics in R. Thank you!