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 17 fields
# Attribute-geometry relationships: constant (12), aggregate (3), identity (2)
# Geometry type: MULTIPOLYGON
# Dimension: XY
# Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.645
# Geodetic CRS: WGS 84
# First 10 features:
# iso_a3 name sovereignt continent area pop_est
# 1 AFG Afghanistan Afghanistan Asia 642393.20 [km^2] 38041754
# 2 ALB Albania Albania Europe 28298.42 [km^2] 2854191
# 3 DZA Algeria Algeria Africa 2312302.79 [km^2] 43053054
# 4 AGO Angola Angola Africa 1249442.03 [km^2] 31825295
# 5 ATA Antarctica Antarctica Antarctica 12258324.81 [km^2] 4490
# 6 ARG Argentina Argentina South America 2783123.36 [km^2] 44938712
# 7 ARM Armenia Armenia Asia 29555.06 [km^2] 2957731
# 8 AUS Australia Australia Oceania 7705931.67 [km^2] 25364307
# 9 AUT Austria Austria Europe 83755.95 [km^2] 8877067
# 10 AZE Azerbaijan Azerbaijan Asia 86152.91 [km^2] 10023318
# pop_est_dens economy income_grp gdp_cap_est
# 1 5.921880e+01 7. Least developed region 5. Low income 507.1007
# 2 1.008604e+02 6. Developing region 4. Lower middle income 5353.1806
# 3 1.861912e+01 6. Developing region 3. Upper middle income 3973.9573
# 4 2.547161e+01 7. Least developed region 3. Upper middle income 2790.7047
# 5 3.662817e-04 6. Developing region 2. High income: nonOECD 200000.0000
# 6 1.614686e+01 5. Emerging region: G20 3. Upper middle income 9912.2779
# 7 1.000753e+02 6. Developing region 4. Lower middle income 4622.4623
# 8 3.291530e+00 2. Developed region: nonG7 1. High income: OECD 55060.3255
# 9 1.059873e+02 2. Developed region: nonG7 1. High income: OECD 50137.6187
# 10 1.163433e+02 6. Developing region 3. Upper middle income 4793.5225
# life_exp well_being footprint HPI inequality gender press
# 1 61.982 2.436034 1.139635 16.21067 NA 0.665 19.09
# 2 76.463 5.255482 3.674262 46.33355 29.4 0.116 54.10
# 3 76.377 5.217018 3.110979 47.44608 27.6 0.460 41.98
# 4 NA NA NA NA 51.3 0.520 52.44
# 5 NA NA NA NA NA NA NA
# 6 75.390 5.908279 6.391003 43.79233 40.7 0.292 63.13
# 7 72.043 5.300569 3.429790 43.00195 27.9 0.198 71.60
# 8 84.526 7.111599 17.703905 39.47980 34.3 0.063 73.42
# 9 81.580 7.079641 12.147266 45.34394 30.7 0.048 74.69
# 10 69.366 5.200000 6.102844 34.37222 26.6 0.329 27.99
# geometry
# 1 MULTIPOLYGON (((66.217 37.3...
# 2 MULTIPOLYGON (((20.605 41.0...
# 3 MULTIPOLYGON (((-4.923 24.9...
# 4 MULTIPOLYGON (((12.322 -6.1...
# 5 MULTIPOLYGON (((-61.883 -80...
# 6 MULTIPOLYGON (((-68.634 -52...
# 7 MULTIPOLYGON (((46.483 39.4...
# 8 MULTIPOLYGON (((147.689 -40...
# 9 MULTIPOLYGON (((16.88 48.47...
# 10 MULTIPOLYGON (((46.144 38.7...À 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
# Attribute-geometry relationships: constant (1), aggregate (2), identity (1)
# Geometry type: MULTIPOLYGON
# Dimension: XY
# Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 42.688
# Geodetic CRS: WGS 84
# name area pop_est economy
# 1 Afghanistan 642393.20 [km^2] 38041754 7. Least developed region
# 2 Albania 28298.42 [km^2] 2854191 6. Developing region
# 3 Algeria 2312302.79 [km^2] 43053054 6. Developing region
# 4 Angola 1249442.03 [km^2] 31825295 7. Least developed region
# 5 Antarctica 12258324.81 [km^2] 4490 6. Developing region
# 6 Argentina 2783123.36 [km^2] 44938712 5. Emerging region: G20
# geometry
# 1 MULTIPOLYGON (((66.217 37.3...
# 2 MULTIPOLYGON (((20.605 41.0...
# 3 MULTIPOLYGON (((-4.923 24.9...
# 4 MULTIPOLYGON (((12.322 -6.1...
# 5 MULTIPOLYGON (((-61.883 -80...
# 6 MULTIPOLYGON (((-68.634 -52...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",
# ENSEMBLE["World Geodetic System 1984 ensemble",
# MEMBER["World Geodetic System 1984 (Transit)"],
# MEMBER["World Geodetic System 1984 (G730)"],
# MEMBER["World Geodetic System 1984 (G873)"],
# MEMBER["World Geodetic System 1984 (G1150)"],
# MEMBER["World Geodetic System 1984 (G1674)"],
# MEMBER["World Geodetic System 1984 (G1762)"],
# MEMBER["World Geodetic System 1984 (G2139)"],
# ELLIPSOID["WGS 84",6378137,298.257223563,
# LENGTHUNIT["metre",1]],
# ENSEMBLEACCURACY[2.0]],
# 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["Horizontal component of 3D system."],
# 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.