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, si les besoins sont ponctuels, consiste à réaliser des cartes directement dans R. L’avantage 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 tableau de bord écrit en R et Shiny, ou d’autres). À 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éalisables avec un résultat très professionnel.
Il existe de très nombreux packages dans R pour traiter les données spatialisées et réaliser des cartes. {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 fonctions pour les convertir en objets sf que nous privilégierons ici. Ainsi, pour avoir une carte du monde, nous pourrons faire ceci :
SciViews::R("spatial") # Charge les packages nécessaires
# Fond de carte mondial de {tmap} converti en objet sf
world <- st_as_sf(read("World", package = "tmap"))
class(world)# [1] "sf" "data.frame"Un objet sf est un tableau de données classiques de type cas par variables data.frame 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: -180 ymin: -89.9 xmax: 180 ymax: 83.64513
# Geodetic CRS: WGS 84
# First 10 features:
# iso_a3 name sovereignt continent
# 1 AFG Afghanistan Afghanistan Asia
# 2 AGO Angola Angola Africa
# 3 ALB Albania Albania Europe
# 4 ARE United Arab Emirates United Arab Emirates Asia
# 5 ARG Argentina Argentina South America
# 6 ARM Armenia Armenia Asia
# 7 ATA Antarctica Antarctica Antarctica
# 8 ATF Fr. S. Antarctic Lands France Seven seas (open ocean)
# 9 AUS Australia Australia Oceania
# 10 AUT Austria Austria Europe
# area pop_est pop_est_dens economy
# 1 652860.000 [km^2] 28400000 4.350090e+01 7. Least developed region
# 2 1246700.000 [km^2] 12799293 1.026654e+01 7. Least developed region
# 3 27400.000 [km^2] 3639453 1.328268e+02 6. Developing region
# 4 71252.172 [km^2] 4798491 6.734519e+01 6. Developing region
# 5 2736690.000 [km^2] 40913584 1.495003e+01 5. Emerging region: G20
# 6 28470.000 [km^2] 2967004 1.042151e+02 6. Developing region
# 7 12259213.973 [km^2] 3802 3.101341e-04 6. Developing region
# 8 7257.455 [km^2] 140 1.929051e-02 6. Developing region
# 9 7682300.000 [km^2] 21262641 2.767744e+00 2. Developed region: nonG7
# 10 82523.000 [km^2] 8210281 9.949082e+01 2. Developed region: nonG7
# income_grp gdp_cap_est life_exp well_being footprint inequality
# 1 5. Low income 784.1549 59.668 3.8 0.79 0.42655744
# 2 3. Upper middle income 8617.6635 NA NA NA NA
# 3 4. Lower middle income 5992.6588 77.347 5.5 2.21 0.16513372
# 4 2. High income: nonOECD 38407.9078 NA NA NA NA
# 5 3. Upper middle income 14027.1261 75.927 6.5 3.14 0.16423830
# 6 4. Lower middle income 6326.2469 74.446 4.3 2.23 0.21664810
# 7 2. High income: nonOECD 200000.0000 NA NA NA NA
# 8 2. High income: nonOECD 114285.7143 NA NA NA NA
# 9 1. High income: OECD 37634.0832 82.052 7.2 9.31 0.08067825
# 10 1. High income: OECD 40132.6093 81.004 7.4 6.06 0.07129351
# HPI geometry
# 1 20.22535 MULTIPOLYGON (((61.21082 35...
# 2 NA MULTIPOLYGON (((16.32653 -5...
# 3 36.76687 MULTIPOLYGON (((20.59025 41...
# 4 NA MULTIPOLYGON (((51.57952 24...
# 5 35.19024 MULTIPOLYGON (((-65.5 -55.2...
# 6 25.66642 MULTIPOLYGON (((43.58275 41...
# 7 NA MULTIPOLYGON (((-59.57209 -...
# 8 NA MULTIPOLYGON (((68.935 -48....
# 9 21.22897 MULTIPOLYGON (((145.398 -40...
# 10 30.47822 MULTIPOLYGON (((16.97967 48...À 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 (s)select(), (s)filter(), *_join()…
# Tableau réduit en éliminant des colonnes superflues
world %>.%
sselect(., 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: -73.41544 ymin: -55.25 xmax: 75.15803 ymax: 42.68825
# Geodetic CRS: WGS 84
# name area pop_est economy
# 1 Afghanistan 652860.00 [km^2] 28400000 7. Least developed region
# 2 Angola 1246700.00 [km^2] 12799293 7. Least developed region
# 3 Albania 27400.00 [km^2] 3639453 6. Developing region
# 4 United Arab Emirates 71252.17 [km^2] 4798491 6. Developing region
# 5 Argentina 2736690.00 [km^2] 40913584 5. Emerging region: G20
# 6 Armenia 28470.00 [km^2] 2967004 6. Developing region
# geometry
# 1 MULTIPOLYGON (((61.21082 35...
# 2 MULTIPOLYGON (((16.32653 -5...
# 3 MULTIPOLYGON (((20.59025 41...
# 4 MULTIPOLYGON (((51.57952 24...
# 5 MULTIPOLYGON (((-65.5 -55.2...
# 6 MULTIPOLYGON (((43.58275 41...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 2024 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("spatial")
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
}
# Convenience functions for north arrow and scale
# Equivalent to ggspatial::annotation_north_arrow but with different defaults
annotation_north <- function(mapping = NULL, data = NULL, ...,
location = "tr", height = unit(0.8, "cm"), width = unit(0.8, "cm"),
style = ggspatial::north_arrow_fancy_orienteering) {
ggspatial::annotation_north_arrow(mapping = mapping, data = data, ...,
location = location, height = height, width = width, style = style)
}
annotation_scale <- ggspatial::annotation_scaleEnsuite, réaliser une carte du monde est enfantin :
Si nous ne voulons pas le fond grisé, nous pouvons utiliser theme_minimal(). Si nous ne voulons pas non plus les méridiens, nous utiliserons alors theme_classic().
À 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 à l’annotation de cartes.
6.1.1 Échelle et orientation
Toute carte qui se respecte contient deux éléments importants : une orientation de la carte avec 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 {ggspatial}.
La fonction annotation_north() oriente le nord sur la carte, tandis que annotation_scale() ajoute une barre d’échelle. Ces deux fonctions prennent comme premier argument, le jeu de données utilisé. La fonction annotation_scale() nous averti si les distances varient sur la carte de plus de 10% (lié à la fois à la projection choisie et à la zone géographique représentée). Par exemple pour l’Italie, cela donne :
italy <- filter(world, name == "Italy")
chart(italy) +
geom_sf() +
theme_minimal() +
annotation_north() +
annotation_scale()# Scale on map varies by more than 10%, scale bar may be inaccurate
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és à 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() renvoie toutes les informations liées 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 plusieurs. 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: EPSG:4326
# wkt:
# GEOGCRS["WGS 84",
# DATUM["World Geodetic System 1984",
# ELLIPSOID["WGS 84",6378137,298.257223563,
# LENGTHUNIT["metre",1]]],
# PRIMEM["Greenwich",0,
# ANGLEUNIT["degree",0.0174532925199433]],
# CS[ellipsoidal,2],
# AXIS["geodetic latitude (Lat)",north,
# ORDER[1],
# ANGLEUNIT["degree",0.0174532925199433]],
# AXIS["geodetic longitude (Lon)",east,
# ORDER[2],
# ANGLEUNIT["degree",0.0174532925199433]],
# USAGE[
# SCOPE["unknown"],
# AREA["World"],
# BBOX[-90,-180,90,180]],
# ID["EPSG",4326]]À 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 utilise la convention EPSG:4326 présentée 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 zone é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.