1 Comment réaliser une carte avec R ?

Voici quelques lignes de code pour vous accompagner à réaliser simplement une première carte à partir de données disponibles sur une plateforme de données.

J’ai choisi de ne traiter ici que des données vectorielles (pas d’images ou raster). Toutes les données sont disponibles sur le portail open data de Montpellier Mediterranée Métropole (MMM: https://data.montpellier3m.fr/).

1.1 Les librairies

Pour cette brève présentation, vous aurez besoin des librariries {sf} et {tmap}, sans oublier l’incontournable {tidyverse}.

Il faut distinguer la manipulation des données de la visualisation des données.

Le package {sf} (Spatial Feature) est dédié à la manipulation, la transformation et l’analyse de données spatiales. A la manière du Tidyverse, il combine les fonctionnalités de {sp}, {rgeos} et {rgdal} dans un package.

Le package {tmap} (thematic map) est quant à lui dédié à la visualisation des données spatiales. Cette fois-ci, à la manière du package {ggplot2}.

library(tidyverse)
library(sf)
library(tmap)

1.2 Les données

Les données brutes peuvent avoir de multiples formats:

  • shapefile
  • geojson
  • csv
  • DBI (PostGIS)
  • gpkg

1.2.1 Import des données

Nous utiliserons la fonction st_read() du package {sf}. Les données des exemples sont disponibles sur la plateforme open data de MMM. Vous pouvez bien sûr remplacer le chemin (l’url) par le chemin de données sur votre machine!

OSM_MMM_BNAC <- read_sf("https://data.montpellier3m.fr/sites/default/files/ressources/OSM_MMM_BNAC.geojson")

1.2.2 Class des données spatiales

class(OSM_MMM_BNAC)
## [1] "sf"         "tbl_df"     "tbl"        "data.frame"

Les données spatiales dans R sont de classe sf. Mais on constate aussi qu’elles sont de classes tbl_df, tbl et data.frame. Des classes d’objets que l’on connait bien et que nous manipulons au quotidien 🙂

1.2.3 Manipuler des données spatiales

Les fonctions du tidyverse permettent de manipuler les tables.

code_commune_select <- "34172" #Montpellier

OSM_MMM_BNAC %>% 
  filter(code_com_g==code_commune_select,ame_d=="VOIE VERTE")  
## Simple feature collection with 113 features and 30 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: 3.818371 ymin: 43.58251 xmax: 3.930179 ymax: 43.64231
## Geodetic CRS:  WGS 84
## # A tibble: 113 x 31
##    id_local  reseau_loc nom_loc id_osm  num_iti code_com_d ame_d regime_d sens_d
##  * <chr>     <chr>      <chr>   <chr>   <chr>   <chr>      <chr> <chr>    <chr> 
##  1 id_loc33~ ""         ""      way/33~ ""      34172      VOIE~ ""       BIDIR~
##  2 id_loc10~ ""         ""      way/10~ ""      34172      VOIE~ ""       BIDIR~
##  3 id_loc69~ ""         ""      way/69~ ""      34172      VOIE~ ""       BIDIR~
##  4 id_loc37~ ""         ""      way/37~ ""      34172      VOIE~ "AIRE P~ BIDIR~
##  5 id_loc55~ ""         ""      way/55~ ""      34172      VOIE~ "AIRE P~ BIDIR~
##  6 id_loc55~ ""         ""      way/55~ ""      34172      VOIE~ "AIRE P~ BIDIR~
##  7 id_loc24~ ""         ""      way/24~ ""      34172      VOIE~ "AIRE P~ BIDIR~
##  8 id_loc50~ ""         ""      way/50~ ""      34172      VOIE~ ""       BIDIR~
##  9 id_loc46~ ""         ""      way/46~ ""      34172      VOIE~ ""       BIDIR~
## 10 id_loc25~ ""         ""      way/25~ ""      34172      VOIE~ ""       BIDIR~
## # ... with 103 more rows, and 22 more variables: largeur_d <chr>,
## #   local_d <chr>, statut_d <chr>, revet_d <chr>, code_com_g <chr>,
## #   ame_g <chr>, regime_g <chr>, sens_g <chr>, largeur_g <chr>, local_g <chr>,
## #   statut_g <chr>, revet_g <chr>, access_ame <chr>, date_maj <date>,
## #   trafic_vit <chr>, lumiere <chr>, d_service <chr>, comm <chr>, source <chr>,
## #   project_c <chr>, ref_geo <chr>, geometry <LINESTRING [°]>

On constate qu’il y a systématiquement une colonne nommée geometry dans nos tables. Cette colonne geometry contient les coordonnées géographiques des individus (lignes). Pour faire simple, elles peuvent être :

  • Un point (POINT) ou plusieurs points (MULTIPOINT)
  • Une ligne (LINESTRING) ou plusieurs points (MULTILINESTRING)
  • Un polygone (POLYGON) ou plusieurs points (MULTIPOLYGON)

1.2.4 Regrouper ou fusionner les geométries

Il suffit d’utiliser la fonction summarise() du package dplyr.

OSM_MMM_BNAC %>% 
  filter(code_com_g==code_commune_select,ame_d=="VOIE VERTE") %>% 
  summarise()
## Simple feature collection with 1 feature and 0 fields
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: 3.818371 ymin: 43.58251 xmax: 3.930179 ymax: 43.64231
## Geodetic CRS:  WGS 84
## # A tibble: 1 x 1
##                                                                         geometry
##                                                            <MULTILINESTRING [°]>
## 1 ((3.842989 43.64178, 3.841227 43.64231, 3.841114 43.64231, 3.840252 43.64189,~

Si on souhaite traiter la table comme une table classique (regrouper sur des caractéristiques), il suffit d’utiliser la fonction as_tibble du package dplyr.

OSM_MMM_BNAC %>% 
  as_tibble() %>% 
  group_by(ame_d) %>% 
  summarise(nb=n())
## # A tibble: 8 x 2
##   ame_d                                            nb
##   <chr>                                         <int>
## 1 AMENAGEMENT MIXTE PIETON VELO HORS VOIE VERTE   941
## 2 AUCUN                                            52
## 3 AUTRE                                           177
## 4 BANDE CYCLABLE                                  315
## 5 COULOIR BUS+VELO                                119
## 6 PISTE CYCLABLE                                 1174
## 7 RAMPE                                             5
## 8 VOIE VERTE                                      241

On peut extraire des informations des géométries:

  • La longueur pour les lignes via la fonction st_length()du package {sf}
  • La surface pour les polygones via la fonction st_area()du package {sf}
OSM_MMM_BNAC %>% 
  mutate(longueur=st_length(geometry)) %>% 
  as_tibble() %>% 
  group_by(ame_d) %>% 
  summarise(longueur=sum(longueur)) 
## # A tibble: 8 x 2
##   ame_d                                         longueur
##   <chr>                                              [m]
## 1 AMENAGEMENT MIXTE PIETON VELO HORS VOIE VERTE  87604. 
## 2 AUCUN                                           8914. 
## 3 AUTRE                                          23746. 
## 4 BANDE CYCLABLE                                 37544. 
## 5 COULOIR BUS+VELO                               19027. 
## 6 PISTE CYCLABLE                                242654. 
## 7 RAMPE                                             24.9
## 8 VOIE VERTE                                     52839.

1.3 Visualiser les données

Au fil du temps, j’ai choisi le package {tmap} car il répond parfaitement à mon besoin de faire des cartes rapidement et plutôt jolies. Il a la même logique que ggplot2 (qui permet aussi de faire des cartes), couvre un large spectre de types de données spatiales et, c’est peut-être pour moi l’essentiel, me permet de réaliser avec le même code des cartes interactives ou statiques.

1.3.1 Une première carte

Pour cette première carte, nous allons afficher l’ensemble des aménagements cyclables de MMM. Il s’agit donc de lignes. On initialise toujours une carte dans un premier temps avec la fonction tm_shape() (un peut comme ggplot) puis on appelle la fonction appropriée aux géométries. Ici on appelle tm_lines().

carte1 <- tm_shape(OSM_MMM_BNAC)+
  tm_lines(col="#246D5F")

carte1

En fonction du type de géométries, on choisira (pour faire simple) la fonction adaptée:

  • tm_dots(), tm_bubbles(), tm_symbols(), tm_markers()les points
  • tm_lines() pour les lignes
  • tm_borders(), tm_fill(), tm_polygons() pour les polygones

1.3.2 Ajouter un titre

carte1 +
  tm_layout(title= 'Aménagements cyclables de Montpellier Méditerranée Métropole')

1.3.3 Ajouter une flêche nord

Les cartographes aiment qu’on leur rappelle le nord. Pas de problèmes avec la fonction tm_compass:

carte1 +
  tm_compass()

1.3.4 Ajouter la source

Afficher la source des données fait partie des bonnes pratiques aussi. Là encore, on assure le coup avec tm_credits() 🙂

carte1 +
  tm_credits("Source: Montpellier Mediterranée Métropole")

1.3.5 Superposer des couches

Vous pourrez superposer plusieurs cartes sans problème. Attention, il faudra vous assurer au préalable qu’elles aient la même projection avec st_crs().

# Contours des communes de Montpellier Méditerranée Métropole 
# https://data.montpellier3m.fr/dataset/contours-des-communes-de-montpellier-mediterranee-metropole
Contours_communes_MMM  <- read_sf("https://data.montpellier3m.fr/sites/default/files/ressources/MMM_MMM_Limites.geojson")

# Contours de Montpellier Méditerranée Métropole 
# https://data.montpellier3m.fr/dataset/contours-de-montpellier-mediterranee-metropole

Contours_MMM <- read_sf("https://data.montpellier3m.fr/sites/default/files/ressources/MMM_MMM_Contours.geojson")

st_crs(OSM_MMM_BNAC)==st_crs(Contours_communes_MMM)
## [1] TRUE

Si elles n’ont pas la même projection, vous pourres toujours utiliser st_transform().

carte2 <- tm_shape(Contours_communes_MMM)+
  tm_polygons(col = "grey",border.col = "white")+
  tm_shape(Contours_MMM)+
  tm_borders(col = "#666666",lwd=3)

carte2+
  carte1+
  tm_layout(main.title ='Aménagements cyclables de Montpellier Méditerranée Métropole',
            main.title.size=1,
            inner.margins=0.1)+
  tm_credits("Source: Montpellier Mediterranée Métropole",position = c(0,0))+
  tm_compass()

1.3.6 Carte thématique

En remplacant, le paramètre col qui jusque là contenait une valeur hexadécimale correspondant à un splendide vert, vous associerez alors une couleur à une modalité. Comme ici, où la variable ame_d de la table OSM_MMM_BNAC qualifie chaque type d’aménagement cyclable:

  • “AUTRE”
  • “AMENAGEMENT MIXTE PIETON VELO HORS VOIE VERTE”
  • “PISTE CYCLABLE”
  • “VOIE VERTE”
  • “BANDE CYCLABLE”
  • “AUCUN”
  • “COULOIR BUS+VELO”
  • “RAMPE”
carte3 <- tm_shape(OSM_MMM_BNAC)+
  tm_lines(col="ame_d") 

carte2+
  carte3+
  tm_layout(title ='Aménagements cyclables de Montpellier Méditerranée Métropole')+
  tm_credits("Source: Montpellier Mediterranée Métropole",position = c(0,0))+
  tm_compass()

La même carte avec une petite mise en forme…

carte3 <- tm_shape(OSM_MMM_BNAC %>% 
                     mutate(ame_d=str_trunc(ame_d,15,side="right",ellipsis="..."))
                     )+
  tm_lines(col="ame_d",title.col = "Type\nd'aménagement",
           palette="viridis") 

carte2+
  carte3+
  tm_layout(main.title ='Aménagements cyclables de Montpellier Méditerranée Métropole',
            main.title.size=1,
            inner.margins=0.1,
            legend.outside = T,
            legend.text.size=0.4)+
  tm_credits("Source: Montpellier Mediterranée Métropole",position = c(0,0))+
  tm_compass()

1.3.7 Les facettes

Comme sur ggplot2, on a la possibilité de réaliser des facettes. C’est parfois utile:

carte2+
  carte1+
  tm_layout(title= 'Aménagements cyclables de Montpellier Méditerranée Métropole',)+
  tm_credits("Source: Montpellier Mediterranée Métropole",position = c(0,0))+
  tm_compass()+
  tm_facets("ame_d",free.coords = F)

1.3.8 Interactivité

leaflet est une librairie javascript qui permet de réaliser des cartes interactives: https://leafletjs.com/

Cette librairie est à l’origine du package {leaflet} qui permet de créer des cartes interactives: https://rstudio.github.io/leaflet/

Vous n’aurez pas besoin de maîtriser ce package pour faire vos premières cartes interactives. Il vous suffira d’appliquer la fonction tmap_leaflet() à une carte tmap créée pour la rendre interactive (à la manière de ggplotly() pour un graphique ggplot).

tmap_leaflet(carte3)

1.4 Transformer un dataframe avec des variables latitude et longitude en sf

Il nous arrive parfois de devoir convertir un data frame en sf. C’est en particulier le cas quand on reçoit un fichier avec des coordonnées géographiques (longitude,latitude).

MMM_MMM_GeolocCompteurs <- read_csv("https://data.montpellier3m.fr/sites/default/files/ressources/MMM_MMM_GeolocCompteurs.csv")

MMM_MMM_GeolocCompteurs
## # A tibble: 46 x 6
##    `Nom du com`             `N° Série`  `N° Sér_1` Latitude Longitude OSM_Line_i
##    <chr>                    <chr>       <chr>         <dbl>     <dbl>      <dbl>
##  1 Compteur Vélo Tanneurs   XTH19101158 XTH191011~     43.6      3.87  188609530
##  2 Compteur Piéton/Vélo Be~ X2H19070220 X2H190702~     43.6      3.90  121403593
##  3 Compteur Vélo Lodève Ce~ Y2H20042633 X2H200426~     43.6      3.83  734202564
##  4 Compteur Vélo Lavérune   X2H20042632 X2H200426~     43.6      3.81   97705885
##  5 Compteur Vélo Vieille p~ Y2H20063161 X2H200631~     43.6      3.91  676645909
##  6 Compteur Vélo Delmas 2   Y2H20063164 X2H200631~     43.6      3.90  105575465
##  7 Compteur Vélo Delmas 1   Y2H20063163 X2H200631~     43.6      3.90  105575465
##  8 Compteur Vélo Gerhardt   Y2H20063162 X2H200631~     43.6      3.87   23231541
##  9 Compteur Vélo Lattes 2   Y2H20042634 X2H200426~     43.6      3.93   25871951
## 10 Compteur Vélo Lattes 1   Y2H20042635 X2H200426~     43.6      3.93  137058167
## # ... with 36 more rows

Il vous suffira d’utiliser la fonction st_as_sf() pour transformer votre table en objet spatial (sf) dans laquelle vous preciser les colonnes contenant les longitudes et latitudes.

MMM_MMM_GeolocCompteurs_sf <- MMM_MMM_GeolocCompteurs %>% 
  st_as_sf(coords = c("Longitude","Latitude"))

Une fois transformer, on peut bien sûr créer une nouvelle carte:

carte_compteur <- tm_shape(MMM_MMM_GeolocCompteurs_sf)+
  tm_dots(col="red",size = 0.1)
  

carte_compteur_finale <- carte2+
  carte3+
  carte_compteur+
  tm_layout(main.title ='Aménagements cyclables de Montpellier Méditerranée Métropole',
            main.title.size=1,
            inner.margins=0.1,
            legend.outside = T,
            legend.text.size=0.4)+
  tm_credits("Source: Montpellier Mediterranée Métropole",position = c(0,0))+
  tm_compass()+
  tm_add_legend('symbol', 
                col = "red",
                size = 0.1,
                labels = "Compteurs")

Et même apporter un peu d’interactions!

tmap_leaflet(carte_compteur_finale)

Alors, on se met en selle 🚴 ?