Dans un très complet billet de blog, Timo Grossenbacher décrit la création d’une carte thématique avec ggplot2.
https://timogrossenbacher.ch/2016/12/beautiful-thematic-maps-with-ggplot2-only/
Il y détaille les étapes nécessaires à la création de cette carte en évoquant par exemple des questions de reproductibilité, l’import et la sélection des données, le choix d’une discrétisation, le choix d’une gamme de couleur, l’utilisation de données raster (pour le relief), les questions d’habillage et de design.

Dans les commentaires quelqu’un propose une version de cette même carte construite avec lattice (beautiful-thematic-maps: alternative code without using ggplot2) et une version interactive a aussi été développée (Thematic Interactive Map).

Voici donc ma contribution sans ggplot2, ni lattice, ni interactivité :-).

# packages nécessaires
library(rgdal)
library(raster)
library(cartography)

# import des données 
gde_15 <- readOGR("input/geodata/gde-1-1-15.shp", layer = "gde-1-1-15")
relief <- raster("input/geodata/02-relief-georef-clipped-resampled.tif")
age <- read.csv("input/avg_age_15.csv", stringsAsFactors = F)

# pour l'export de la carte
# png("chuiche.png", width = 600, height = 440, res = 96)

# j'utilise la police "liberation" qui est une police libre 
# https://en.wikipedia.org/wiki/Liberation_fonts
par(mar = c(3,0.2,1.2,0.2), bg = "#f5f5f2", family = "liberation")

# affichage du relief
bb <- extent(relief)
plot.new()
plot.window(xlim = bb[1:2], ylim = bb[3:4], asp = TRUE)
plot(relief, legend = FALSE, axes = FALSE, box = FALSE,
     col = colorRampPalette(c("grey50","white"))(25),
     add = TRUE)

# discrétisation de la variable
# j'utilise la méthode "q6" qui permet de bien 
# distinguer les valeurs extrèmes
# ?getBreaks
bks <- getBreaks(v = age$avg_age_15, method = "q6")

# choix d'une gamme de couleur
cols <- c("#FCC696", "#F88F79", "#DA6688", 
          "#AB5596", "#7D4495", "#4D3E71")

# affichage de l'age moyen
choroLayer(spdf = gde_15, df = age, 
           spdfid = "BFS_ID", dfid = "bfs_id", 
           var = "avg_age_15", method = "q6", col = cols, 
           border = "white", lwd = 0.2, legend.pos = "n",
           add = TRUE )

# ajout d'une échelle et d'une flèche nord
north(col = "#4D3E71")
barscale(size = 20)

# ajout du titre, sous titre et sources.
mtext(text = "Switzerland's regional demographics", 
      side = 3, line = 0,adj = 0, col = "#4D3E71", font = 2)
mtext(text = "Average age in Swiss municipalities, 2015", 
      side = 3, line = - 1,adj = 0, cex = 0.8, col = "#4D3E71" )
mtext(text = "Map CC-BY-SA; Author: Timothée Giraud, Geometries: ThemaKart, BFS; Data: BFS, 2016; Relief: swisstopo, 2016",
      side = 1, line = 2,adj = 1, cex = 0.5, col = "grey40" )


# J'utilise la méthode décrite ici (https://rgeomatic.hypotheses.org/878)
# pour insérer un histogramme dans la carte. 
par(fig = c(grconvertX(c(480000, 840000), from="user", to="ndc"),
            grconvertY(c(55000,  77000) , from="user", to="ndc")),
    mar = c(0,0,0,0),  
    new = TRUE)

# ajout d'un histogramme qui servira de légende
hist(age$avg_age_15, probability = T, breaks = bks, 
     col = cols, axes = F, xlab="", ylab = "",
     main = "", border = cols)
axis(side = 1, at = bks[c(1,3,5,7)], labels = round(bks[c(1,3,5,7)], 1),
     cex.axis = 0.5, tick = FALSE, line = -1.1, col.axis = "grey10")
axis(side = 1, at = bks[c(2,4,6)], labels = round(bks[c(2,4,6)], 1),
     cex.axis = 0.5, tick = FALSE, line = -1.1, col.axis  = "grey20")
text(x = 33.1, y = 0.16, label = "Average age", cex = 0.8, adj = c(0,1), 
     font = 3, col = "#4D3E71")

Une bonne partie de la discussion autour de la carte initiale porte sur le choix d’une bonne discrétisation et sur la création d’une légende rendant compte de cette discrétisation. La légende que je propose rend plutôt bien compte de la distribution de la variable et de sa discrétisation, certes au prix d’un poids visuel assez important…

Vous pouvez retrouver ce script et les données initiales dans un fork du dépôt de Timo Grossenbacher.