Dans de nombreuses situations, il peut être utile de représenter graphiquement des disparités géographique, par le biais d’une carte. Même si R n’est pas un outil cartographique dédié (comme peut l’être ArcGIS ou QGIS), de nombreux packages existent pour permettre de manipuler des données spatiales et pour fabriquer des cartes.

L’utilisation de R peut s’avérer judicieuse quand il s’agit d’intoduire la production de carte dans un workflow d’analyse reposant déjà sur R. On évite ainsi de multiples importations/exportations des données avec leur lot de problèmes et erreurs, et on profite de la flexibilité du langage R et ses packages.

Ce document a pour but :

Les problématiques liées à l’analyse spatiale sous R ne seront pas traitées ici.

1 Cartographie sous R : théorie

Pour réaliser une carte, on a généralement recours à plusieurs ingrédients :

1.1 Le fond de carte

Un fond de carte permet de représenter graphiquement les caractéristiques d’une région géographique d’intérêt. Le fond de carte permet de “planter le décor” et d’apporter les éléments contextuels utiles. C’est sur ce fond de carte que vont venir se greffer nos données géolocalisées.

Il existe une quantité infinie de fond de carte différent, qui peuvent grossièrement être rangés dans deux catégories :

Raster et vector layer, exemple

A droite, une représentation matricielle de la zone. A gauche, une représentation vectorielle.
Crédits : Pau Fonseca i Casas, CC BY-NC-ND 3.0

  • les représentations vectorielles (ou vector layer). Dans ce cas, les caractéristiques de la région d’intérêt sont représentées sous la forme d’un ensemble de formes (ou features) comme des points ou des tracés (lignes, polygones…). Ce format est particulièrement adapté à la représentation de zones géographiques bien délimitées, telles que les limites administratives (pays, régions, départements, communes), les limites géologiques, l’emplacement de villes, etc…

  • les représentations matricielles (ou raster layer). Dans ce cas, les caractéristiques de la région d’intérêt sont stockées sous la forme d’une image (ensemble de pixel rangés en lignes/colonnes, d’où le nom). A chaque pixel est associé une valeur, permettant de décrire les caractéristiques. Par exemple, l’altitude peut être représentée sous le forme de pixels allant du vert au rouge, les plans d’eau peuvent être représentés par des pixels de couleur bleu, etc…

Bien évidemment, il est possible de combiner plusieurs représentations pour créer un fond de carte plus complexe (ex : combiner un fond raster représentant la densité de médecins au Km² et des données vectorielles représentant les limites des communes de France).

EN PRATIQUE :

Pour obtenir des fonds de cartes reprenant les limites des régions / départements / communes de France (sous forme de vector layer), le plus simple est d’utiliser les données GEOFLA mises à disposition par l’IGN. Ces données sont distribuées sont au format ESRI Shapefile (SHP), standard très largement répandu et quasi-universel. A noter que le système cartographique utilisé dans ces fichiers est le Lambert-93, le format légal en France.

Dans certains cas, ces tracés peuvent ne pas être suffisants. D’autres versions, avec des niveaux de détail plus grossier ou plus fin sont distribuées sur internet, en particulier via la fondation OpenStreetMap, ou OSM pour les intimes. La plupart des jeux de données OSM directement réutilisables (fond de carte et cie) sont publiées sur le site de l’Open Data. A titre d’exemple, les données de tracé administratifs issus d’OSM sont disponibles ici. Attention, les données issues d’OSM sont diffusées dans le système WGS84.

1.2 Données géolocalisées

De nos jours, les données géolocalisées ne sont pas ce qui manque. Les données géolocalisées peuvent être rattachées à une zone particulière (ex : à une région ou un département français), ou bien être rattachées à un point (ex : coordonnées GPS votre domicile).

En fonction du type de données (zone/point), l’exploitation qui peut en être faite va varier. Avec des données de zone, il est possible de créer des cartes à bases d’aplat de couleurs, ou carte choroplèthe. Avec des données de points, on utilisera plutôt des représentations ponctuelles (ex : un point de taille variable, en fonction de la valeur) ou des représentation de densité/intensité (carte de chaleur, ou heatmap).

Pour les données de zone, on a généralement des données au format suivant :

  • une colonne pour l’identifiant de la zone concernée
  • une ou plusieurs colonnes contenant les valeurs des variables d’intérêt (ex : population/m², nombre d’hospitalisations…)

L’enjeux sera alors généralement de fusionner sous R ces informations de zones à un objet décrivant ces zones (lui-même issu d’un fichier SHP le plus souvent). Par exemple, rattacher les informations sur le nombre d’hospitalisation par départements à un object spatial contenant les contours des départements… A noter que dans ce cas là, les problèmes de coordonnées/projection ne se posent qu’à la création de l’objet spatial, pas à l’import des données.

Pour les données ponctuelles, on a généralement les données au format suivant :

  • une colonne longitude
  • une colonne latitude
  • une ou plusieurs colonnes contenant les valeurs des variables d’intérêt

Dans cette situation, il est impératif de connaître le système de coordonnées utilisé pour le calcul des coordonnées géographiques, sans quoi la représentation sera fausse. Le plus souvent, c’est le système WGS84 qui sera utilisé (si on parle de GPS quelque part ou que ce n’est pas précisé, c’est le WGS84…). Si on souhaite fusionner ces données à un objet spatial (raster ou vectoriel), il faudra s’assurer que les systèmes de coordonnées et la projection est la même.

L’enjeux sera ici de créér un objet spatial à partir des coordonnées longitude/latitude, puis d’ajouter/fusionner les données d’intérêt.

1.3 Système de coordonnées et projection cartographique

Imaginons qu’on souhaite expliquer à quelqu’un où se situe précisement Paris sur la Terre. Il y a une quantité infinie de façon de le faire, en fonction du référentiel utilisé. Prenons par exemple comme point d’origine Rio de Janeiro, au Brésil. Dans ce cas là, on peut situer Paris en précisant que la ville est à X km vers l’Est, et à Y km vers le Nord. A partir du point d’origine (Rio de Janeiro) et des coordonnées (X,Y), vous avez la position de Paris. Cette position sous-entend que la position de Rio de Janeiro est déjà connue, et qu’il est admis que les coordonnées sont données en Est puis Nord, en kilomètres.

Mais, si cela peut sembler suffisant de première abord, il existe un grand nombre de paramètres qui ne sont pas pris en compte ici : la dérive des continents (la distance Rio/Paris change constamment), la déformation de la Terre (qui n’est pas une sphère parfaite), l’altitude

En pratique, pour positionner des éléments sur le globe, on a recours au coordonnées à base de latitude/longitude.

Longitude et latitude : exemple

Coordonnées sphériques : latitude δ, longitude θ, and rayon ρ.
Crédits : Wikimedia Commons

Ces coordonnées vont dépendre de plusieurs éléments, les principaux étant :

  • l’ellipsoïde de référence, c’est-à-dire la forme de la sphère qui va représenter la Terre. Comme dit juste au dessus, la terre n’est pas une sphère parfaite : elle présente des irrégularités plus ou moins importantes selon les endroits (en gros, c’est une patatoïde !). Pour les calculs, on utilise des approximations, des formes simplifiées qui s’approche de la forme de la Terre : des ellipsoïdes de révolution.
    En fonction des utilisations et de la représentation voulue, le choix de l’ellipsoïde doit être adapté. Par exemple, pour une représentation globale de la Terre, on choisira l’ellipsoïde globalement le plus proche de la forme générale de la Terre. A l’inverse, pour une représentation d’une zone particulière (ex : la France) on peut utiliser un ellipsoïde plus proche de la forme de la Terre à cet zone, pour des raisons de précisions.
    Parmi les ellipsoïdes fréquemment utilisés, on peut citer le WGS84 (celui utilisé dans le système du même nom) ou l’IAG-GRS80 (celui utilisé dans le système Lambert-93). A noter que, malgré quelques différences subtiles, ces deux ellipsoïdes sont quasi-identiques.

  • un point d’origine. Selon les cas, l’origine sera déterminée soit par la position du centre de l’ellipsoïde, soit par un point de référence à la surface de la terre.

On regroupe l’ensemble de ces éléments sous le terme de système de coordonnées géographique (coordinate reference system ou CRS pour nos amis anglais) ou de datum. Il en existe une très grande quantité.

Pour faire de la cartographie, un paramètre supplémentaire est nécessaire. En effet, une carte est une représentation plane du globe terrestre, il faut donc appliquer une projection aux coordonnées pour pouvoir les afficher dans un plan. Là aussi, de multiples possibilités s’offre à nous : projection cylindrique, projection conique… avec des propriétés distinctes. Le choix de la projection va fortement influencer le rendu sur la carte (ex : taille et forme des pays).

L’ensemble datum + projection est désigné par le terme système de référence spatiale (spatial reference system ou SRS). Ils sont référencés sous la forme d’un code EPSG à 4 chiffres, pour simplifier leur utilisation.

EN PRATIQUE :

Des coordonnées en latitude/longitude ne veulent rien dire si le datum utilisé n’est pas donné. Pour une carte, cette information est nécessairement couplée au type de projection utilisée. Les principaux codes EPSG à connaître sont le EPSG:4326 (basé sur le système WGS84, très usité puisque le GPS se base dessus) et le EPSG:2154 (basé sur le système Lambert-93).
Pour trouver un code EPSG, on pourra se référer aux sites internet spatialrefrence.org ou epsg.io.

2 Cartographie sous R : pratique

Nous allons utiliser les packages sp et rgdal pour importer un fond de carte au format SHP reprenant le tracé des communes de France. Nous travaillerons ensuite sur l’Ile-de-France : des données extérieures sur le temps de parcours pour hospitalisation seront jointes à notre object spatial, et une carte basée sur ces données sera produite.
Les packages RColorBrewer et classInt seront utilisés pour sélectionner des groupes de couleurs adaptées, et pour catégoriser simplement des variables, respectivement.
Le rendu des carte se fera par l’intermédiaire du package base, partie intégrante de la distribution standard de R. D’autres packages, comme ggplot2 , rCarto , maptools ou ggmap (un dérivé de ggplot2) peuvent être utilisés. Des exemples utilisant ces packages seront ajoutés à l’avenir.

Première chose à faire : installer et charger les packages nécessaires

# Installation des packages de base
install.packages("sp")
install.packages("rgdal")

# Installation des packages accessoires
install.packages("RColorBrewer")
install.packages("classInt")

# Installation des packages de rendu
install.packages("ggplot2")
install.packages("ggmap")
install.packages("maptools")
# Import des librairies nécessaires
library("sp")
library("rgdal")

library("RColorBrewer")
library("classInt")

library("ggplot2")
library("ggmap")
library("maptools")

On récupère sur le site de l’IGN la carte et ses données associées (fichier .PROJ contenant les informations de projection et fichier .DBF contenant les informations géolocalisées fournies par l’IGN). L’ensemble des données utilisées est disponible ici dans l’onglet Téléchargements. Nous utiliserons pour cette exemple les données 2013 des communes de France métropolitaine.

Le fichier 7z (format 7zip, lisible par le logiciel du même nom) contient un ensemble de dossiers/sous-dossiers. Les informations qui nous intéressent sont dans un sous-dossier du type DONNEES_LIVRAISON. Décompressez l’ensemble des fichiers de ce sous-dossier. Nous allons enfin pouvoir passer à la fabrication de la carte dans R. De manière général, il est recommandé de vérifier que les données qu’on va manipuler sont correctement reconnues dans R. Ceci s’applique à toutes les situations, et la cartographie ne fait pas exception.

2.1 Chargement des données spatiales

Il existe plusieurs façons de charger des données spatiales dans R, faisant intervenir différents packages. Ici, nous utilserons uniquement les fonctions de chargement issus du packages rgdal. Ce package est une interface entre R et la “Geospatial Abstraction Library” (GDAL), une bibliothèque de fonctions dédiées aux manipulations de données spatiales très largement utilisée.

La fonction ogrInfo permet de scanner les données que l’on veut manipuler, de façon à vérifier qu’elles seront bien lues et reconnues par la suite. Cette fonction ogrInfo (tout comme readOGR, qu’on va voir plus loin) prend 2 arguments indispensables :

  • le chemin pointant vers le fichier SHP et ses données associées. Si vos données sont dans le répertoire courant de travail sous R, il est quand même nécessaire de remplir cet argument, avec “.” (le point . désignant sous R le répertoire de travail courant)
  • le nom du layer, c’est-à-dire le nom du fichier .SHP (sans mettre l’extension au bout). Il est communément admis que l’ensemble des fichiers relatifs au fichier SHP portent le même nom, seule l’extension change. Ex : COMMUNE.SHP est accompagné de COMMUNE.PROJ, de COMMUNE.DBF, etc. Ainsi, avec une commande, R va charger le SHP, reconnaitre la projection et charger les données associées…

On commence par déterminer l’emplacement de notre fichier SHP, et on scan son contenu :

# Emplacement de l'archive décompressée, à remplacer par le votre
pathToShp <- "./data/shp"
# Description des données via orgInfo
# Attention à ne pas mettre l'extension à la fin du nom
ogrInfo(dsn = pathToShp,layer="COMMUNE")
## Source: "./data/shp", layer: "COMMUNE"
## Driver: ESRI Shapefile number of rows 36594 
## Feature type: wkbPolygon with 2 dimensions
## Extent: (99226 6049647) - (1242375 7110480)
## CRS: +proj=lcc +lat_1=44 +lat_2=49 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs  
## LDID: 0 
## Number of fields: 18 
##          name type length typeName
## 1   ID_GEOFLA    0     10  Integer
## 2   CODE_COMM    4      3   String
## 3   INSEE_COM    4      5   String
## 4    NOM_COMM    4     50   String
## 5      STATUT    4     20   String
## 6  X_CHF_LIEU    0      6  Integer
## 7  Y_CHF_LIEU    0      6  Integer
## 8  X_CENTROID    0      6  Integer
## 9  Y_CENTROID    0      6  Integer
## 10    Z_MOYEN    0      4  Integer
## 11 SUPERFICIE    0     10  Integer
## 12 POPULATION    2      6     Real
## 13  CODE_CANT    4      2   String
## 14   CODE_ARR    4      1   String
## 15  CODE_DEPT    4      2   String
## 16   NOM_DEPT    4     30   String
## 17   CODE_REG    4      2   String
## 18 NOM_REGION    4     30   String

Ici, on voit que les données sont correctement reconnues (driver adapté au format ESRI Shapefile, le SRS est reconnu, les données associées sont bien présentes). A noter que pour accéder directement aux données de coordonnées/projection (pour vérifier ou modifier ces informations), il est plus simple d’utiliser directement la fonction proj4string().

Les données étant OK, on peut importer le fichier SHP et ses données associées via la fonction readOGR(), afin de créer un objet nommé ici comm. On regarde ensuite la structure de l’objet crée par rgdal ;

# Import via la fonction readOGR de rgdal
# Pour info, d'autres outils existent (ex : fonction readShapeSpatial() du package maptools)
comm <- readOGR(dsn = pathToShp, layer="COMMUNE", stringsAsFactors=FALSE)
## OGR data source with driver: ESRI Shapefile 
## Source: "./data/shp", layer: "COMMUNE"
## with 36594 features and 18 fields
## Feature type: wkbPolygon with 2 dimensions
# Description de la structure globale
str(comm, 2)
## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  36594 obs. of  18 variables:
##   ..@ polygons   :List of 36594
##   .. .. [list output truncated]
##   ..@ plotOrder  : int [1:36594] 6722 3141 2093 4934 3161 1246 35885 14600 13687 5719 ...
##   ..@ bbox       : num [1:2, 1:2] 99226 6049647 1242375 7110480
##   .. ..- attr(*, "dimnames")=List of 2
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots

Les objets créés via readOGR ont des structures particulières. Ils sont divisés en plusieurs sous-objets, ou slots, préfixé par le symbole @. Chaque sous-objet possède sa propre structure, et ses propres informations. Pour y accéder, la commande suit le même motif que pour accéder à une variable dans un dataframe, mais via le symbole @ au lieu du symbole $.

Les informations de géométries, relatives au fichier .SHP, sont stockées dans les slots @polygons, @plotOrder et @bbox. Les données associées, issues du fichier .DBF, sont dans le slot @data. Enfin, les données de projection, issues du fichier .PROJ, sont dans un slot à part, @proj4string. Pour plus d’informations, vous pouvez regarder la documentation accompagnant cette classe particulière avec la commande class?SpatialPolygons.

Jetons par curiosité un oeil au slot @data :

# Description de la structure des données associées 
str(comm@data, 2)
## 'data.frame':    36594 obs. of  18 variables:
##  $ ID_GEOFLA : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ CODE_COMM : chr  "013" "152" "294" "044" ...
##  $ INSEE_COM : chr  "65013" "30152" "61294" "47044" ...
##  $ NOM_COMM  : chr  "ANSOST" "LES MAGES" "MORTREE" "CAHUZAC" ...
##  $ STATUT    : chr  "Commune simple" "Commune simple" "Chef-lieu canton" "Commune simple" ...
##  $ X_CHF_LIEU: int  4660 7931 4847 5065 5098 5004 7706 4823 9856 10482 ...
##  $ Y_CHF_LIEU: int  62635 63483 68416 63986 64012 64024 63949 62670 68140 68312 ...
##  $ X_CENTROID: int  4659 7938 4849 5067 5097 5012 7686 4815 9845 10493 ...
##  $ Y_CENTROID: int  62635 63483 68415 63986 64009 64016 63934 62673 68140 68307 ...
##  $ Z_MOYEN   : int  197 222 183 89 86 100 1139 231 366 145 ...
##  $ SUPERFICIE: int  221 1265 2373 807 1189 1357 4597 1100 2730 1164 ...
##  $ POPULATION: num  0.1 1.9 1.1 0.3 0.2 0.3 0.2 0.3 2.5 4.8 ...
##  $ CODE_CANT : chr  "18" "24" "24" "09" ...
##  $ CODE_ARR  : chr  "3" "1" "2" "3" ...
##  $ CODE_DEPT : chr  "65" "30" "61" "47" ...
##  $ NOM_DEPT  : chr  "HAUTES-PYRENEES" "GARD" "ORNE" "LOT-ET-GARONNE" ...
##  $ CODE_REG  : chr  "73" "91" "25" "72" ...
##  $ NOM_REGION: chr  "MIDI-PYRENEES" "LANGUEDOC-ROUSSILLON" "BASSE-NORMANDIE" "AQUITAINE" ...

On voit ici le nom des variables importées, et leur format (numérique, charactère…).

De même pour le slot @polygons :

# Description de la structure des données vectorielles 
head(comm@polygons, n=1)
## [[1]]
## An object of class "Polygons"
## Slot "Polygons":
## [[1]]
## An object of class "Polygon"
## Slot "labpt":
## [1]  465929.9 6263547.0
## 
## Slot "area":
## [1] 2124547
## 
## Slot "hole":
## [1] FALSE
## 
## Slot "ringDir":
## [1] 1
## 
## Slot "coords":
##        [,1]    [,2]
## [1,] 465804 6262421
## [2,] 465527 6262754
## [3,] 465237 6264352
## [4,] 465913 6264609
## [5,] 466532 6264260
## [6,] 466394 6262456
## [7,] 465804 6262421
## 
## 
## 
## Slot "plotOrder":
## [1] 1
## 
## Slot "labpt":
## [1]  465929.9 6263547.0
## 
## Slot "ID":
## [1] "0"
## 
## Slot "area":
## [1] 2124547

Nos données associées, comme le code INSEE de la commune, son nom, son statut administratif, sa superficie, sa population, etc… sont correctement importées dans notre slot @data.

On peut dès à présent fabriquer une carte très simple représentant les contours des communes, à partir des données spatiales extriates. A titre d’exemple, on va utiliser ici la nouvelle méthode de plot, spécifique à la classe SpatialPolygons et ajouté par le package sp. Concrètement, cet ajout permet d’utiliser la fonction classique plot pour représenter graphiquement des objets spatiaux.

# On représente les contours des communes de France avec la fonction plot
plot(comm)

Il y vraiment beaucoup de communes en France ! Du coup, une telle carte est difficilement lisible et particulièrement longue à créer… Et si on se limitait à l’Ile-de-France pour nos tests ?

Pour représenter uniquement les communes d’Ile-de-France (astuce : c’est la région numéro 11), il suffit d’appliquer une sélection “classique” à l’objet comm grâce aux variables à notre disposition, comme on pourrait le faire avec un dataframe :

# On plot les communes contenues dans l'Ile-de-France (CODE_REG égal à 11)
plot(comm[comm$CODE_REG=="11",],)

On obtient donc bien la carte des contours des communes d’Ile-de-France.

Supposons maintenant qu’on veut faire ressortir de cette carte les communes dont la population est de moins de 10 000 habitants. On peut chaîner les critères de sélections, comme on le fait habituellement via & ou | :

# Communes d'IDF avec moins de 10 000 habitants
plot(comm[comm$CODE_REG=="11" & comm$POPULATION<1,],)

Notre carte ressemble à du gruyère ! C’est normal, puisque R n’a représenté que les communes remplissant nos critères. Les autres communes n’ont donc pas été traitées et sont absentes de la carte…

Si on souhaite juste faire ressortir les communes d’intérêt, on peut commencer par représenter toutes les communes d’Ile-de-France, puis représenter en rouge les communes de moins de 10 000 habitants uniquement. Par exemple :

# Représentation de toutes les communes d'IDF
plot(comm[comm$CODE_REG=="11",])
# Ajout sur le graphique précédent (via add = TRUE) des communes d'intérêt, remplies en rouge (via col="red") 
plot(comm[comm$CODE_REG=="11" & comm$POPULATION<1,], col="red", add=TRUE)

C’est déjà mieux !

2.2 Utiliser des données externes dans notre carte

Nous avons vu comment créer une carte très simple à partir de données spatiales (issues du .SHP) et des données associées (issues du .DBF). Les données associées sont généralement très limitées, et ne contiennent pas le plus souvent les informations qui nous intéresse réellement. Il est donc nécessaire de marier nos données spatiales avec des données externes.

Pour cet exemple, j’ai choisi un jeu de données décrivant la distance parcourue par un patient pour se rendre sur son lieu d’hospitalisation. Cette distance est calculée ici en terme temps moyen de trajet (en minutes) pour accéder en voiture au lieu d’hospitalisation. Le jeu de donnée original est disponible en accès libre sur le site de l’OpenData. Comme d’habitude, le fichier est au format XLS et comporte des mentions (titre, précision, etc…) rendant son utilisation directe impossible. Pour gagner du temps, vous trouverez ici une version utilisable de la table de données.

Commençons par importer nos données :

# Emplacement des données
pathToData <- "./data/distance_hospit_Etalab.csv"
# Import des données en respectant le format "Character" des codes communes
tempsTrajet <- read.csv2(file = pathToData, stringsAsFactor = FALSE, colClasses = c("character", "integer")) 
# Verification import
str(tempsTrajet)
## 'data.frame':    36206 obs. of  2 variables:
##  $ commune: chr  "01001" "01002" "01004" "01005" ...
##  $ temps  : int  31 32 36 45 17 29 11 10 40 19 ...

Nous avons donc d’un coté notre objet comm comprenant nos données spatiales et les données associées, et d’un autre côté notre objet tempsTrajet comprenant nos données externes. L’idée est de rajouter nos données externes aux données associées déjà présentes dans notre slot @data. De cette façon, on conserve les informations de base (le nom des communes, leur numéro de département, etc…) tout en ajoutant les informations pertinentes pour nos analyses (les temps de trajet).

Plusieurs façons de réaliser cet ajout sont possibles. De façon générale, il faut faire attention à ne pas bousculer l’ordre des données du slot @data car ces données sont jointes ligne-à-ligne avec les données spatiale. On peut utiliser la fonction match(),la classique fonction merge(), voir la fonction left_join() issue du package dplyr (les habitués du SQL seront en terrain familier).

Pour la beauté de l’exemple, la fonction match() est utilisée ici. On constitue un nouveau dataframe, composé d’une part des données originales, et d’autre part des données de temps de trajet, rangées dans le même ordre que les données originales, grâce à match(). Les variables clefs permettant d’identifier uniquement chaque commune n’étant pas les mêmes dans les deux jeux de données, elles sont précisées manuellement.

# On remplace le slot @data par un nouveau tableau de données, contenant les anciennes données + les nouvelles
comm@data <- data.frame(comm@data, tempsTrajet[match(comm@data[, "INSEE_COM"],tempsTrajet[, "commune"]), ])
head(comm@data)
##   ID_GEOFLA CODE_COMM INSEE_COM                NOM_COMM           STATUT
## 0         1       013     65013                  ANSOST   Commune simple
## 1         2       152     30152               LES MAGES   Commune simple
## 2         3       294     61294                 MORTREE Chef-lieu canton
## 3         4       044     47044                 CAHUZAC   Commune simple
## 4         5       272     47272  SAINT-QUENTIN-DU-DROPT   Commune simple
## 5         6       373     24373 SAINT-AUBIN-DE-CADELECH   Commune simple
##   X_CHF_LIEU Y_CHF_LIEU X_CENTROID Y_CENTROID Z_MOYEN SUPERFICIE
## 0       4660      62635       4659      62635     197        221
## 1       7931      63483       7938      63483     222       1265
## 2       4847      68416       4849      68415     183       2373
## 3       5065      63986       5067      63986      89        807
## 4       5098      64012       5097      64009      86       1189
## 5       5004      64024       5012      64016     100       1357
##   POPULATION CODE_CANT CODE_ARR CODE_DEPT        NOM_DEPT CODE_REG
## 0        0.1        18        3        65 HAUTES-PYRENEES       73
## 1        1.9        24        1        30            GARD       91
## 2        1.1        24        2        61            ORNE       25
## 3        0.3        09        3        47  LOT-ET-GARONNE       72
## 4        0.2        09        3        47  LOT-ET-GARONNE       72
## 5        0.3        12        1        24        DORDOGNE       72
##             NOM_REGION commune temps
## 0        MIDI-PYRENEES   65013    30
## 1 LANGUEDOC-ROUSSILLON   30152    23
## 2      BASSE-NORMANDIE   61294    20
## 3            AQUITAINE   47044    39
## 4            AQUITAINE   47272    40
## 5            AQUITAINE   24373    26

On a bien fusionné nos infos, mais on remarque au passage que nous avons récupéré toutes les colonnes de tempsTrajet, y compris la variable commune qui fait doublon, mais passons.

En étant plus pragmatique, nous aurions pu faire un simple :

# Chargement de dplyr
library(dplyr)
# On remplace le slot @data par un nouveau tableau de données, contenant les anciennes données + les nouvelles
comm@data <- left_join(comm@data, tempsTrajet, by = c("INSEE_COM" = "commune"))

Imaginons maintenant qu’on souhaite produire la carte des temps de trajet, pour la région Ile-De-France. Nous allons réduire notre jeu de données pour ne garder que la région Ile-De-France. Nous allons aussi utiliser classInt pour discrétiser notre temps de trajet, et RColorBrewer pour choisir des couleurs adaptées à notre nouvelle variable discrète.

# Filtrage des données pour ne garder que la région Ile-De-France
idf <- comm[comm@data$CODE_REG=="11",]
# Vérification
str(idf, 2)
## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  1300 obs. of  20 variables:
##   ..@ polygons   :List of 1300
##   .. .. [list output truncated]
##   ..@ plotOrder  : int [1:1300] 131 80 563 163 100 31 1253 374 976 36 ...
##   ..@ bbox       : num [1:2, 1:2] 586509 6780072 740960 6904934
##   .. ..- attr(*, "dimnames")=List of 2
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
# Découpage du temps de trajet en 5 classes via la méthodes des quantiles : idenfication des bornes (breaks, ou brks)
classTemps <- classIntervals(idf@data$temps, 5, style = "quantile")
# Choix d'une palette de couleur pour les 5 catégories
palette <- brewer.pal(n = 5, name = "YlOrRd")


# Application de ce découpage à variable temps, sauvegarde dans temps_cat
# On stocke, pour chaque observation, la valeur de la couleur correspondante
idf@data$temps_cat <- as.character(cut(idf@data$temps, breaks = classTemps$brks, labels = palette, include.lowest = TRUE))

# On stocke l'information des classes pour créer une légende 
legende <- as.character(levels(cut(idf@data$temps, breaks = classTemps$brks, include.lowest = TRUE, right = FALSE)))

# Vérification
str(idf@data)
## 'data.frame':    1300 obs. of  21 variables:
##  $ ID_GEOFLA : int  101 211 328 345 388 555 562 684 686 687 ...
##  $ CODE_COMM : chr  "307" "175" "486" "001" ...
##  $ INSEE_COM : chr  "78307" "91175" "77486" "77001" ...
##  $ NOM_COMM  : chr  "HERMERAY" "CORBREUSE" "VAUDOY-EN-BRIE" "ACHERES-LA-FORET" ...
##  $ STATUT    : chr  "Commune simple" "Commune simple" "Commune simple" "Commune simple" ...
##  $ X_CHF_LIEU: int  6033 6232 7058 6677 5960 6550 6269 7015 7027 7246 ...
##  $ Y_CHF_LIEU: int  68396 68225 68431 68051 68633 68602 68258 68857 68580 68531 ...
##  $ X_CENTROID: int  6033 6231 7060 6677 5963 6574 6257 7018 7037 7245 ...
##  $ Y_CENTROID: int  68400 68221 68441 68054 68630 68596 68266 68855 68585 68543 ...
##  $ Z_MOYEN   : int  154 151 120 116 131 49 138 107 120 155 ...
##  $ SUPERFICIE: int  1835 1584 2750 1257 393 1636 3071 1909 1670 1037 ...
##  $ POPULATION: num  0.9 1.7 0.8 1.2 0.3 ...
##  $ CODE_CANT : chr  "18" "07" "27" "03" ...
##  $ CODE_ARR  : chr  "2" "1" "3" "4" ...
##  $ CODE_DEPT : chr  "78" "91" "77" "77" ...
##  $ NOM_DEPT  : chr  "YVELINES" "ESSONNE" "SEINE-ET-MARNE" "SEINE-ET-MARNE" ...
##  $ CODE_REG  : chr  "11" "11" "11" "11" ...
##  $ NOM_REGION: chr  "ILE-DE-FRANCE" "ILE-DE-FRANCE" "ILE-DE-FRANCE" "ILE-DE-FRANCE" ...
##  $ commune   : chr  "78307" "91175" "77486" "77001" ...
##  $ temps     : int  55 8 23 19 42 14 0 22 8 32 ...
##  $ temps_cat : chr  "#BD0026" "#FFFFB2" "#FD8D3C" "#FECC5C" ...

Il ne reste plus qu’à faire notre carte :

# Avec le package base
plot(idf, col = idf@data$temps_cat, border = "black")
legend("bottomleft", legend = legende, fill = palette, cex=0.6, title = "Temps de parcours moyen (minutes)")

3 Annexe

3.1 classInt

Si la fonction cut() (inclue par défaut dans R) permet de découper très facilement une variable continue en classes, elle nécessite une liste de valeurs seuils (ou breaks). La fonction classIntervals() du package classInt permet de calculer facilement ces valeurs seuils, en appliquant différentes méthodes. Par défaut, cette fonction retourne un objet contenant les valeurs seuils dans une variable nommée brks.

Pour plus d’information les différentes méthodes applicables pour le calcul des seuils, référez-vous à la documentation de ce package.

3.2 RBrewerColor

Le choix d’une palette de couleur pour un graphique n’est pas forcément facile. En fonction de la palette choisie, il est possible de faire ressortir (ou pas) différentes informations. Par exemple, l’utilisation d’une palette Kallant du jaune au rouge implique l’idée de progression entre les catégories : ce mode de représentation ne sera pas vraiment adapté à la représentation d’une variable qualitative… Dernier détail d’importance, certaines couleurs s’associent naturellement et vont sembler harmonieuses, d’autres non, ce qui peut fortement influencer la perception d’un graphique. Le package RBrewerColor facilite cette tâche, en proposant différentes palettes adaptées au nombre de classes (donc de couleurs) voulues, et adaptée au mode de représentation désiré (données séquentielles, divergeantes ou qualitatives). Les palettes séquentielles vont être adaptée à des données ordonnées et traduisent l’idée de progression (ex : jaune = faible valeur, rouge = forte valeur). Les palettes divergeantes permettent de mettre en avant les valeurs extrêmes, et en retrait les valeurs centrales. Enfin, les palettes qualitatives sont utilisables pour des variables où les catégories n’ont pas lien de valeur entre elles (ex : blond ou roux).

Pour avoir un aperçu des palettes disponibles, vous pouvez utiliser la fonction display.brewer.all(). Pour plus d’information (sur le nom des palettes, etc…), référez-vous à la documentation du package.