IdentifiantMot de passe
Loading...
Mot de passe oublié ?Je m'inscris ! (gratuit)

R et espace

Apprendre le traitement de l'information géographique


précédentsommairesuivant

X. CHAPITRE 10 Introduction aux objets spatiaux et à la cartographie

Objectifs
Ce chapitre présente les principales techniques de cartographie d'objets vectoriels (points et polygones) et d'images raster. Il mobilise les fonctions graphiques de base, mais aussi le package ggplot2 introduit dans le chapitre précédent.

Avertissement
Compte tenu de la contrainte d'impression en noir et blanc, le code utilisé est reproduit, mais les sorties graphiques ne sont pas toujours affichées. C'est également cette contrainte qui amène, dans l'ensemble du manuel, à ne pas toujours respecter les règles de sémiologie graphique concernant les couleurs (cf. Section 9.2Gestion des couleurs).

Prérequis
Notions de cartographie thématique : discrétisation et sémiologie graphique. Connaissance des formats d'objets spatiaux.

Description des packages utilisés
L'importation et la manipulation de données spatiales dans R s'appuient sur un certain nombre de packages spécialisés et de types d'objets spatiaux. Les fonctions et les objets spatiaux de base sont définis dans le package sp (spatial). Tous les autres packages spatiaux dépendent du package sp. Ce dernier, ainsi qu'un certain nombre d'autres packages plus spécialisés, sont développés par Roger Bivand, Edzer Pebesma et Virgilio Gómez-Rubio qui sont également les auteurs du manuel de référence en analyse spatiale avec R.

Avant de commencer cette section, il faut donc installer et charger les packages suivants :

  • sp : la base de tous les autres packages spatiaux, dans lequel sont définis les types d'objets spatiaux ;
  • rgdal : implémentation dans R de la bibliothèque GDAL (Geospatial Data Abstraction Library). Ce package permet d'importer de nombreux formats de données spatiales, raster et vectorielles ;
  • raster : fonctions pour la manipulation des fichiers raster ;
  • rgeos : interface avec la bibliothèque Geometry Engine-Open Source (GEOS) qui permet de manipuler la géométrie des objets spatiaux ;
  • mapproj : fonctions de conversion des coordonnées géographiques en coordonnées projetées ;
  • fields : ensemble de fonctions dédiées aux objets spatiaux, ce package est utilisé ici pour calculer des matrices de distance entre des semis de points ;
  • classInt : permet de discrétiser des variables continues. Il est utilisé ici pour faire des cartes choroplèthes;
  • ggplot2 : package utilisé pour les représentations graphiques, il est aussi utile pour faire des représentations cartographiques (cf. Chapitre 9CHAPITRE 9 Focus sur la visualisation graphique) ;
  • scales : package complémentaire de ggplot2 ;
  • RColorBrewer : implémentation dans R de Color Brewer, un ensemble de palettes de couleurs ;
  • animation : permet de faire des représentations graphiques ou cartographiques animées ;
  • reshape2 : contient un ensemble de fonctions permettant de restructurer l'information d'un tableau (cf. Section 2.4.3Transposition variables-observations).

Description des données utilisées
La distinction entre format vecteur et format raster est essentielle pour le traitement informatique d'images(33). Le format raster, ou bitmap, représente l'image par une matrice de pixels dont les valeurs se traduisent par des couleurs. Les photographies en sont l'exemple le plus courant. Le format vecteur représente des objets géométriques par leurs caractéristiques (points, segments, polygones) et leurs coordonnées.

Dans un usage géographique, le fichier raster ou vectoriel est plus qu'une image, il contient une information géoréférencée. Dans le cas de données raster, cette information est contenue dans les valeurs de la matrice de pixels ; dans le cas de données vectorielles, l'information intéressante peut être la géométrie des objets elle-même (la structure radioconcentrique d'un ensemble de lignes), mais il s'agit souvent des attributs attachés à ces objets (la population d'une commune).

Image non disponible
FIGURE 10.1 - Formats vecteur et raster

Plusieurs types de fichiers sont mobilisés dans cette section pour montrer les manipulations courantes dans un travail de cartographie : introduire des données externes stockées dans un tableau et avoir à gérer différents formats de données spatiales. Ici seuls deux formats courants sont importés, un format ESRI et un format MapInfo, mais la fonction readOGR() permet d'importer tous les formats de la bibliothèque OGR (PostGIS, kml, Idrisi, etc.). Les fichiers nécessaires sont les suivants :

  • popCom3608.csv : fichier des populations dans les communes de Paris et la petite couronne de 1936 à 2008 ;
  • parispc_com.shp : limites communales de Paris et petite couronne, données vectorielles (polygones) au format ArcGIS constitué de quatre fichiers avec les extensions .shp, .dbf, .shx et .prj ;
  • parispc_bnd.shp : contour de l'espace d'étude (Paris et petite couronne), données vectorielles (polygones) au format ArcGIS constitué de quatre fichiers avec les extensions .shp, .dbf, .shx et .prj ;
  • parispc_hop.tab : établissements hospitaliers de Paris et petite couronne, données vectorielles (points) au format MapInfo constitué de deux fichiers avec les extensions .mif, .mid ;
  • Paris_dens.tif : grille de population, données raster au format .tif.

X-A. R pour la cartographie

L'utilisation de R pour la cartographie n'est pas un choix évident. L'utilisateur qui souhaite faire une simple carte choroplèthe aura plus vite fait de la produire avec un logiciel de cartographie classique. L'utilisation de R pour faire de la cartographie se justifie si on envisage l'ensemble du flux de travail (workflow), c'est-à-dire l'intégration de la cartographie dans une chaîne de traitements.

On peut envisager (au moins) quatre situations dans lesquelles R s'avèrera avantageux. D'abord, pour obtenir rapidement des sorties massives : les logiciels de cartographie et de SIG sont des logiciels à interface graphique (même si certains permettent d'automatiser les traitements) et il devient rapidement coûteux de produire manuellement des ensembles de nombreuses cartes.

Ensuite, il sera intéressant d'utiliser R pour la cartographie dans le cadre d'un flux de travail unifié, c'est-à-dire pour intégrer la cartographie dans un flux de travail « classique » (analyse de données classiques, i.e. non spatiales) entièrement fait avec R et se passer ainsi d'importations et d'exportations parfois plus nombreuses que prévues, et qui sont source d'erreurs. En effet, il arrive fréquemment de faire des calculs dans un logiciel de traitement de données, puis d'exporter les résultats pour les cartographier sur un autre logiciel, puis de refaire les calculs suite à une erreur ou une modification, puis de réexporter les résultats à cartographier, etc.

Unifier la chaîne de traitement est d'autant plus utile dans un flux de travail qui mobilise des fonctions d'analyse spatiale implémentées sous R. En effet, R dispose de plusieurs packages d'analyse spatiale et d'analyse de graphes (dont certaines sont abordées dans ce manuel) et il est très pratique de pouvoir manipuler des tableaux de données classiques, des objets spatiaux et/ou des graphes dans un flux de travail unifié.

Enfin, dans certains cas, l'utilisateur cherche à combiner les traitements et l'écriture des contenus : il s'agit d'intégrer des sorties numériques, graphiques et cartographiques dans des supports qui sont également écrits avec R (cf. 9.1). Dans ce cas, la cartographie, comme le reste des traitements, s'intègrera dans la production d'ensemble.

X-B. Prise en main des données spatiales

Certains packages spatiaux renvoient des messages d'avertissement lorsqu'on les charge, en particulier des messages sur les versions des bibliothèques externes utilisées : c'est le cas de rgdal qui, lors du chargement, donne une précision de ce type sur les bibliothèques GDAL et PROJ.4(34).

 
Sélectionnez
1.
2.
3.
4.
5.
6.
library(sp)
library(rgdal)
library(rgeos)
library(mapproj)
library(maptools)
library(raster)

X-B-1. Construction d'un objet spatial

Dans cette première section, un objet spatial est créé à partir de zéro, il s'agira ici de deux sections de route (géométrie linéaire) accompagnées et à la cartographie de deux attributs, la vitesse sur la section et le nombre de voies. Il est très rare qu'un utilisateur ait à créer lui-même des objets spatiaux, la plupart du temps, il récupère les données produites par des organismes tels que l'IGN. Il est cependant important de comprendre comment un tel objet est construit.

Un objet spatial vectoriel est la combinaison d'une composante géométrique et d'une composante sémantique. La composante géométrique comporte des points avec des identifiants et des coordonnées. Ces points peuvent être accompagnés d'un attribut de regroupement s'ils forment des géométries plus compliquées, comme des lignes ou des polygones (cf. Section 9.3.5Variables de regroupement). La composante sémantique est un tableau relié aux géométries par un identifiant et qui contient des attributs correspondant aux objets géométriques.

L'exemple qui suit reproduit une structure de données fréquente pour les données spatiales linéaires telles que les routes ou les voies ferrées. La route est composée d'un ensemble de sections avec des attributs variables (type de route, vitesse, nombre de voies, etc.), chaque section est composée d'un ensemble de lignes, chaque ligne est composée d'un ensemble de points.

Image non disponible
FIGURE 10.2 - Construction d'un objet linéaire

Dans un premier temps, trois tableaux sont créés qui contiennent des points avec des coordonnées.

 
Sélectionnez
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
coordPoints1 <- data.frame(id = c("P1", "P2", "P3"),
                           coordx = c(1, 1, 2),
                           coordy = c(1, 2, 2))

coordPoints2 <- data.frame(id = c("P3", "P4", "P5"),
                           coordx = c(2, 3, 3),
                           coordy = c(2, 2, 1))

coordPoints3 <- data.frame(id = c("P5", "P6", "P7"),
                           coordx = c(3, 4, 5),
                           coordy = c(1, 1, 1))

Dans un deuxième temps, les trois tableaux sont utilisés pour créer trois lignes avec la fonction Line() du package sp. Ces lignes ne contiennent pour le moment rien d'autre qu'un ensemble de points avec leurs coordonnées.

 
Sélectionnez
1.
2.
3.
line1 <- Line(coordPoints1[ , 2: 3])
line2 <- Line(coordPoints2[ , 2: 3])
line3 <- Line(coordPoints3[ , 2: 3])

Dans un troisième temps, les trois lignes sont regroupées dans des sections avec la fonction Lines() et un identifiant leur est attribué. Les deux objets section1 et section2 seront les formes géométriques (features) de l'objet spatial. La géométrie finale de l'objet spatial est créée avec la fonction SpatialLines() qui contient toutes les formes géométriques ainsi qu'un système de projection (ce point est détaillé dans la section suivante).

 
Sélectionnez
1.
2.
3.
4.
5.
6.
7.
8.
section1 <- Lines(slinelist = list(line1, line2),
                  ID = "Sec1")

section2 <- Lines(slinelist = list(line3),
                  ID = "Sec2")

geoLines <- SpatialLines(list(section1, section2),
                         CRS("+init=epsg:2154"))

Enfin, une table attributaire est créée avec deux champs, la vitesse et le nombre de voies. Ce tableau a deux lignes correspondant aux deux formes géométriques. Le lien se fait entre l'identifiant des géométries et les noms de ligne du tableau. La fonction SpatialLinesDataFrame() crée l'objet spatial final qui combine la composante géométrique et la composante sémantique.

 
Sélectionnez
1.
2.
3.
4.
5.
6.
7.
8.
attrTable <- data.frame(SPEED = c(30, 50),
                        LANES = c(2, 3))

row.names(attrTable) <- c("Sec1", "Sec2")

geoDataLines <- SpatialLinesDataFrame(sl = geoLines,
                                      data = attrTable,
                                      match.ID = TRUE)

L'objet spatial sans composante sémantique (geoLines) peut déjà être manipulé et cartographié. L'objet spatial final contient une information sémantique (geoDataLines), il peut être cartographié en utilisant les variables de la table attributaire, en traduisant la variation des attributs numériques par des variations visuelles, de taille ou de couleur par exemple.

 
Sélectionnez
1.
2.
3.
plot(geoDataLines,
     lwd = geoDataLines@data$SPEED,
     axes = TRUE)

Dans cet exemple, un objet spatial de type vectoriel et linéaire a été construit en assemblant des objets géométriques avec des attributs sémantiques. Dans un flux de travail réel, il est plus fréquent de travailler avec des objets spatiaux déjà construits sur lesquels il faut pouvoir extraire et manipuler des éléments géométriques et/ou sémantiques.

Dans les chapitres précédents, on a vu que certains types d'objets, data.frame ou list par exemple, sont composés de plusieurs éléments (slots) auxquels on accède avec l'opérateur $. Les objets spatiaux de type vectoriel comportent deux niveaux d'éléments : on accède au premier avec l'opérateur @ et au second avec l'opérateur $.

Les principaux éléments de premier niveau sont : @data qui contient la table attributaire, @lines (@points ou @polygons selon le type de géométrie) qui contient la géométrie des objets et @proj4string qui contient l'indication du système de projection. Ensuite, il est possible avec l'opérateur $ d'accéder à un champ de la table attributaire de la même façon qu'à un champ de tableau classique. Lors de l'utilisation des opérateurs @ et $, ne pas oublier la touche Tab qui complète automatiquement les noms des objets (cf. Section 1.3Utilisation de RStudio).

 
Sélectionnez
geoDataLines@data$TYPE <- "Highway"
attrTable <- geoDataLines@data

Pour récupérer et manipuler la composante géométrique, plusieurs fonctions sont disponibles, en particulier bbox() qui renvoie l'étendue spatiale des données (bounding box) et coordinates() qui renvoie les coordonnées des objets géométriques. Si l'objet spatial est ponctuel, la fonction renvoie les coordonnées de ces points ; si l'objet est linéaire, elle renvoie les coordonnées des points qui forment les lignes ; si l'objet est de type polygonal, elle renvoie les centroïdes des polygones.

 
Sélectionnez
bboxGeolines <- bbox(geoDataLines)
coordPoints <- coordinates(geoDataLines)

Cette courte introduction a montré comment construire et manipuler un objet spatial à partir d'un exemple très simple. Les sections qui suivent poursuivent l'étude des communes de Paris et de la petite couronne. Deux types de géométries y sont manipulées : des points et des polygones.

X-B-2. Importation et prise en main des objets spatiaux

La fonction readOGR() du package rgdal est utilisée pour charger les données Communes (format ArcGIS) et Hopitaux (format MapInfo). Cette fonction importe les données et les stocke dans un objet de type spatial, selon leur géométrie : SpatialPolygonsDataFrame, SpatialLinesDataFrame ou SpatialPointsDataFrame.

La fonction readOGR() peut lire des formats contenant plusieurs couches, c'est pourquoi il ne suffit pas de préciser le nom du fichier (dsn), mais également le nom de la couche (layer). Dans le cas des fichiers ESRI (shp) et MapInfo (mif), ils ne contiennent qu'une seule couche (du même nom que le fichier) qu'il faut préciser avec l'argument layer.

L'argument (encoding) permet de préciser l'encodage, de la même façon que dans l'importation de simples tableaux (cf. Section 2.2.3Codage des caractères et des séparateurs), et l'argument StringsAsFactors sert à spécifier que les champs alphanumériques ne doivent pas être transformés en champs de type factor. Pour importer l'image raster, la fonction readGDAL() du package raster est utilisée.

 
Sélectionnez
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
muniBound <- readOGR(dsn = "data/parispc_com.shp",
                      layer = "parispc_com",
                      encoding = "utf8",
                      stringsAsFactors = FALSE)

regBoundary <- readOGR(dsn = "data/parispc_bnd.shp",
                      layer = "parispc_bnd",
                      encoding = "utf8",
                      stringsAsFactors = FALSE)

publicHosp <- readOGR(dsn = "data/parispc_hop.MIF",
                      layer = "parispc_hop",
                      encoding = "latin1",
                      stringsAsFactors = FALSE)

densGrid <- readGDAL("data/paris_dens.tif")

Lors de l'importation des fichiers, certaines informations sont renvoyées : le driver utilisé (ESRI et MapInfo dans le cas présent), le nombre d'objets (features) et de variables (fields) pour les fichiers vectoriels, le type de géométries (points et polygones dans cet exemple), la taille de la matrice de pixels dans le cas du fichier raster (392*380 dans ce cas).

Il est intéressant d'examiner le type et la structure de ces données spatiales grâce aux fonctions class() et str(). La fonction str() renvoie les attributs, le type de données spatiales, elle précise l'étendue des données (bbox) et le système de coordonnées (proj4string).

 
Sélectionnez
1.
2.
3.
4.
5.
class(muniBound)
class(publicHosp)

str(muniBound)
str(publicHosp)

Une fois les données importées, elles peuvent être visualisées avec la fonction plot(). Cette fonction générique s'utilise de la même façon pour les données spatiales que pour les données classiques : l'argument col précise la couleur de remplissage, l'argument border la couleur du trait, l'argument pch le type de point, et finalement l'argument add pour superposer chaque nouvelle couche à la précédente.

 
Sélectionnez
1.
2.
3.
4.
5.
6.
plot(muniBound,
     col = "grey",
     border = "white",
     axes = TRUE)

plot(publicHosp, pch = 20, col = "black", add = TRUE)

Image non disponible

Pour visualiser le raster, on commence par créer une palette de couleurs (cf. Section 9.2Gestion des couleurs) qui est donnée comme argument de la fonction image() du package raster. Pour travailler directement sur les valeurs des pixels, il est possible de transformer l'objet raster (SpatialGridDataFrame) en une matrice de pixels avec la fonction as.matrix().

 
Sélectionnez
1.
2.
3.
4.
5.
colPalHeat <- c("white", rev(heat.colors(12)))
image(densGrid, col = colPalHeat)
plot(muniBound, add = TRUE)

pixGrid <- as.matrix(densGrid)

La conception de l'échelle et de la légende est détaillée dans la suite du chapitre.

X-B-3. Le système de projection

Jusqu'à présent, les données spatiales ont été importées et visualisées sans se soucier de leur système de projection qui est pourtant essentiel lorsqu'on manipule des données spatiales. Le système de projection est la transformation mathématique qui assure le passage d'une portion d'espace terrestre en trois dimensions à une représentation en deux dimensions. Cette transformation se traduit nécessairement par une déformation, ce qui a un impact en termes de visualisation (cartographie) et en termes de calculs de surfaces ou de distances (statistique spatiale).

Il existe un ensemble de normes d'identification des systèmes de projection. Le code qui identifie le système de projection suit un standard international (SRID - Spatial Reference System Identifier) et le standard le plus répandu est le code EPSG (European Petroleum Survey Group). Pour connaître le système de projection des objets spatiaux, la fonction proj4string() est utilisée, elle renvoie une information contenue dans les formats de fichiers classiques. Par exemple, pour un format ESRI shape, la projection est donnée dans le fichier d'extension .prj.

 
Sélectionnez
proj4string(muniBound)

Le résultat obtenu est une chaîne de caractères contenant plusieurs informations, principalement le système de projection et l'étendue des données (bounding box). Ici, le système est de type Lambert Conformal Conic (lcc) et l'étendue est donnée par les coordonnées géographiques de deux points définissant la « boîte ». Il est utile de trouver le code qui identifie cette projection dans la liste de codes EPSG.

Certains sites web proposent une telle fonctionnalité(35) en récupérant l'information contenue dans le fichier d'extension .prj du fichier ESRI (shapefile) d'origine. Ceci dit, l'information et la modification du système de projection peuvent se faire intégralement avec R. La projection de l'objet muniBound semble mal référencée et sera convertie en « RGF93- Lambert93 » qui est le système de référence proposé par l'IGN(36).

Pour repérer le code correspondant à cette projection, il faut utiliser la fonction make_EPSG() du package rgdal qui renvoie un data.frame contenant l'ensemble des 4211 systèmes de projection avec leurs identifiants. Une sélection est effectuée sur ce tableau grâce à la fonction grep(), qui recherche les expressions régulières, c'est-à-dire les chaînes de caractères répétés (motifs). En recherchant tous les systèmes qui contiennent la chaîne de caractères « RGF » (Référentiel Géodésique Français), on obtient un tableau de 14 lignes qui contient le code EPSG correspondant au système voulu : il s'agit du code 2154.

 
Sélectionnez
projCodes <- make_EPSG()
selecRgf <- projCodes[ grep("RGF", projCodes$note), ]

Il s'agit maintenant de définir le système de projection des données spatiales, limites communales et hôpitaux, grâce à la fonction spTransform() qui ajoute au système initial le nouveau code EPSG. Avec la fonction proj4string() utilisée plus haut, on vérifie que le code EPSG a bien été modifié.

 
Sélectionnez
1.
2.
3.
4.
5.
muniBound <- spTransform(x = muniBound,
                         CRSobj = CRS("+init=epsg:2154"))

publicHosp <- spTransform(x = publicHosp,
                         CRSobj = CRS("+init=epsg:2154"))

X-B-4. Exportations

Au cours du travail cartographique, on peut souhaiter exporter deux types d'objets : des résultats cartographiques (des images) ou bien des objets spatiaux. Les dispositifs graphiques pour exporter les résultats sous forme d'images sont présentés dans la Section 9.1Exportation des tableaux et des graphiques.

Pour exporter des objets spatiaux, on peut utiliser la fonction writeOGR(), pendant de la fonction readOGR() utilisée précédemment, et créer des fichiers dans les formats les plus courants, ESRI ou MapInfo par exemple. Le site de la bibliothèque GDAL/OGR(37) recense tous les formats disponibles. Dans l'exemple qui suit, le fichier de communes est exporté au format ESRI (.shp) dans un dossier intitulé « MyFolder ».

 
Sélectionnez
1.
2.
3.
4.
writeOGR(muniBound,
         dsn = "MyFolder",
         layer = "muniBound",
         driver = "ESRI Shapefile")

Si le dossier n'existe pas, il est créé au moment de l'exportation. Si seul le nom du dossier est donné sans préciser le chemin qui y mène, ce dossier sera créé dans le répertoire de travail (cf. Section 2.2.1Le répertoire de travail).

X-C. Cartographie des objets ponctuels

Le fichier des hôpitaux est doté d'une variable quantitative de stock intitulée Capacite. Il s'agit de la capacité totale de l'hôpital mesurée en nombre de lits. Cette variable peut être traduite par une variable visuelle de taille pour produire une carte en cercles proportionnels.

La fonction plot() est utilisée en faisant varier la taille des symboles (cex, character expansion factor) selon la capacité des hôpitaux. Par défaut, la variation de la valeur numérique est traduite par une variation proportionnelle du rayon du cercle. Ce paramétrage est corrigé : pour une visualisation appropriée, la proportionnalité doit affecter l'aire du cercle et non son rayon.

Il faut également ajouter une légende ainsi qu'une échelle. La légende est créée puis localisée sur le graphique soit par une localisation relative, soit en précisant les coordonnées x et y. Dans le premier cas, on indique la localisation par des arguments prédéfinis (left, right, topleft, etc.). Dans le second cas, on peut s'aider de la fonction locator() qui permet de récupérer les coordonnées d'un point sur le graphique en cliquant dessus.

Finalement, on utilise la fonction legend() en veillant à reprendre le style exact des symboles des hôpitaux choisis précédemment, notamment la forme, la couleur et la taille. Les arguments pch, col et pt.cex sont utilisées pour spécifier respectivement le type du point, la couleur et la taille des points, et l'argument bty (box type) permet de dessiner ou non un encadré de la légende.

 
Sélectionnez
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
plot(muniBound, border = "darkgrey")

plot(publicHosp,
     pch = 16,
     cex = 0.2 * sqrt(publicHosp$Capacite / pi),
     col = "black", add = TRUE)

title("Capacité des hôpitaux à Paris-Petite couronne")

legendHosp <- c(100, 400, 700, 1000)

legend("bottomleft",
       legend = legendHosp,
       pch = 16, col = "black",
       pt.cex = 0.2 * sqrt(legendHosp / pi),
       bty = "n", title = "Nombre de lits")

Image non disponible

Pour ajouter une échelle, il y a deux solutions, qui demandent toutes deux de spécifier manuellement la localisation du trait, la longueur du trait, l'épaisseur et le texte. Dans le cas ci-dessous sont utilisées les fonctions arrows() pour le trait de l'échelle, puis text() pour le texte de l'échelle. Il est également possible de placer une échelle avec les fonctions locator() et SpatialPolygonsRescale(). La fonction locator() est interactive, elle demande de cliquer sur la carte pour obtenir les coordonnées du point où l'on souhaite placer l'échelle. Puis, la fonction SpatialPolygonsRescale() permet de placer l'échelle à partir de ce point.

Finalement, pour visualiser correctement les cercles proportionnels qui se chevauchent, il est nécessaire de spécifier une couleur de bordure qui les différencie (ici le blanc) et de trier les hôpitaux selon leur taille pour que les cercles les plus gros se retrouvent en arrière-plan.

Attention, l'argument col désigne deux choses différentes selon que l'on représente un objet zonal ou un objet ponctuel. Pour les objets zonaux, col désigne la couleur de remplissage des polygones et border la couleur du trait. Pour les objets ponctuels, col désigne la couleur du symbole, donc du trait, et l'argument bg (background) la couleur du fond, donc le remplissage du symbole ponctuel. Le type de point choisi (pch) est le numéro 21, parce qu'on cherche à distinguer la couleur de remplissage de la couleur de bordure. Il faut donc un type de symbole avec bordure, ce qui n'est pas le cas des types de symboles 15 à 20. On peut obtenir la liste des symboles ponctuels dans l'aide (?pch).

Le fonctionnement de la fonction arrows() nécessite une explication. Il s'agit de situer la droite représentant l'échelle par rapport aux marges du graphique. La fonction par() permet de manipuler les paramètres graphiques généraux, l'argument usr est un vecteur de quatre valeurs : x1, x2, y1, y2 qui désignent les bornes de la sortie graphique. La fonction arrows() demande également quatre valeurs, les coordonnées x et y du premier point de la droite et les coordonnées x et y du second point de la droite. Ici le premier point de la droite a pour coordonnées x1+1000 unités (dans l'unité de mesure de l'objet spatial) et y1+1000 unités (désigne le coin en bas à gauche de la sortie graphique). Le second point de la droite a pour coordonnées x1+10000 unités (soit 10 km) et y1+1000 (bien sûr, y reste constant puisque la droite est horizontale). La fonction text() fonctionne de la même façon, en situant le point à partir duquel commence le texte et en précisant le contenu de ce texte (10 km).

 
Sélectionnez
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.
31.
32.
33.
34.
35.
publicHosp <- publicHosp[ order(publicHosp$Capacite,
                                decreasing = TRUE), ]

plot(muniBound, border = "darkgrey")

plot(publicHosp,
     pch = 21,
     cex = 0.2 * sqrt(publicHosp$Capacite / pi),
     col = "white",
     bg = "black",
     add = TRUE)

title("Capacité des hôpitaux à Paris-Petite couronne")

legendHosp <- c(100, 400, 700, 1000)

legend("topleft",
       legend = legendHosp,
       pch = 21,
       col = "white",
       pt.bg = "black",
       pt.cex = 0.2 * sqrt(legendHosp / pi),
       bty = "n",
       title = "Nombre de lits")

arrows(par() $usr[ 1] + 1000,
       par() $usr[ 3] + 1000,
       par() $usr[ 1] + 10000,
       par() $usr[ 3] + 1000,
       lwd = 2, code = 3,
       angle = 90, length = 0.05)

text(par() $usr[ 1] + 5050,
     par() $usr[ 3] + 2700,
     "10 km", cex = 1.2)

Image non disponible

Ce code semble exagérément long pour produire une simple carte de cercles proportionnels. Sa longueur vient du fait que chacun des éléments est paramétré de façon peu automatisée. Cependant, le code peut être encapsulé dans une fonction beaucoup plus simple pour un usage plus rapide. Le même résultat aurait aussi pu être produit avec le package ggplot2 présenté dans le chapitre précédent.

X-D. Cartographie des objets zonaux

À la date de préparation de ce manuel, il existe principalement deux façons de produire des cartes avec R : avec la fonction générique plot() qui permet de représenter directement des objets spatiaux, ou bien avec la fonction ggplot() du package ggplot2 très utile pour tous types de représentations graphiques. Ces deux méthodes sont expliquées ici ainsi que leurs avantages et leurs inconvénients respectifs. L'exemple d'application est la production d'une carte choroplèthe de la densité de population par commune en 2008.

Pour chacune des deux fonctions plot() et ggplot(), deux critères de discrétisation sont testés : en quartiles et selon l'algorithme de Jenks. Dans ce dernier cas, le package ClassInt, est utilisé, comme dans le Chapitre 4CHAPITRE 4 Analyse univariée. L'algorithme de Jenks est très utilisé pour produire des cartes choroplèthes parce qu'il a été conçu dans ce but et qu'il est implémenté dans la plupart des logiciels de cartographie. Cependant, d'autres algorithmes existent pour discrétiser une variable continue selon des seuils dits « naturels » et ils sont souvent plus efficaces. Certains sont implémentés dans le package ClassInt, comme celui de Fisher ou les algorithmes de classification multivariée (k-means, CAH).

X-D-1. Cartographie avec la fonction plot()

On commence par créer des palettes de couleurs à partir du package RColorBrewer. Ce n'est pas la seule façon de construire une palette, mais ce package présente l'avantage de prédéfinir des palettes adaptées à différents types de représentations (cf. Section 9.2Gestion des couleurs). La palette créée est une palette continue avec cinq couleurs de même ton et d'intensité variable (montée de valeurs).

 
Sélectionnez
library(RColorBrewer)
greyPal <- brewer.pal(n = 5, name = "Greys")

Ensuite, une jointure est faite entre les données attributaires de l'objet spatial et les données externes contenues dans le tableau popCom3608. Il y a deux précautions à prendre pour faire la jointure, l'une concerne la variable de jointure, l'autre concerne la fonction de jointure. La variable de jointure, ici le code de la commune, doit être du même type dans les deux tableaux. La fonction class() permet de s'en assurer. Concernant la fonction de jointure, la fonction merge() présentée dans la Section 2.4.1Superposition, concaténation, jointure peut être source d'erreurs. Il faut en effet s'assurer que la table attributaire n'est modifiée ni dans l'ordre des lignes ni dans le nombre de lignes. Si, au cours de la jointure, la table attributaire bouge (l'ordre des lignes change par exemple) la correspondance entre la composante sémantique (i.e. la table attributaire) de la couche à cartographier et sa composante géométrique n'est plus assurée.

La fonction match() garantit cette intégrité de la jointure. Les lignes qui suivent se lisent de la façon suivante : la table attributaire finale est la combinaison de la table attributaire initiale et du tableau externe (popCom3608). Les lignes du tableau externe sont ajoutées quand il y a correspondance entre les identifiants. L'ajout conserve l'ordre et le nombre de lignes du premier argument de la fonction match(), à savoir la table attributaire de l'objet spatial.

 
Sélectionnez
1.
2.
3.
4.
5.
6.
7.
8.
str(muniBound@data$INSEE_COM)
muniBound$INSEE_COM <- as.character(muniBound$INSEE_COM)

muniBound@data <- data.frame(
    muniBound@data,
    popCom3608[ match(muniBound$INSEE_COM,
                      popCom$CODGEO), ]
)

Ce code nécessaire à la jointure est difficilement lisible et il sera fastidieux de le réécrire à chaque jointure. Il serait intéressant de l'encapsuler dans une fonction (cf. Section 3.3Les fonctions) pour pouvoir faire des jointures plus facilement. Cette fonction de jointure sur la table attributaire est nommée AttribJoin(), ses arguments sont le tableau externe (df), l'objet spatial (spdf) et les variables de jointure correspondantes (df.field et spdf.field).

 
Sélectionnez
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
AttribJoin <- function(df, spdf, df.field, spdf.field) {
    if(is.factor(spdf@data[ , spdf.field]) == TRUE) {
        spdf@data[ , spdf.field] <- as.character(
            spdf@data[ , spdf.field]
        )
    }
    if(is.factor(df[ , df.field]) == TRUE) {
        df[ , df.field] <- as.character(df[ , df.field])
    }
    spdf@data <- data.frame(
        spdf@data,
        df[ match(spdf@data[ , spdf.field], df[ , df.field]), ]
    )
    return(spdf)
}

La fonction de jointure commence par tester le type des variables d'entrée : s'il s'agit de facteurs, ils sont transformés en champs alphanumériques. La fonction de jointure est la fonction match() présentée plus haut. La fonction qui vient d'être créée est simple et générique, elle peut être utilisée sur n'importe quel type de données spatiales.

 
Sélectionnez
1.
2.
3.
4.
muniBound <- AttribJoin(df = popCom3608,
                        spdf = muniBound,
                        df.field = "CODGEO",
                        spdf.field = "INSEE_COM")

Une fois la jointure effectuée, on peut travailler directement sur les données attributaires, les discrétiser et les représenter. Un vecteur de seuils (breaks) est créé avec la fonction quantile(). Celle-ci découpe la variable selon les probabilités données en argument avec la fonction seq(), qui correspondent ici à des quintiles. La discrétisation et l'assignation de la palette de couleurs se font dans la même étape, avec la fonction cut(). Cette fonction renvoie un objet de type factor, qui désigne chaque catégorie par un entier (ici de 1 à 5) et un label (étiquette de valeur) qui prend ici la couleur de la palette créée précédemment. L'argument break est utilisé pour préciser les seuils, l'argument labels pour préciser les étiquettes des valeurs, les arguments include.lowest et right pour créer des intervalles fermés à gauche et ouverts à droite de type [valeur1 ; valeur2[.

Finalement la fonction as.character() transforme le facteur en un vecteur alphanumérique qui contient la couleur correspondant à chaque catégorie. Le nouveau champ créé (DENSQ5) dans la table attributaire contiendra donc le code couleur utilisé par la suite avec la fonction plot(). La légende est construite de la même façon, la différence étant dans la fonction levels() qui renvoie la liste de modalités de la variable. Les modalités sont ici les intervalles de discrétisation.

 
Sélectionnez
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
quintBrks <- quantile(muniBound$DENSITE,
                      probs = seq(0, 1, 0.2),
                      names = TRUE)

muniBound$DENSQ5 <- as.character(cut(muniBound$DENSITE,
                                     breaks = quintBrks,
                                     labels = greyPal,
                                     include.lowest = TRUE,
                                     right = FALSE))

legendQ5 <- as.character(levels(cut(muniBound$DENSITE,
                                    breaks = quintBrks,
                                    include.lowest = TRUE,
                                    right = FALSE)))

La représentation cartographique proprement dite se fait comme dans l'exemple de carte en cercles proportionnels de la section précédente. La variable créée plus haut contient des étiquettes de couleur : en précisant cette variable dans l'argument col, elle devient une variable visuelle.

 
Sélectionnez
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
plot(muniBound, col = muniBound$DENSQ5, border = "black")

legend("bottomleft",
       legend = legendQ5,
       bty = "n",
       fill = greyPal,
       cex = 0.8,
       title = "Densité (hab/km²)")

title("Densité de population à Paris et en petite couronne")

Image non disponible

La même manipulation peut être effectuée avec un autre algorithme de discrétisation, celui de Jenks par exemple. Il suffira de substituer le vecteur de seuils (breaks) en quintiles par un autre vecteur de seuils.

 
Sélectionnez
1.
2.
3.
4.
5.
jenksDiscret <- classIntervals(var = muniBound$DENSITE,
                               n = 5,
                               style = "jenks")

jenksBrks <- jenksDiscret$brks

X-D-2. Cartographie avec la fonction ggplot()

La fonction ggplot() du package ggplot2 permet de personnaliser les aspects de la représentation plus finement que la fonction plot(). Elle a aussi certains inconvénients : elle ne prend pas comme argument un objet de type spatial, il faut donc passer par une transformation préalable expliquée par la suite.

L'usage de ggplot() pour la cartographie se justifie quand le résultat graphique attendu est trop complexe pour pouvoir être fait avec plot() et quand les données spatiales ne sont pas trop lourdes. En effet, la transformation des données et la fonction de représentation sont plus lentes que la fonction plot(). Pour le présent jeu de données, composé de 143 communes et arrondissements, cela ne pose pas de problème, mais il est par exemple impossible de faire, sur un ordinateur de bureau, des cartes des 36 000 communes françaises avec ggplot().

En premier lieu, la variable de densité est discrétisée selon la même méthode que dans la section précédente (quintiles). En revanche, la variable est discrétisée dans le tableau externe (popCom3608) et non dans le tableau de données attributaires de l'objet spatial.

 
Sélectionnez
1.
2.
3.
4.
5.
6.
7.
popCom3608$DENSQ5 <- cut(
    popCom3608$DENSITE,
    breaks = quintBrks,
    include.lowest = TRUE,
    right = FALSE,
    labels = c("Q1", "Q2", "Q3", "Q4", "Q5")
)

Pour visualiser les communes avec ggplot(), il faut transformer l'objet spatial (SpatialDataframe) en un objet que ggplot() puisse lire, à savoir un data.frame classique. La fonction fortify() du package ggplot2 effectue ce travail : pour un objet à géométrie polygonale, comme ici le fond communal, l'objet résultant est un tableau contenant les sommets des polygones avec leurs coordonnées (long, lat) ainsi qu'une variable de regroupement qui identifie le polygone avec les sommets le constituant. La représentation graphique dessinera les polygones en regroupant ces sommets, de la même façon que dans l'exemple présenté au chapitre précédent (cf. Section 9.3.5Variables de regroupement).

C'est à cause de cette transformation du format spatial au format « fortifié » qu'il faut nécessairement discrétiser la variable au préalable et non sur l'objet fortifié. En effet, dans cet objet un polygone est représenté par n lignes du tableau, n étant le nombre de sommets le constituant. Si la discrétisation était faite sur ce tableau, la valeur de la variable associée à chaque polygone serait pondérée par le nombre de sommets, ce qui fausserait le résultat.

Dans un deuxième temps, une jointure est effectuée entre l'objet fortifié et le tableau externe contenant la variable discrétisée. À la différence de la jointure de la section précédente, la fonction merge() peut être utilisée ici, à condition de préciser les arguments all et sort. Cette jointure est l'équivalent d'un left join dans le langage SQL.

 
Sélectionnez
1.
2.
3.
4.
5.
6.
7.
8.
fortMuni <- fortify(muniBound, region = "INSEE_COM")

fortMuni <- merge(x = fortMuni,
                  y = popCom3608,
                  by.x = "id",
                  by.y = "CODGEO",
                  all.x = TRUE,
                  sort = FALSE)

La fonction geom_polygon() sert à graphier des objets zonaux : l'argument group désigne l'identifiant des polygones et l'argument fill désigne la couleur de remplissage, contenue dans le champ DENSQ5 créé précédemment. Finalement on modifie l'objet graphique (mapDens) avec la fonction coord_equal() qui assure la conservation de la même échelle x et y et avec la fonction scale_fill_brewer() qui permet de donner un titre à la légende et de désigner une palette de couleurs.

 
Sélectionnez
1.
2.
3.
4.
5.
6.
7.
8.
9.
mapDens <- ggplot(fortMuni) +
 geom_polygon(aes(x = long,
                  y = lat,
                  group = group,
                  fill = DENSQ5),
              colour = "black") +
 scale_fill_manual(name = "Densité 2008",
                   values = greyPal) +
 coord_equal()

L'objet mapDens peut maintenant être visualisé, dans ce cas avec le thème défini par défaut theme_bw().

 
Sélectionnez
mapDens + theme_bw()

Image non disponible

Il est également possible de définir un thème vide, sans grille, sans axes, sans titres d'axes, etc.

 
Sélectionnez
1.
2.
3.
4.
5.
6.
7.
8.
9.
theme_update(axis.ticks = element_blank(),
             axis.text.x = element_blank(),
             axis.title.x = element_blank(),
             axis.text.y = element_blank(),
             axis.title.y = element_blank(),
             panel.grid.minor = element_blank(),
             panel.grid.major = element_blank(),
             panel.background = element_rect(),
             legend.position = "bottom")

X-E. Cartes take away

Il faut maintenant revenir aux constats posés en début de chapitre (cf. Section 10.1R pour la cartographie). Si l'on met dans la balance la longueur du code à écrire pour produire une simple carte, l'utilisation de R n'est pas toujours avantageuse vis-à-vis de logiciels à interface graphique. R est avantageux si on replace la sortie cartographique dans l'ensemble du flux de travail, et en particulier lorsqu'il faut produire des cartes à la chaîne.

Il reste cependant possible de produire rapidement une carte prête à l'emploi, prête à être intégrée dans un document par exemple. Le package rCarto présenté ci-dessous(38) travaille avec un fichier spatial au format ESRI (.shp) et un tableau (data.frame) existant dans l'espace de travail.

Deux fonctions permettent de reproduire les cartes présentées dans ce chapitre, mapCircles() pour la carte en cercles proportionnels et mapChoropleth() pour la carte choroplèthe. Voici en exemple qui reproduit la carte de la densité de population discrétisée en quintiles.

 
Sélectionnez
1.
2.
3.
4.
5.
6.
7.
8.
9.
library(rCarto)
mapChoropleth(
  shpFile = "data/parispc_com.shp", shpId = "INSEE_COM",
  df = popCom3608, dfId = "CODGEO", var = "DENSITE",
  nclass = 5, palCol = "Greys", style = "quantile",
  posLeg = "topleft", lgdRnd = 0, legend = "Hab./km2",
  title = "Densité de population (2008)",
  author = "Groupe ElementR",
  sources = "Insee, Recensement de la population 2008")

Image non disponible

X-F. Cartographie à la chaîne

L'un des intérêts majeurs de faire de la cartographie avec R est de pouvoir faire des traitements automatisant la production de cartes, ce qui est coûteux à faire manuellement sur un logiciel à interface graphique. Il s'agit ici de comparer les densités à plusieurs époques (de 1936 à 2008), ce qui pose deux questions : la question du format de sortie et la question de la méthode de discrétisation. Concernant la première question, il y a au moins trois façons d'intégrer dans une seule sortie un ensemble de plusieurs cartes :

  • afficher toutes les cartes sur une même page ;
  • afficher toutes les cartes sur plusieurs pages d'un même document ;
  • intégrer toutes les cartes dans une animation qui affiche chaque carte successivement.

L'exemple suivant montre la marche à suivre pour la première option. La deuxième option se fait avec les fonctions pdf(), png() ou autres suivant le format souhaité (cf. Section 9.1Exportation des tableaux et des graphiques). Ces fonctions prendraient comme argument une boucle qui cartographie les données à chaque date du recensement.

La troisième option nécessite l'utilisation de packages spécifiques, par exemple animation qui permet d'intégrer toutes les cartes dans un format animé tel que HTML (fonction saveHTML() ), GIF (fonction saveGIF() ) ou autre.

Concernant la discrétisation, il faut choisir une méthode qui découpe la série selon des valeurs de centralité et de dispersion (moyenne, écart-type, quantiles, etc.) et non une méthode comme celle de Jenks qui ne s'appuie pas sur des résumés numériques fixes. On peut dès lors discrétiser chaque variable de densité de 1936 à 2008 séparément, ou bien discrétiser sur l'ensemble des valeurs de 1936 à 2008. C'est cette seconde option qui est montrée ici.

Concernant la fonction de représentation graphique, cette planche sera construite avec ggplot(), comme dans l'exemple présenté à la Section 9.3.4Construction d'une planche.

Dans un premier temps, les densités de population de 1936 à 2008 sont calculées. Ensuite, ces variables sont discrétisées ensemble. Il n'est pas nécessaire de joindre toutes ces données dans un seul champ pour calculer les quantiles : la fonction quantile() peut être appliquée à l'ensemble du tableau, à condition de le transformer en objet de type matrix, qui est un vecteur bidimensionnel.

 
Sélectionnez
1.
2.
3.
4.
evolDens <- popCom3608[ , 3: 11] / popCom3608$SURF
evolQuint <- quantile(as.matrix(evolDens),
                      names = TRUE,
                      probs = seq(0, 1, 0.2))

Une fois créé le vecteur de seuils, chacune des colonnes du tableau est successivement discrétisée grâce à la fonction cut() appliquée colonne à colonne avec la fonction apply() (cf. Section 3.4L'application de fonctions sur des ensembles). Le résultat final est transformé en tableau et une jointure est effectuée entre ce tableau et l'objet fortifié.

 
Sélectionnez
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
discrDens <- apply(evolDens, 2, cut,
                  breaks = evolQuint,
                  include.lowest = TRUE,
                  right = FALSE,
                  labels = c("Q1", "Q2", "Q3", "Q4", "Q5"))

discrDens <- as.data.frame(discrDens)
discrDens$CODGEO <- popCom3608$CODGEO

colnames(discrDens) <- c(
  "DENS1936", "DENS1954", "DENS1962", "DENS1968", "DENS1975",
  "DENS1982", "DENS1990", "DENS1999", "DENS2008", "CODGEO")

fortMuni <- merge(x = fortMuni,
                  y = discrDens,
                  by.x = "id",
                  by.y = "CODGEO",
                  all.x = TRUE,
                  sort = FALSE)

L'objet fortMuni est un tableau « large » contenant neuf champs de densité, un champ pour chaque année du recensement. Pour construire une planche en utilisant une variable désignant des sous-groupes de données (ici des années), il faudrait en faire une tableau « long » qui n'aurait que deux champs à cartographier : un champ contenant une variable qualitative indiquant l'année du recensement et un champ contenant la variable à représenter (ici des classes de densités). L'objet est transformé en format long avec la fonction melt() (cf. Section 2.4.3Transposition variables-observations).

 
Sélectionnez
1.
2.
3.
4.
5.
fortMuniLong <- melt(
    data = fortMuni,
    id.vars = colnames(fortMuni)[ 1: 7],
    measure.vars = colnames(fortMuni)[ 22: 30]
)

Cet objet peut être cartographié, en utilisant l'argument facet_wrap pour intégrer l'ensemble des cartes dans une seule sortie. Cette fonction constuit une planche de n cartes, une pour chaque modalité de la variable donnée comme argument (facets).

 
Sélectionnez
1.
2.
3.
4.
5.
6.
7.
8.
9.
mapEvolDens <- ggplot(fortMuniLong) +
geom_polygon(aes(long, lat, group = group, fill = value),
             colour = "black", size = 0.3)

mapEvolDens +
  facet_wrap(facets = ~variable) +
  scale_fill_manual(name = "Densité (quintiles)",
                    values = greyPal) +
  coord_equal()

Image non disponible

Un rapide commentaire thématique s'impose : la population de l'agglomération parisienne, comme c'est le cas dans la plupart des grandes métropoles européennes, a tendu à s'étaler au cours du XXe. D'autres résultats concordants ont été obtenus dans les chapitres précédents : le cœur de l'agglomération, c'est-à-dire Paris intramuros, s'est dépeuplé durant la seconde moitié du XXe ; parallèlement, le peuplement de la proche banlieue s'accélérait.

Au terme de ce chapitre, l'utilisateur sait construire et manipuler des objets spatiaux, il sait travailler sur les tables attributaires et faire des jointures de données externes, il dispose enfin des outils pour construire différents types de cartes thématiques. Ces bases seront nécessaires pour aborder le chapitre suivant sur l'analyse spatiale.


précédentsommairesuivant
Le terme image est employé ici dans son acception du langage courant. Dans le domaine de la géomatique, image est synonyme de raster.
Pour installer le package rgdal sur les systèmes Debian ou dérivés, il faut installer au préalable les bibliothèques de développement GDAL et PROJ.4, soit les paquets libgdal-dev et libproj-dev.
Voir par exemple le site http://prj2epsg.org.
Le site de l'IGN est riche en informations sur la géodésie : http://geodesie.ign.fr/index.php?page=geodesie.
Package développé par Timothée Giraud de l'UMS RIATE.

Copyright 2014 : Groupe ElementR, Framasoft (coll. Framabook). R et espace est placé sous licence CC By-SA 3.0.