Chapitre 3 Cartographie thématique avancée

3.1 Les cartes de discontinuités

Ce type de représentation permet de souligner cartographiquement les discontinuités territoriales d’un phénomène. L’accent est porté sur ce qui distingue des territoires. Pour chaque frontière nous calculons le rapports ou la différence des valeurs des polygones de part et d’autre. Puis nous représentons la frontière par un figuré d’autant plus épais que la différence est forte. Il est souvent bénéfique de coupler ce type de représentation à une représentation choroplèthe (pour comprendre le sens des discontinuités).

Ces cartes ne sont pas évidentes à paramétrer. Le choix des critères (seuil, type de différences…) influence fortement la représentation. En fonction du maillage utilisé la lisibilité de la carte peut être faible.
Ces représentations sont très puissantes pour montrer les inégalités.

La fonctions getBorder() du package cartography permet de construire une couche des frontières terrestres. La fonction discLayer() permet d’afficher les discontinuités.

library(sf)
library(cartography)
mtq <- st_read(system.file("gpkg/mtq.gpkg", package="cartography"), 
               quiet = TRUE)
# Get borders
mtq_bord <- getBorders(x = mtq)
# Plot polygons
plot(st_geometry(mtq), border = NA, col = "grey60")
# Plot borders
plot(
  st_geometry(mtq_bord), 
  col = sample(x = rainbow(nrow(mtq_bord))), 
  lwd = 3, 
  add = TRUE
)
layoutLayer("Frontières inter-communales",tabtitle=TRUE, north = TRUE,
            author= "INSEE 2016", sources="", frame=FALSE, scale = 5)

# Compute the population density (inhab./km2) using sf::st_area()
mtq$POPDENS <- as.numeric(1e6 * mtq$POP / st_area(mtq))
# plot municipalities (only the backgroung color is plotted)
plot(st_geometry(mtq), col = NA, border = NA, bg = "lightblue1", 
     xlim = c(690574, 745940))
# Plot the population density with custom breaks
choroLayer(x = mtq, var = "MED",
           breaks = c(min(mtq$MED), seq(13000, 21000, 2000), max(mtq$MED)),
           col = carto.pal("green.pal", 6),border = "white", lwd = 0.5, 
           legend.pos = "topright", legend.title.txt = "Median Income\n(euros)",
           add = TRUE)
# Plot discontinuities
discLayer(
  x = mtq_bord, 
  df = mtq, 
  var = "MED",
  type = "rel", 
  method = "geom", 
  nclass = 3,
  threshold = 0.4,
  sizemin = 0.7, 
  sizemax = 6, 
  col = "red4",
  legend.values.rnd = 1, 
  legend.title.txt = "Relative\nDiscontinuities", 
  legend.pos = "right",
  add = TRUE
)
# Layout
layoutLayer(title = "Wealth Disparities in Martinique, 2015", 
            author =  paste0("cartography ", packageVersion("cartography")),
            sources = "Sources: Insee and IGN, 2018",
            frame = FALSE, scale = 5, tabtitle = TRUE,theme = "grey.pal")
# north arrow
north(pos = "topleft")

3.2 Les grilles régulières

La méthode du carroyage consiste à découper l’espace géographique en un maillage formé de carrés réguliers dans une projection donnée. La donnée est répartie sur ce quadrillage régulier au prorata de la surface représentée. Le quadrillage permet ainsi de s’affranchir des mailles administratives.

Ces représentation induisent une perte de précision. Les maillages produit n’ont pas de signification. La version simple (les valeurs sont redistribuées au prorata de la surface), implique une equirépartition du phénomène dans chaque unités.
La comparaison de maillages différents, à plusieurs dates ou de différentes sources est rendue possible.

La fonction getGridLayer() du package cartography permet de construire ces grilles régulières.

library(sf)
library(cartography)
mtq <- st_read(system.file("gpkg/mtq.gpkg", package="cartography"), 
               quiet = TRUE)
# Plot dentsity of population 
mtq$dens <- mtq$POP / as.numeric(st_area(mtq) / (1000 * 1000)) 
bks <- getBreaks(v = mtq$dens, method = "q6")
cols <- carto.pal(pal1 = "taupe.pal", n1 = 6)
choroLayer(
  x = mtq, 
  var = "dens", 
  breaks = bks, 
  border = "burlywood3", 
  col = cols, 
  legend.pos = "topright", 
  legend.values.rnd = 1,
  legend.title.txt = "Densité de population\n(hab/km2)"
)
layoutLayer("Population en Martinique",tabtitle=TRUE, 
            author= "INSEE 2016", sources="", frame=FALSE, scale = 5)
north(pos = "topleft")

# Création de la grille
mygrid <- getGridLayer(
  x = mtq, 
  cellsize = 10000 * 10000, 
  type = "hexagonal", 
  var = "POP"
)
## conversion from square meter to square kilometers
mygrid$densitykm <- mygrid$POP / (mygrid$gridarea / (1000 * 1000)) 
choroLayer(
  x = mygrid, 
  var = "densitykm", 
  breaks = bks,
  border = "burlywood3", 
  col = cols, 
  legend.pos = "topright", 
  legend.values.rnd = 1,
  legend.title.txt = "Densité de population\n(hab/km2)"
)
plot(st_geometry(mtq), lwd = 0.2, add=TRUE, border = "#ffffff75")
layoutLayer(
  title = "Population en Martinique",
  tabtitle=TRUE, 
  author= "INSEE 2016", 
  sources="", 
  frame=FALSE, 
  scale = 5
)
north(pos = "topleft")

3.3 Le lissage spatial

L’idée principale du lissage est de filtrer l’information pour révéler des structures spatiales sous-jacentes. C’est un ensemble de méthodes qui consistent à affecter aux points que l’on observe une valeur prenant en compte les valeurs de leur voisinage. Il existe plusieurs méthodes de lissage (kde, potentiels…) plus ou moins paramétrables. Cette méthode permet de passer d’une représentations de données ponctuelles vers la représentation d’une surface continue.

Il est difficile de paramétrer correctement les fonctions de lissages.
Elles doivent s’appuyer sur des hypothèses de comportement dans l’espace.
La compréhension par un public large n’est pas évidente, il faut alors simplifier les légendes, la présentation de la méthode.
Permet de faire ressortir des phénomènes spatiaux sous-jacents invisibles directement.
Les cartes produites attirent l’oeil par leur originalité.
Cette méthode permet de passer d’une représentation ponctuelle ou discontinue (dans un maillage) à une représentation continue s’affranchissant des maillages existants.

La méthode utilisée ici est celle de l’estimation par noyau (KDE).

library(sf)
library(spatstat)
library(maptools)
library(raster)
# Import des données
mtq <- st_read(system.file("gpkg/mtq.gpkg", package="cartography"), 
               quiet = TRUE)
# resto <- st_read(dsn = "data/resto.gpkg", quiet = TRUE)
sigma = 4000
res = 250
# Define an observation window
w <- as.owin(as(st_union(st_buffer(mtq, 100)), "Spatial"))
# sf to coords
pts <- st_coordinates(st_centroid(st_geometry(mtq)))
# Coords to ppp
p <- ppp(pts[,1], pts[,2], window = w, marks = mtq$POP)
# Compute KDE
dens <- density.ppp(p, weights = p$marks, sigma = sigma, eps = res)
# Image to raster (+ proj & km2)
result <- raster(dens, crs = st_crs(mtq)[[2]]) * 1000000
result <- mask(result, mtq)
# compute breaks
bks <- unique(getBreaks(values(result), nclass = 10, method = "q6"))
# Color ramp
cols <- hcl.colors(10, "viridis")[2:8]
# Plot the map
plot(st_geometry(mtq), col = NA, border = NA, bg = "lightblue1")
plot(result, breaks = bks, col=cols, add = T,legend=F)
plot(st_union(st_geometry(mtq)), add = T, lwd = 1.5, border = "lightblue1")
legendChoro(
  pos = "topright",
  title.txt = "Densité de\npopulation\n(hab./km2)",
  breaks = bks, 
  nodata = FALSE,
  values.rnd = -1,
  col = cols
)
layoutLayer(title = "Densité de population", tabtitle = TRUE,
            author = "INSEE, 2016", sources = "", 
            frame = FALSE, scale = 5)
north(pos = "topleft")

3.4 Cartes en 3D

3.4.1 linemap

Le package linemap (Giraud 2017a) permet de réaliser des cartes composées de lignes.

library(linemap)
library(sf)
data("popOcc")
data("occitanie")
opar <- par(mar=c(0,0,0,0), bg = "ivory2")
bb <- st_bbox(occitanie)
plot(st_geometry(occitanie), col="ivory1", border = NA)
linemap(
  x = popOcc, 
  var = "pop", 
  k = 2.5, 
  threshold = 50,
  col = "ivory1", 
  border = "ivory4", 
  lwd = 0.6, 
  add = TRUE
)
text(x = bb[1], y = bb[4],adj = c(0,1),
     labels = "Répartition de la\npopulation\nen Occitanie",  
     col = "ivory4", font = 2,  cex = 1.8)
# add sources
mapsources <-"Timothée Giraud\nlinemap 0.1.0\nDonnées carroyées à 1 kilomètre, INSEE 2010"
text(x = bb[3], y = bb[2],labels = mapsources,  
     col = "ivory4", font = 3, adj = c(1,0), cex = 0.6 )

3.4.2 Relief Tanaka

Nous utilisons le package tanaka (Giraud 2020) qui met à disposition une méthode (Tanaka 1950) utilisée pour améliorer la perception du relief.

library(raster)
library(cartography)
library(sf)
library(tanaka)
mtq <- st_read(system.file("gpkg/mtq.gpkg", package="cartography"), 
               quiet = TRUE)
# use WGS84 proj
mtq_latlon <- st_transform(mtq, 4326)
# import raster
ras <- raster("data/srtm_24_10.tif")
# crop on martinique area
mtq_ras <- crop(ras, st_bbox(mtq_latlon)[c(1,3,2,4)])
# aggregate the raster
mtq_ras <- aggregate(mtq_ras, fact=4, fun=mean)
mtq_ras <- projectRaster(mtq_ras, crs=st_crs(mtq)$proj4string)
# break values
bks <- c(seq(0,1300,100),1302)
# custom palette
cols <- c("#5D9D52", "#8DBC80", "#B8D9A9", "#FDEBBE", "#F7E0AC", "#F2D69B", 
          "#EDCC8A", "#E8C279", "#E2B563", "#DBA84C", "#D49B36", "#BA8428", 
          "#9A6A1E", "#7B5114")
# map
plot(st_geometry(mtq), col = NA, border = NA, bg = "lightblue1")
tanaka(
  x = mtq_ras,
  breaks = bks,
  col = cols, 
  legend.pos = "NA", 
  mask = mtq,
  add = TRUE
)
legendChoro(pos = c(689000,1598000 ), breaks = bks, col = cols, nodata = F,
            title.txt = "Elevation\n(metres)", cex = 1)
layoutLayer(title = "Martinique Relief", north = T,
            sources = 'T. Giraud, 2019', author = "SRTM, 2018", 
            col = "lightblue", 
            tabtitle = T, coltitle = "black")

The tanaka package

3.4.3 Rayshader

Le package rayshader (Morgan-Wall 2018) permet de réaliser de belles cartes en relief. L’export des images n’est pas évident, il s’agit ici d’une simple capture d’écran.

library(sf)
library(raster)
library(rayshader)
mtq <- st_read(system.file("gpkg/mtq.gpkg", package="cartography"), 
               quiet = TRUE)
st_geometry(mtq) <- st_buffer(st_geometry(mtq), 5000)
mtq_latlon <- st_transform(mtq, 4326)
ras <- raster("data/srtm_24_10.tif")
mtq_ras <- crop(ras, st_bbox(mtq_latlon)[c(1,3,2,4)])
mtq_ras <- projectRaster(mtq_ras, crs=st_crs(mtq)$proj4string)
elmat = matrix(extract(mtq_ras,extent(mtq_ras),buffer=1000),
               nrow=ncol(mtq_ras),ncol=nrow(mtq_ras))
elmat[is.na(elmat)] <- 0
raymat = ray_shade(elmat,lambert = TRUE,anglebreaks = 85,sunangle = 125)
ambmat = ambient_shade(elmat,anglebreaks =  85)
elmat %>%
  sphere_shade(texture = "imhof1",sunangle = 125) %>%
  add_water(detect_water(elmat), color="desert") %>%
  add_shadow(raymat,0.5) %>%
  add_shadow(ambmat,0.5) %>%
  plot_3d(elmat,zscale=25,fov=10,theta=-15,phi=70, background="black",
          zoom=.5, windowsize =c(900, 900))

3.5 Les cartogrammes

L’anamorphose classique est une représentation des États (ou de mailles quelconques) par des rectangles ou des polygones quelconques en fonction d’une quantité qui leur est rattaché. (…) On s’efforce de garder l’arrangement général des mailles ou la silhouette du continent."
Brunet, Ferras, and Théry (1993)

3 types d’anamorphoses ou cartogrammes sont ici présentés :

Un cours complet sur les anamorphoses : Les anamorphoses cartographiques (Lambert 2015).
Faire des cartogrammes avec R

Pour réaliser les cartogrammes nous utilisons le package cartogram (Jeworutzki 2018).

3.5.1 Les cartogrammes de Dorling

Les territoires sont représentés par des figurés (cercles, des carrés ou des rectangles) ne se recouvrant pas dont les surfaces sont proportionnelles à une variable. Les positions des figurés sont définie selon les positions de départ.

(McCormick 2017)

On identifie assez mal l’espace.
On peut nommer les cercles pour se repérer et/ou s’aider de la couleur pour faire apparaitre des clusters et mieux identifier les blocks géographiques.
La perception des quantités est très bonne.
Les tailles de cercles sont vraiment comparables. 
library(cartography)
library(cartogram)
library(sf)
mtq <- st_read(system.file("gpkg/mtq.gpkg", package="cartography"), 
               quiet = TRUE)
mtq_dorling <- cartogram_dorling(x = mtq, weight = "POP", k = 12) 
plot(st_geometry(mtq_dorling), col = "#940000", border= "white", bg = "lightblue")
labelLayer(x = mtq_dorling, txt = "LIBGEO", overlap = FALSE, show.lines = FALSE, 
           halo = TRUE, r = .15)
layoutLayer("Population en Martinique - Cartogramme de Dorling",tabtitle=TRUE, 
            author= "INSEE 2016", sources="", frame=FALSE, scale = NULL)

Le paramètre k permet de faire varier le facteur d’expansion des cercles.

3.5.2 Les cartogrammes non continus

La taille des polygones est proportionnelle à une variable. L’agencement des polygones les uns par rapport aux autres est conservée. La forme des polygones est ressemblante.

(Cauvin, Escobar, and Serradj 2013)

La topologie des régions est perdue.
La conservation de la forme des polygones est optimale. 
library(cartography)
library(cartogram)
library(sf)
mtq <- st_read(system.file("gpkg/mtq.gpkg", package="cartography"), 
               quiet = TRUE)
mtq_ncont <- cartogram_ncont(x = mtq, weight = "POP", k = 1.5) 
plot(st_geometry(mtq), col = NA, border = "white", lwd = 0.5, bg = "lightblue")
plot(st_geometry(mtq_ncont), col = "#940000", border= "white", add=TRUE)
layoutLayer("Population en Martinique - Cartogramme de Olson",tabtitle=TRUE, 
            author= "INSEE 2016", sources="", frame=FALSE, scale = NULL)

Le paramètre k permet de faire varier le facteur d’expansion des polygones.

3.5.3 Les cartogrammes continus

La taille des polygones est proportionnelle à une variable. L’agencement des polygones les uns par rapport aux autres est conservée. Pour conserver la contiguité, la forme des polygones est fortement transformée.

(Paull and Hennig 2016)

La forme des polygones est fortement distordue.
C’est une “vraie carte de géographie” : la topologie et la contiguité sont conservées. 
library(cartography)
library(cartogram)
library(sf)
mtq <- st_read(system.file("gpkg/mtq.gpkg", package="cartography"), 
               quiet = TRUE)
mtq_cont <- cartogram_cont(x = mtq, weight = "POP", prepare = "none") 
plot(st_geometry(mtq_cont), col = "#940000", border= "white", bg = "lightblue")
layoutLayer("Population en Martinique - Cartogramme de Dougenik",tabtitle=TRUE, 
            author= "INSEE 2016", sources="", frame=FALSE, scale = NULL)

3.5.4 Forces et faiblesses des cartogrammes

Les cartogrammes sont des représentations cartographiques perçues comme innovante (bien que la methode date de 40 ans). Ces images très généralisées rendent bien compte des quantités et des gradiants. Ce sont de vraies images de communication qui provoquent, suscitent l’intérêt, véhiculent un message fort, interpellent.

Mais les cartogrammes induisent une perte des repères visuels (difficile de retrouver son pays, ou sa région sur la carte), demandent un effort de lecture qui peut être important et ne permettent pas de gérer des données manquantes.

3.6 Cartographie interactive

Le package mapview (Appelhans et al. 2018) utilise le package leaflet (Cheng, Karambelkar, and Xie 2018) basé sur la librairie javascript Leaflet (Agafonkin 2015) pour créer des cartes interactives. La création de carte est assez aisée, la documentation est un peu dense.

Utilisation simple

library(sf)
library(leaflet.providers)
library(mapview)
mtq <- st_read(system.file("gpkg/mtq.gpkg", package="cartography"), 
               quiet = TRUE)
resto <- st_read(dsn = "data/resto.gpkg", quiet = TRUE)
mapview(resto) + mapview(mtq)


Utilisation personnalisée

mapview(
  resto, 
  map.types = "OpenStreetMap", 
  col.regions = "#940000", 
  label = resto$name,
  alpha.regions=90,
  color = "white", 
  legend = TRUE, 
  layer.name = "Restaurants", 
  homebutton = FALSE, 
  lwd = 0.5, 
  popup = NA
) +
  mapview(
    mtq, 
    col.regions = "lightblue", 
    color = "white", 
    legend = TRUE, 
    label = mtq$LIBGEO, 
    alpha.regions= 0.5, 
    map.types = "OpenStreetMap",
    lwd = 0.5, 
    layer.name = "Communes",
    homebutton = FALSE, 
    popup = leafpop::popupTable(mtq, 
                       zcol = 1:4, 
                       row.numbers = FALSE, 
                       feature.id = FALSE)
  )