Lorsque des cartes de la France sont diffusées, je suis toujours surpris par l’absence des départemtens d’outre-mer (ex1, ex2 ex3).
Le but ici est de créer une visualisation de tous les départements français, sans que ceux-ci ne soient à leur véritable position géographique.
Tous les fichiers au format sf
peuvent être récupérés sur le site Global Administrative Areas.
Les données sont géoréférencées au format WGS 84.
# Récupération des données départementales
fra <- readRDS("gadm36_FRA_2_sf.rds")
mtq <- readRDS("gadm36_MTQ_0_sf.rds")
glp <- readRDS("gadm36_GLP_0_sf.rds")
reu <- readRDS("gadm36_REU_0_sf.rds")
myt <- readRDS("gadm36_MYT_0_sf.rds")
guy <- readRDS("gadm36_GUF_0_sf.rds")
# Modifier le nom des départements d'outre-mer (en français)
guy$NAME_0 <- "Guyane"
reu$NAME_0 <- "La Réunion"
head(fra)
st_crs(fra)
## Coordinate Reference System:
## EPSG: 4326
## proj4string: "+proj=longlat +datum=WGS84 +no_defs"
La projection utilisée sur ces cartes est de type “long/lat”, qui est un système non projeté.
Pour une visualisation de cartes, il est généralement préférable d’utiliser un système projeté.
On peut modifier les systèmes de projection de coordonnées, pour obtenir des affichages de cartes meilleurs.
Pour cela, on va utiliser le système utilisé par OpenStreetMap entre autres.
# Mofifier le système de coordonnées pour avoir un système "WGS 84 / Pseudo-Mercator"
# (facilite la visualisation des cartes)
list(fra, mtq, glp, reu, myt, guy) %>%
purrr::map(.f = st_transform, crs = 3857) %>%
setNames(nm = c("fra", "mtq", "glp", "reu", "myt", "guy")) %>%
list2env(envir = .GlobalEnv)
## <environment: R_GlobalEnv>
Quelques mises en carte.
# France métropolitaine
plot(st_geometry(fra))
# Martinique
plot(st_geometry(mtq))
Si on essaie de joindre ensemble toutes ces données et de les afficher, elles s’afficheront avec leurs “vraies” coordonnées, et donc espacées sur le monde entier.
# Exemple avec la métropole et la Martinique
fra %>%
select(colnames(mtq)) %>%
rbind(mtq) %>%
st_geometry() %>%
plot(main = "Métropole et Martinique à leur position géographique")
Le but est de visualiser les départements d’outre-mer avec un affichage “proche” de la métropole.
ggplot
)Dans un premier temps, on crée une carte de la métropole (colorée par régions).
fra_map <- ggplot(fra) +
geom_sf(aes(fill = NAME_1), show.legend = FALSE, color = "white", size = 0.2) +
# enlever l'affichage des coordonnés et de la grille
coord_sf(datum = NA, expand = FALSE) +
scale_fill_viridis(discrete = TRUE, option = "D") +
ggtitle("Carte des départements français") +
theme(panel.background = element_blank(),
plot.title = element_text(hjust = 0.5, size = 16))
fra_map
De la même façon, on va créer une fonction permettant de créer une carte pour chaque DOM.
plot_map <- function(sf) {
ggplot(sf) +
geom_sf(fill = viridis(n = 1, option = "D"), alpha = 0.7, col = "white", size = 0.4) +
# enlever l'affichage des coordonnés
coord_sf(datum = NA, expand = FALSE) +
ggtitle(sf$NAME_0) +
theme(panel.background = element_blank(),
plot.title = element_text(hjust = 0.5, size = 10))
}
Essai avec la Martinique.
plot_map(mtq)
On va maintenant rapprocher les cartes des DOM près de la métropole.
Je me suis inspiré de ce site, qui a procédé de même pour afficher ensemble les Etats-Unis, l’Alaska et Hawaii.
Afin de faciliter le positionnement des cartes des DOM, on peut extraire les coordonnées du contour de la carte de la métropole.
(fra_bbox <- st_bbox(fra))
## xmin ymin xmax ymax
## -572599.8 5061695.3 1064260.7 6637122.4
On se sert de ces valeurs de fra_bbox
pour ajuster les cartes au mieux.
Les coordonnées xmin
, xmax
, ymin
et ymax
ont été trouvées en tâtonnant pour obtenir un bon ratio d’affichage.
fra_map +
annotation_custom(
grob = ggplotGrob(plot_map(mtq)),
xmin = fra_bbox[1],
xmax = fra_bbox[1] - 240000,
ymin = fra_bbox[2],
ymax = fra_bbox[2] + 400000
) +
annotation_custom(
grob = ggplotGrob(plot_map(glp)),
xmin = fra_bbox[1],
xmax = fra_bbox[1] - 240000,
ymin = fra_bbox[2],
ymax = fra_bbox[2] + 1000000
) +
annotation_custom(
grob = ggplotGrob(plot_map(reu)),
xmin = fra_bbox[1],
xmax = fra_bbox[1] - 240000,
ymin = fra_bbox[2],
ymax = fra_bbox[2] + 1550000
) +
annotation_custom(
grob = ggplotGrob(plot_map(myt)),
xmin = fra_bbox[1],
xmax = fra_bbox[1] - 240000,
ymin = fra_bbox[2],
ymax = fra_bbox[2] + 2200000
) +
annotation_custom(
grob = ggplotGrob(plot_map(guy)),
xmin = fra_bbox[1],
xmax = fra_bbox[1] - 240000,
ymin = fra_bbox[2],
ymax = fra_bbox[2] + 2900000
)
Avec cette méthode, on obtient la carte souhaitée.
On peut aussi remarquer que les échelles des DOM ne sont pas respectées et il est difficile de mettre en place les bonnes échelles avec les fonctions annotation_custom()
.
De plus, si on veut modifier les différentes cartes, il faut le faire en amont de la carte finale, puisque cet objet est de type ggplot
(donc un graphe), et non plus sf
(objet spatial).
tmap
)De la même façon, on peut construire l’assemblage de cartes en utilisant la librarie tmap
(comme expliqué ici).
library(tmap)
fra_map <- tm_shape(fra, projection = 4326) +
tm_borders(col = "white", lwd = 0.5) +
tm_polygons("NAME_1", palette = "viridis", option = "D") +
tm_layout(frame = FALSE, legend.show = FALSE,
title = "Carte des départements français")
mtq_map <- tm_shape(mtq) +
tm_polygons(col = viridis(n = 1, option = "D"), alpha = 0.7) +
tm_layout(frame = FALSE)
glp_map <- tm_shape(glp) +
tm_polygons(col = viridis(n = 1, option = "D"), alpha = 0.7) +
tm_layout(frame = FALSE)
reu_map <- tm_shape(reu) +
tm_polygons(col = viridis(n = 1, option = "D"), alpha = 0.7) +
tm_layout(frame = FALSE)
myt_map <- tm_shape(myt) +
tm_polygons(col = viridis(n = 1, option = "D"), alpha = 0.7) +
tm_layout(frame = FALSE)
guy_map <- tm_shape(guy) +
tm_polygons(col = viridis(n = 1, option = "D"), alpha = 0.7) +
tm_layout(frame = FALSE)
fra_map
## Warning: One tm layer group has duplicated layer types, which are omitted.
## To draw multiple layers of the same type, use multiple layer groups (i.e.
## specify tm_shape prior to each of them).
print(mtq_map, vp = grid::viewport(x = 0.12, y = 0.13, width = 0.15, height = 0.15))
print(glp_map, vp = grid::viewport(x = 0.12, y = 0.3, width = 0.15, height = 0.15))
print(reu_map, vp = grid::viewport(x = 0.12, y = 0.45, width = 0.15, height = 0.13))
print(myt_map, vp = grid::viewport(x = 0.12, y = 0.6, width = 0.15, height = 0.15))
print(guy_map, vp = grid::viewport(x = 0.12, y = 0.76, width = 0.15, height = 0.15))
Je trouve que le résultat est moins efficace qu’avec ggplot
, notamment car les titres ne sont pas ajustés (on peut ajouter le nom de départements d’outre-mer, mais ceux-ci chevauchent les cartes).
Le problème d’échelle des départements n’est pas non plus réglé.
En revanche, le code est beaucoup plus succinct.
sf
Différentes transformations sont possibles sur les objets de classe sf
:
Des explications sont disponibles ici.
#######
# Déplacer une géométrie
#######
mtq2 <- mtq
# extraire la géométrie
mtq2_sfc <- st_geometry(mtq)
# ici, déplacement de 80km vers l'est et 20km vers le nord
st_geometry(mtq2) <- mtq2_sfc + c(80000, 20000)
# redéfinir le système de coordonnées
st_crs(mtq2) <- st_crs(mtq)
# pour affichage
mtq2$NAME_0 <- "Martinique_move"
ggplot() +
geom_sf(data = mtq) +
geom_sf(data = mtq2) +
geom_sf_text(data = mtq, aes(label = NAME_0)) +
geom_sf_text(data = mtq2, aes(label = NAME_0)) +
coord_sf(datum = st_crs(mtq))
#######
# homothétie / changement d'échelle d'une géométrie
#######
mtq2 <- mtq
# extraire la géométrie
mtq2_sfc <- st_geometry(mtq)
# extraire le centroîde de la géométrie
mtq2_centroid <- st_centroid(mtq2_sfc)
# ici, réduction d'échelle avec un facteur 0.3
st_geometry(mtq2) <- (mtq2_sfc - mtq2_centroid) * 0.3 + mtq2_centroid
# redéfinir le système de coordonnées
st_crs(mtq2) <- st_crs(mtq)
# pour affichage
mtq2$NAME_0 <- "Martinique_scale"
ggplot() +
geom_sf(data = mtq, fill = "lightsalmon") +
geom_sf(data = mtq2, fill = "lightblue") +
geom_sf_label(data = mtq, aes(label = NAME_0), nudge_x = -10000, nudge_y = 15000, fill = "lightsalmon") +
geom_sf_label(data = mtq2, aes(label = NAME_0), fill = "lightblue")
Pour combiner les 2 transformations :
#######
# translation + homothétie d'une géométrie
#######
mtq2 <- mtq
# extraire la géométrie
mtq2_sfc <- st_geometry(mtq)
# extraire le centroîde de la géométrie
mtq2_centroid <- st_centroid(mtq2_sfc)
# combinaision des 2 transformations
st_geometry(mtq2) <- (mtq2_sfc - mtq2_centroid) * 0.3 + mtq2_centroid + c(80000, 20000)
# redéfinir le système de coordonnées
st_crs(mtq2) <- st_crs(mtq)
# pour affichage
mtq2$NAME_0 <- "Martinique_move+scale"
ggplot() +
geom_sf(data = mtq, fill = "lightsalmon") +
geom_sf(data = mtq2, fill = "lightblue") +
geom_sf_label(data = mtq, aes(label = NAME_0), fill = "lightsalmon") +
geom_sf_label(data = mtq2, aes(label = NAME_0), nudge_x = -12000, nudge_y = -10000, fill = "lightblue")
On peut combiner le tout dans une fonction, afin de déplacer les objets (en modifiant éventuellement leur échelle pour avoir un meilleur affichage).
(le code suivant est inspiré de : https://github.com/Nowosad/us-map-alternative-layout)
# Fonction pour déplacer une géométrie
place_geometry <- function(geometry, position, scale = 1) {
# prend en entrée une géométrie existante : 'geometry'
# déplace cette géométrie au point 'position'
# par défaut, pas de changement d'échelle
# Nouvelle géométrie
output_geometry <- (geometry - st_centroid(geometry)) * scale + st_centroid(geometry) +
# translation
position
# Ajouter le système de coordonnées
st_crs(output_geometry) <- st_crs(geometry)
return(output_geometry)
}
Essai avec la Martinique
On peut utiliser les “bounding box” de l’île et de la métropole pour calculer la position finale sur la carte de la métropole.
L’idée est de faire correspondre les xmin
et ymin
de l’île avec ceux de la métropole.
On retire ensuite une certaine valeur sur l’axe X pour décaler les DOM sur la gauche, et on ajoute une certaine valeur que l’axe Y pour “empiler” les cartes des DOM.
(note : ces valeurs sont trouvées par tâtonnement pour trouver un affichage correct).
# position finale de la Martinique que la carte de la métropole
mtq_position2 <- c(st_bbox(fra)$xmin - st_bbox(mtq)$xmin # correspondance des 'xmin'
- 150000, # décalage axe X
st_bbox(fra)$ymin - st_bbox(mtq)$ymin # correspondance des 'ymin'
+ 100000) # décalage axe Y
# Modifier la géométrie de l'ojet 'sf' :
# on utilise la fonction précédente, en utilisant
# la géométrie originale de la Martinique, et en la
# déplaçant au nouveau point défini, avec un changement
# d'échelle.
mtq2 <- mtq %>%
mutate(geometry = place_geometry(geometry = st_geometry(mtq),
position = mtq_position2,
scale = 2.5))
ggplot() +
geom_sf(data = fra) +
geom_sf(data = mtq2) +
geom_sf_text(data = mtq2, aes(label = NAME_0), nudge_y = (st_bbox(mtq2)[4] - st_bbox(mtq2)[2]) / 2 + 50000) +
coord_sf(datum = st_crs(fra), expand = TRUE)
On peut ainsi procéder pour tous les DOM.
Plutôt que d’utiliser à chaque fois la valeur ymin
de la carte de la métropole, on peut utiliser à la place la valeur ymax
de la carte affichée dessous pour empiler les cartes les unes au-dessus des autres.
# Définition des positions finales des cartes et modification de la géométrie
# Les valeurs à ajouter ou retirer des coordonnées brutes sont trouvées en tâtonnant.
# On utilise la valeur 'ymax' de la nouvelle géométrie 'mtq2' comme point de départ
# Guadeloupe
glp2_position <- c(st_bbox(fra)$xmin - st_bbox(glp)$xmin - 150000, # position X
st_bbox(mtq2)$ymax - st_bbox(glp)$ymin + 1.25 * 130000) # position Y
glp2 <- glp %>%
mutate(geometry = place_geometry(geometry = st_geometry(glp),
position = glp2_position,
scale = 2.5))
# Réunion
reu2_position <- c(st_bbox(fra)$xmin - st_bbox(reu)$xmin - 150000, # position X
st_bbox(glp2)$ymax - st_bbox(reu)$ymin + 1.25 * 130000) # position Y
reu2 <- reu %>%
mutate(geometry = place_geometry(geometry = st_geometry(reu),
position = reu2_position,
scale = 2.5))
# Mayotte
myt2_position <- c(st_bbox(fra)$xmin - st_bbox(myt)$xmin - 150000, # position X
st_bbox(reu2)$ymax - st_bbox(myt)$ymin + 1.25 * 130000) # position Y
myt2 <- myt %>%
mutate(geometry = place_geometry(geometry = st_geometry(myt),
position = myt2_position,
scale = 3))
# Guyane
guy2_position <- c(st_bbox(fra)$xmin - st_bbox(guy)$xmin - 260000, # position X
st_bbox(myt2)$ymax - st_bbox(guy)$ymin + 1.3 * 50000) # position Y
guy2 <- guy %>%
mutate(geometry = place_geometry(geometry = st_geometry(guy),
position = guy2_position,
scale = 0.4))
Carte finale :
ggplot() +
geom_sf(data = fra) +
geom_sf(data = mtq2) +
geom_sf_text(data = mtq2, aes(label = NAME_0), nudge_y = (st_bbox(mtq2)[4] - st_bbox(mtq2)[2]) / 2 + 50000) +
geom_sf(data = glp2) +
geom_sf_text(data = glp2, aes(label = NAME_0), nudge_y = (st_bbox(glp2)[4] - st_bbox(glp2)[2]) / 2 + 50000) +
geom_sf(data = reu2) +
geom_sf_text(data = reu2, aes(label = NAME_0), nudge_y = (st_bbox(reu2)[4] - st_bbox(reu2)[2]) / 2 + 50000) +
geom_sf(data = myt2) +
geom_sf_text(data = myt2, aes(label = NAME_0), nudge_y = (st_bbox(myt2)[4] - st_bbox(myt2)[2]) / 2 + 50000) +
geom_sf(data = guy2)+
geom_sf_text(data = guy2, aes(label = NAME_0), nudge_y = (st_bbox(guy2)[4] - st_bbox(guy2)[2]) / 2 + 50000) +
coord_sf(xlim = c(st_bbox(guy2)$xmax - 300000, st_bbox(fra)$xmax))
Là où cette méthode est intéressante, c’est qu’on peut maintenant afficher les DOM avec leur vraie échelle.
# Recopie du code précédent, avec une échelle de 1
#Martinique
mtq3_position <- c(st_bbox(fra)$xmin - st_bbox(mtq)$xmin - 150000, # position X
st_bbox(fra)$ymin - st_bbox(mtq)$ymin + 100000) # # position Y
mtq3 <- mtq %>%
mutate(geometry = place_geometry(geometry = st_geometry(mtq),
position = mtq3_position,
scale = 1))
# Guadeloupe
glp3_position <- c(st_bbox(fra)$xmin - st_bbox(glp)$xmin - 150000, # position X
st_bbox(mtq3)$ymax - st_bbox(glp)$ymin + 1.25 * 130000) # position Y
glp3 <- glp %>%
mutate(geometry = place_geometry(geometry = st_geometry(glp),
position = glp3_position,
scale = 1))
# Réunion
reu3_position <- c(st_bbox(fra)$xmin - st_bbox(reu)$xmin - 150000, # position X
st_bbox(glp3)$ymax - st_bbox(reu)$ymin + 1.25 * 130000) # position Y
reu3 <- reu %>%
mutate(geometry = place_geometry(geometry = st_geometry(reu),
position = reu3_position,
scale = 1))
# Mayotte
myt3_position <- c(st_bbox(fra)$xmin - st_bbox(myt)$xmin - 150000, # position X
st_bbox(reu3)$ymax - st_bbox(myt)$ymin + 1.25 * 130000) # position Y
myt3 <- myt %>%
mutate(geometry = place_geometry(geometry = st_geometry(myt),
position = myt3_position,
scale = 1))
# Guyane
guy3_position <- c(st_bbox(fra)$xmin - st_bbox(guy)$xmin - 300000, # position X
st_bbox(myt3)$ymax - st_bbox(guy)$ymin + 1.3 * 130000) # position Y
guy3 <- guy %>%
mutate(geometry = place_geometry(geometry = st_geometry(guy),
position = guy3_position,
scale = 1))
Carte finale, avec les départements à leur échelle :
ggplot() +
geom_sf(data = fra) +
geom_sf(data = mtq3) +
geom_sf_text(data = mtq3, aes(label = NAME_0), nudge_y = (st_bbox(mtq3)[4] - st_bbox(mtq3)[2]) / 2 + 50000) +
geom_sf(data = glp3) +
geom_sf_text(data = glp3, aes(label = NAME_0), nudge_y = (st_bbox(glp3)[4] - st_bbox(glp3)[2]) / 2 + 50000) +
geom_sf(data = reu3) +
geom_sf_text(data = reu3, aes(label = NAME_0), nudge_y = (st_bbox(reu3)[4] - st_bbox(reu3)[2]) / 2 + 50000) +
geom_sf(data = myt3) +
geom_sf_text(data = myt3, aes(label = NAME_0), nudge_y = (st_bbox(myt3)[4] - st_bbox(myt3)[2]) / 2 + 50000) +
geom_sf(data = guy3)+
geom_sf_text(data = guy3, aes(label = NAME_0), nudge_y = (st_bbox(guy3)[4] - st_bbox(guy3)[2]) / 2 + 50000)
Cela permet de voir que la Guyane est de loin le plus grand département français.
sf
Afin de faciliter les opérations d’affichage, on peut créer un nouvel objet sf
regroupant toutes ces données. On affichera les cartes des DOM voulues en filtrant (ou pas) ce nouveau dataframe.
Regrouper les cartes ensemble, en ne gardant que les variables d’intérêts :
# Changer le nom des colonnes pour correspondre au dataframe 'fra'
# NAME_0 -> NAME_2
# GID_O -> HASC_2
list(mtq, glp, reu, myt, guy) %>%
purrr::map(.f = rename, HASC_2 = GID_0, NAME_2 = NAME_0) %>%
setNames(nm = c("mtq", "glp", "reu", "myt", "guy")) %>%
list2env(envir = .GlobalEnv)
## <environment: R_GlobalEnv>
# Regrouper tous les objets 'sf' en un seul
all_sf <- plyr::rbind.fill(fra, mtq, glp, reu, myt, guy)
tail(all_sf)
L’objet obtenu est de classe matrix
.
On peut ajouter d’autres informations, en complétant les données manquantes.
(les départments d’outre-mer sont également des régions ; pour simplifier l’affichage, on va regrouper le code région GID_1
des DOM en un seul code “FRA.14_1”).
all_sf$GID_0 <- "FRA"
all_sf$NAME_0 <- "France"
all_sf$GID_1[is.na(all_sf$GID_1)] <- "FRA.14_1"
all_sf$TYPE_2 <- "Département"
all_sf$CC_2[is.na(all_sf$CC_2)] <- c(971, 972, 974, 976, 973)
tail(all_sf)
On peut maintenant créer de nouveaux dataframes, en modifiant la géométrie (nouvelle position) selon que l’on veut une carte avec les DOM à l’échelle ou non.
# Cartes à l'échelle : modifier les cartes avec les géométries voulues
all_sf_scale <- all_sf
all_sf_scale$geometry[all_sf_scale$NAME_2 == "Martinique"] <- mtq3$geometry
all_sf_scale$geometry[all_sf_scale$NAME_2 == "Guadeloupe"] <- glp3$geometry
all_sf_scale$geometry[all_sf_scale$NAME_2 == "La Réunion"] <- reu3$geometry
all_sf_scale$geometry[all_sf_scale$NAME_2 == "Mayotte"] <- myt3$geometry
all_sf_scale$geometry[all_sf_scale$NAME_2 == "Guyane"] <- guy3$geometry
# Transformer en objet 'sf'
all_sf_scale <- st_as_sf(all_sf_scale)
Visuslisation de la carte des départements avec les départements d’outre-mer à l’échalle.
all_sf_scale %>%
ggplot() +
geom_sf(aes(fill = GID_1), show.legend = FALSE, col = "white", size = 0.2) +
scale_fill_viridis(discrete = TRUE, option = "D") +
geom_sf_text(data = all_sf_scale %>% filter(NAME_2 %in% c("Martinique", "Guadeloupe", "La Réunion", "Mayotte")), aes(label = NAME_2), nudge_y = 80000) +
geom_sf_text(data = all_sf_scale %>% filter(NAME_2 %in% c("Guyane")), aes(label = NAME_2), nudge_y = 250000) +
coord_sf(datum = NA, expand = FALSE) +
labs(title = "Carte des départemens français",
subtitle = "(échelle des DOM respectée)") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
panel.background = element_blank(),
axis.title = element_blank())
Puis la carte avec les échelles des DOM non respectée.
# Cartes avec échelle non respectée : modifier les cartes avec les géométries voulues
all_sf_unscale <- all_sf
all_sf_unscale$geometry[all_sf_unscale$NAME_2 == "Martinique"] <- mtq2$geometry
all_sf_unscale$geometry[all_sf_unscale$NAME_2 == "Guadeloupe"] <- glp2$geometry
all_sf_unscale$geometry[all_sf_unscale$NAME_2 == "La Réunion"] <- reu2$geometry
all_sf_unscale$geometry[all_sf_unscale$NAME_2 == "Mayotte"] <- myt2$geometry
all_sf_unscale$geometry[all_sf_unscale$NAME_2 == "Guyane"] <- guy2$geometry
# Transformer en objet 'sf'
all_sf_unscale <- st_as_sf(all_sf_unscale)
Visuslisation de la carte des départements avec l’échelle des départements d’outre-mer non respectée.
all_sf_unscale %>%
ggplot() +
geom_sf(aes(fill = GID_1), show.legend = FALSE, col = "white", size = 0.2) +
scale_fill_viridis(discrete = TRUE, option = "D") +
geom_sf_text(data = all_sf_unscale %>% filter(NAME_2 %in% c("Martinique", "Guadeloupe", "La Réunion", "Mayotte", "Guyane")), aes(label = NAME_2), nudge_y = 125000) +
coord_sf(datum = NA) +
scale_x_continuous(expand = expand_scale(mult = c(0.2, 0))) +
labs(title = "Carte des départemens français",
subtitle = "(échelle des DOM non respectée)") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
panel.background = element_blank(),
axis.title = element_blank())
# sauvegarder les 2 nouveaux dataframes pour réutilisation
save(all_sf_scale, all_sf_unscale, file = "carte_dep_fra.RData")