6.1 Réalisation de cartes
Les cartes sont des outils de présentation très utiles dès lors que les données à représenter sont liées à des coordonnées géographiques. Des logiciels spécialisés dans le traitement de ce type de données et dans la réalisation de cartes existent : ce sont les logiciels de SIG (Systèmes d’Information Géographique, ou GIS en anglais). Le plus connu est sans aucun doute ArcGIS. Des alternatives Open source existent comme GRASS ou QGIS. Ces outils sont assez complexes et requièrent une formation minimale non négligeable pour pouvoir les utiliser.
Une alternative simple est de réaliser des cartes directement dans R. L’avantage de réaliser une carte avec R est de pouvoir associer des représentations graphiques avec les outils d’analyses statistiques. Il s’agit, à la fois, d’une option viable pour l’utilisateur occasionnel, et un choix qui se justifie dans des contextes où la carte n’est pas le seul élément important du projet (carte à intégrer dans un document R Markdown, ou à inclure dans un contenu dynamique tel un “dashboard” écrit en R et Shiny, par exemple). A titre d’exemple, vous pouvez explorer ces deux applications : COVID-16 tracker et ZIP explorer qui vous montrent que la cartographie et la gestion de données spatialisées dans R est parfaitement réalisable avec un résultat très professionnel.
Il existe de très nombreux packages dans R pour réaliser des cartes, cependant, {sf} est un package relativement récent (2017) dont l’objectif est de regrouper et optimiser les fonctionnalités qui sont proposées dans divers packages plus anciens tels que {sp}, {rgdal} et {rgeos}. Nous utiliserons donc ici principalement {sf}, combiné aux fonctionnalités de cartographie de {ggplot2} par l’intermédiaire de notre fonction habituelle de création de graphiques chart().
D’autres packages seront nécessaires, notamment ceux qui contiennent des fonds de cartes comme {tmap}. Dans R, il existe de nombreux objets différents pour représenter des données spatialisées. {sf} propose donc toute une série de fonction pour les convertir en objets sf que nous privilégierons ici. Ainsi, pour avoir une carte du monde, nous pourrons faire ceci :
# Linking to GEOS 3.8.1, GDAL 3.1.4, PROJ 6.3.1# Fond de carte mondial de {tmap} converti en objet sf
world <- st_as_sf(read("World", package = "tmap"))
class(world)# [1] "sf" "dataframe" "tbl_df" "tbl" "data.frame"Un objet sf est tout simplement un tableau de données classiques de type cas par variables que vous avez l’habitude de manipuler, mais qui contient une colonne décrivant la zone géographique concernée. Dans notre objet world, la colonne geometry décrit des polygones grâce au type sfc_MULTIPOLYGON. Le type de la variable geometry varie en fonction des coordonnées géographiques avec par exemple sfc_POINT, sfc_LINESTRING ou encore sfc_POLYGON.
# Simple feature collection with 177 features and 15 fields
# Geometry type: MULTIPOLYGON
# Dimension: XY
# Bounding box: xmin: -16656120 ymin: -8460601 xmax: 16656120 ymax: 8375779
# CRS: +proj=eck4 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
# # A tibble: 177 x 16
# iso_a3 name sovereignt continent area pop_est pop_est_dens economy
# <fct> <fct> <fct> <fct> [km^2] <dbl> <dbl> <fct>
# 1 AFG Afghan… Afghanistan Asia 65286… 2.84e7 43.5 7. Leas…
# 2 AGO Angola Angola Africa 124670… 1.28e7 10.3 7. Leas…
# 3 ALB Albania Albania Europe 2740… 3.64e6 133. 6. Deve…
# 4 ARE United… United Ara… Asia 7125… 4.80e6 67.3 6. Deve…
# 5 ARG Argent… Argentina South Amer… 273669… 4.09e7 15.0 5. Emer…
# 6 ARM Armenia Armenia Asia 2847… 2.97e6 104. 6. Deve…
# 7 ATA Antarc… Antarctica Antarctica 1225921… 3.80e3 0.000310 6. Deve…
# 8 ATF Fr. S.… France Seven seas… 725… 1.4 e2 0.0193 6. Deve…
# 9 AUS Austra… Australia Oceania 768230… 2.13e7 2.77 2. Deve…
# 10 AUT Austria Austria Europe 8252… 8.21e6 99.5 2. Deve…
# # … with 167 more rows, and 8 more variables: income_grp <fct>,
# # gdp_cap_est <dbl>, life_exp <dbl>, well_being <dbl>, footprint <dbl>,
# # inequality <dbl>, HPI <dbl>, geometry <MULTIPOLYGON [m]>À vous de jouer !
Nous voyons qu’outre les frontières des pays dans la variable geometry, {tmap} nous renvoie aussi une série d’autres données utiles, comme la superficie de chaque pays dans area (utile pour rapporter des observations par unité de surface, par exemple). De votre côté, vous pourrez donc combiner également vos données avec leur localisation dans un tableau sf.
De plus, les objets sf sont compatibles avec les fonctions de manipulation de tableaux de données qui sont familières pour vous comme select(), filter(), *_join(),…
# Tableau réduit en éliminant des colonnes superflues
world %>.%
select(., name, area, pop_est, economy, geometry) -> world2
head(world2)# Simple feature collection with 6 features and 4 fields
# Geometry type: MULTIPOLYGON
# Dimension: XY
# Bounding box: xmin: -6237303 ymin: -6576291 xmax: 6474206 ymax: 5306915
# CRS: +proj=eck4 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
# # A tibble: 6 x 5
# name area pop_est economy geometry
# <fct> [km^2] <dbl> <fct> <MULTIPOLYGON [m]>
# 1 Afghanis… 652860.… 2.84e7 7. Least d… (((5310471 4514216, 5408503 4470114, …
# 2 Angola 1246700.… 1.28e7 7. Least d… (((1531585 -773954.7, 1553843 -871831…
# 3 Albania 27400.… 3.64e6 6. Develop… (((1729835 5215797, 1722474 5178345, …
# 4 United A… 71252.… 4.80e6 6. Develop… (((4675864 3138289, 4691298 3144343, …
# 5 Argentina 2736690.… 4.09e7 5. Emergin… (((-5017766 -6571677, -5088441 -65762…
# 6 Armenia 28470.… 2.97e6 6. Develop… (((3677241 5131625, 3791201 5148878, …Pour réaliser une carte, nous utiliserons chart() comme d’habitude, mais en ne lui indiquant aucune “esthétique” à utiliser à l’aide de formule. Les “aesthetics” dans {ggplot2} sont les variables utilisées pour construire le graphique. Ensuite, nous utiliserons geom_sf() qui sait où trouver les polygones à utiliser pour créer le fond de carte.
Malheureusement, la version de chart() dans la SciViews Box 2021 ne supporte pas encore les objets sf. Le code suivant remédie à cette lacune (ajoutez ce code dans chacun de vos documents avant d’utiliser chart() pour faire des cartes). Il sera naturellement intégré dans la prochaine version du package {chart} et donc, dans la prochaine SciViews Box.
SciViews::R
chart.sf <- function(data, specif = aes(), formula = NULL, mapping = aes(), type = NULL, theme = theme_svmap, env = parent.frame()) {
if (!missing(specif))
if (inherits(specif, "formula")) {
formula <- specif
} else if (inherits(specif, "uneval")) {
mapping <- specif
} else rlang::abort("'specif' must be either a formula or aes()/f_aes()")
if (!is.null(formula)) {
args <- chart:::as_list(match.call())[-1]
args$data <- NULL
args$specif <- NULL
args$formula <- NULL
args$mapping <- NULL
args$type <- NULL
args$auto.labs <- NULL
args$env <- NULL
mapping <- chart:::.rename_aes(chart:::.f_to_aes(formula, args, with.facets = TRUE))
# Special case ~0
if (is.numeric(mapping$x) && mapping$x == 0)
mapping$x <- NULL
}
facets <- mapping$facets
mapping$facets <- NULL
p <- ggplot(data = data, mapping = mapping, environment = env) +
theme()
if (!is_null(facets)) {
if (is_null(rlang::f_lhs(facets))) {
p <- p + facet_wrap(facets)
} else {
p <- p + facet_grid(facets)
}
}
if (inherits(p, "ggplot"))
class(p) <- unique(c("Chart", class(p)))
p
}Ensuite, réaliser une carte du monde est enfantin :
Si nous ne voulons pas le fond grisé, ni les méridiens, nous pouvons utiliser le thème ggsn::blank(). Nous verrons plus loin que le package {ggsn} permet d’autres annotations pour orienter et indiquer l’échelle sur la carte ensuite.
À vous de jouer !
Là où cela devient intéressant, c’est que nous avons la possibilité d’utiliser les autres variables contenues dans le tableau sf pour annoter la carte. Par exemple, nous pouvons remplir les polygones (%fill=%) représentant les différents pays à l’aide d’une couleur en fonction de la variable economy. Une petite astuce toutefois : une formule doit obligatoirement contenir une variable x après le tilde. Or ici, nous ne voulons fournir que fill, mais ni x, ni y dans la formule. Pour cela, nous indiquons ~ 0 qui sera compris par chart() comme aucun x et aucun y. Ensuite, nous pourrons rajouter toutes les autres esthétiques comprises par {ggplot2} à la suite dans la formule sous la forme %aes=% variable, donc ici, %fill=% economy.
À vous de jouer !
Pour présenter de l’information sur une carte, nous aurons donc besoin d’une ou plusieurs variables qualitatives ou quantitatives associées à des coordonnées spatiales. Les données et les coordonnées spatiales seront, en pratique, souvent dans deux tableaux différents. Vous aurez donc à les combiner en un tableau unique sf avant de pouvoir réaliser la carte voulue.
Dans la suite de cette partie, nous verrons successivement comment ajouter une échelle et orienter une carte, comment importer des données géographiques et des fonds de cartes, et ce que sont les CRS ou systèmes de coordonnées géodésiques. Ces éléments sont indispensables pour tracer sa carte correctement. Dans la partie suivante, nous nous attaquerons à la partie plus intéressante qui concerne l’annotation de vos cartes.
6.1.1 Échelle et orientation
Toute carte qui se respecte contient deux éléments importants : une orientation de la carte via l’indication du nord géographique, et une échelle qui permet d’avoir une idée des distances. Nous allons ajouter ces éléments à notre carte grâce au package {ggsn}.
La fonction ggsn::north() oriente le nord sur la carte, tandis que ggsn::scalebar() ajoute une barre d’échelle. Ces deux fonctions prennent comme premier argument, le jeu de données utilisé. Pour ggsn::scalebar() il faut aussi indiquer pas mal de choses pour paramétrer l’échelle, en particulier, la distance et l’unité. Cela demande quelques essais et erreurs avant d’avoir une barre d’échelle correcte. Par exemple pour l’Italie, cela donne :
world %>.%
filter(., name == "Italy") %>.%
chart(data = .) +
geom_sf() +
ggsn::north(., location = "topright", symbol = 3) +
ggsn::scalebar(., dist = 150, dist_unit = "km", location = "bottomleft",
transform = FALSE, height = 0.015,
st.dist = 0.02, st.size = 3, border.size = 0.5)
6.1.2 Importation de “shapefiles”
Les fonds de cartes sont principalement stockés dans des fichiers au format shapefile de ESRI. Comme nous l’avons déjà vu, certains packages R comme {tmap} en proposent déjà directement intégrés à R mais souvent, vous devrez les importer vous-même. Votre première tâche sera donc de trouver des données spatiales qui correspondent à la zone géographique étudiée. Vous pouvez les trouver dans des sites spécialisés comme :
- Géoportail des institutions fédérales belges pour la Belgique
- GEODATA sur eurostat pour l’Europe
À vous de jouer !
La fonction st_read() permet de lire ces données dans R. Attention, les fichiers au format shapefile (dont l’extension est .shp) sont associé à trois autres fichiers indispensables à la réalisation de votre carte (avec les extensions .shx, .dbf et .prj). Si, par exemple, vous avez placé ces fichiers dans le sous-répertoire data avec le nom europe, vous ferez :
Les fichiers au format shapefile sont souvent assez volumineux. La précision de la carte aura une très grande importance sur la taille du fichier. Le tableau ci-dessous met en avant la différence de la taille du fichier en fonction de l’échelle de la carte des pays d’Europe https://ec.europa.eu/eurostat/web/gisco/geodata/reference-data/administrative-units-statistical-units/countries#countries20
Cartes de l’Europe au format shapefile en fonction de l’échelle de la carte.
| Echelle | Taille du fichier |
|---|---|
| 1:1 million | 168 Mo |
| 1:3 million | 42.3 Mo |
| 1:10 million | 9.6 Mo |
| 1:20 million | 4.5 Mo |
| 1:60 million | 0.57 Mo |
6.1.3 Systèmes de coordonnées
La terre n’est pas plate, vous le savez (elle n’est pas tout à fait ronde non plus). Son axe de rotation n’est pas fixe. La gravité y est variable. Nous devons donc projeter des objets tridimensionnels en deux dimensions sur notre carte. Une projection va inévitablement plus ou moins déformer notre visualisation en deux dimensions. Il existe plusieurs projections différentes avec des caractéristiques particulières (respect des distances, des volumes, etc.). Vous devrez sélectionner et utiliser celle qui convient le mieux à votre carte.
La fonction st_crs() vous renvoie toutes les informations lié au système géodésique (Coordinate Reference System, CRS) utilisé. Ces notions sont très complexes et sortent du cadre de ce cours. Nous pouvons néanmoins citer le WGS84 qui est le système utilisé par le GPS que nous connaissons tous. Chaque pays a un système qui lui est propre et qui a évolué au cours du temps. Le tout rend le traitement des données spatiales encore plus complexe.
Afin de clarifier et unifier ces différents systèmes de référence, un code est attribué à chacun d’eux dans le standard EPSG. Par exemple, le EPSG:4326 correspond au système WGS84 cité plus haut et utilisé par les GPS. Lors de la combinaison de différentes données spatiales, il sera indispensable de vous assurer que ces données sont compatibles et utilisent le même CRS. La Belgique a employé et continue d’employer plusieurs systèmes différents. Le site epsg.io en recense huit ! Nous vous conseillons vivement de lire maintenant comprendre les CRS pour tout savoir ou presque sur ces systèmes de coordonnées.
La carte du monde présenté précédemment suit les conventions de la norme EPSG:4326, alias WGS84.
# Coordinate Reference System:
# User input: +proj=eck4 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
# wkt:
# PROJCRS["unknown",
# BASEGEOGCRS["unknown",
# DATUM["World Geodetic System 1984",
# ELLIPSOID["WGS 84",6378137,298.257223563,
# LENGTHUNIT["metre",1]],
# ID["EPSG",6326]],
# PRIMEM["Greenwich",0,
# ANGLEUNIT["degree",0.0174532925199433],
# ID["EPSG",8901]]],
# CONVERSION["unknown",
# METHOD["Eckert IV"],
# PARAMETER["Longitude of natural origin",0,
# ANGLEUNIT["degree",0.0174532925199433],
# ID["EPSG",8802]],
# PARAMETER["False easting",0,
# LENGTHUNIT["metre",1],
# ID["EPSG",8806]],
# PARAMETER["False northing",0,
# LENGTHUNIT["metre",1],
# ID["EPSG",8807]]],
# CS[Cartesian,2],
# AXIS["(E)",east,
# ORDER[1],
# LENGTHUNIT["metre",1,
# ID["EPSG",9001]]],
# AXIS["(N)",north,
# ORDER[2],
# LENGTHUNIT["metre",1,
# ID["EPSG",9001]]]]À vous de jouer !
Intéressons-nous maintenant à l’Afrique australe définie par l’ONU comme constituée des pays suivants : Afrique du Sud, Namibie, Botswana, Lesotho et Swaziland.
Nous commençons avec une carte du continent africain tout entier :
world %>.%
filter(., continent == "Africa") %>.%
mutate(., austral = name %in% southern_africa) -> africa
chart(data = africa, ~0 %fill=% austral) +
geom_sf()
La fonction qui permet de transformer les coordonnées géographiques d’un CRS à un autre est st_transform(). Nous pouvons simplement comparer les mêmes données spatiales avec deux CRS différents. Ce premier graphique utile la convention EPSG:4326 présenté précédemment.
Ce second graphique utilise le CRS EPSG:22235 qui est centré sur l’Afrique du Sud. Nous n’avons aucunement changé les données. Il s’agit simplement de deux représentations différentes des mêmes données géographiques.
africa %>.%
st_transform(., 22235) -> africa_sa
chart(data = africa_sa, ~0 %fill=% austral) +
geom_sf()
La carte dont le système géodésique est adapté à la zones étudiée est plus intéressante, visuellement parlant et aussi parce que les distances et les surfaces dans cette zone sont mieux respectées. En effet, cette zone est moins déformée par la projection en deux dimensions.