Données sources

On crée un data.frame qui contient les lieux que l’on souhaite géocoder, ie. quelques lieux de la géographie parisienne. On utilise ici le package tibble et sa fonction data_frame() plutôt que la fonction de base data.frame(). Les objets de type data_frame sont des data.frame améliorés, plus rapides, à la syntaxe moins verbeuse (plus besoin de faire appel aux sempiternels stringsAsFactors = FALSE etc.), et qui sont dotés d’un aperçu plus condensé et lisible. Ce format data_frame (ou tbl_df) a été apporté par le package dplyr et est extensivement utilisé dans l’ensemble des packages du Hadleyverse²

library(tibble)

geoplaces <- data_frame(Nom = c("Institut de Géographie",
                                   "Géographie-cités",
                                   "PRODIG",
                                   "Centre PMF",
                                   "Centre Montreal",
                                   "Olympe de Gouges"),  
                           
                           Adresse = c("Institut de Géographie, 75005 Paris, France",
                                      "Rue du Four, 75006 Paris, France",
                                      "Rue Valette, 75005 Paris, France",
                                      "90 rue de Tolbiac, 75013 Paris, France",
                                      "105 rue de Tolbiac, 75013 Paris, France",
                                      "Rue Albert Einstein, 75013 Paris, France")
                           )

Récuperer les coordonnées via une API de géocodage

On utilise pour cette étape le package photon³, créé par Timothée Giraud et décrit sur son blog. La fonction geocode() prend en entrée un vecteur des adresses et renvoie un data.frame contenant les adresses et coordonnées probables correspondant aux entrées. On choisit ici de limiter la sortie au meilleur résultat pour chaque adresse (paramètre limit), afin de pouvoir simplifier l’enrichissement des données initiales. Le paramètre lang définir la langue de l’adresse, et améliore ainsi la précision du géocodage.

N.B. : La fonction kable() (knit table), du package knitr permet d’afficher un data.frame, dans un document markdown tel que celui-ci, de plus élégante manière qu’un simple print().

library(photon)  # devtools::install_github(repo = 'rCarto/photon')

geoCodingResults <- geocode(geoplaces$Adresse, limit = 1, lang = "fr")

[1] "Institut de Géographie, 75005 Paris, France"
[1] "Rue du Four, 75006 Paris, France"
[1] "Rue Valette, 75005 Paris, France"
[1] "90 rue de Tolbiac, 75013 Paris, France"
[1] "105 rue de Tolbiac, 75013 Paris, France"
[1] "Rue Albert Einstein, 75013 Paris, France"

kable(geoCodingResults)

On enrichit le tableau initial

La fonction renvoie un nouveau data.frame contenant l’ensemble des informations liées au géocodage, ce qui permet d’évaluer la qualité de celui-ci et de vérifier que des erreurs n’ont pas été commises. On ne gardera ici que les colonnes correspondant aux coordonnées (exprimées en latitudes et longitudes, donc en degrés inscrits dans le système WGS84⁴). On va alors pouvoir enrichir le tableau initial en lui concaténant ces deux nouvelles colonnes résultantes. On isole les deux colonnes avec la fonction select() du package dplyr, et on les concatène au tableau initial avec bind_cols() issu du même package.

library(dplyr)

geogeoplaces <- geoplaces %>% bind_cols(select(geoCodingResults, lat, lon))

kable(geogeoplaces)

Conversion en format R

On a désormais un data.frame contenant les coordonnées, qu’il va falloir convertir dans un format spatial afin de pouvoir le cartographier. On réalise cette opération avec le package sp, qui va convertir le data.frame en objet SpatialPointsDataFrame, à partir des colonnes lon et lat, et lui associer le système de coordonnées géographique WGS 84⁵.

library(sp)
geogeogeoplaces <- as.data.frame(geogeoplaces, stringsAsFactors = FALSE)
coordinates(geogeogeoplaces) <- ~lon + lat
proj4string(geogeogeoplaces) <- CRS("+init=epsg:4326")

Cartographie statique

On peut alors afficher l’objet via la fonction plot(), ce qui donne une idée de la répartition des points, puis en faire une rapide carte statique avec le package ggmap et sa fonction éponyme, en allant chercher un fond de carte OpenStreetMap stylisé par Stamen avec get_map().

plot(geogeogeoplaces)

library(ggmap)
library(ggplot2)
baseMap <- get_map(location = geogeogeoplaces@bbox, source = "stamen", crop = FALSE, 
    maptype = "toner-lite")

ggmap(baseMap) + geom_point(data = geogeoplaces, aes(lon, lat), shape = 25, 
    fill = "blue", size = 3) + coord_map() + labs(x = "", y = "")

Cartographie dynamique

Afin de mener une exploration des données plus poussées, on peut faire appel à des outils de cartographie dynamiques, qui permettront par exemple de zoomer/dézoomer, ou encore de se déplacer pour la carte. Pour un usage “rapide” et interactif, on peut utiliser la fonction mapView() du package mapview, qui se révèle très pratique mais ne permet pas une personnalisation très poussée. Notons que cette fonction ne fonctionne pas dans le cadre de la création d’un document interactif tel que ce support, et que le résultat qu’elle produit n’est consultable que sur son propre ordinateur.

library(mapview)

mapView(geogeogeoplaces)

Afin de mieux personnaliser la sortie, on utilise ici le package leaflet, qui permet de produire une carte via l’outil du même nom.

library(leaflet)

leaflet(data = geogeogeoplaces) %>% 
addProviderTiles("CartoDB.Positron") %>% 
addMarkers(popup = ~as.character(Nom), icon = makeIcon("https://cdn2.iconfinder.com/data/icons/solar_system_png/512/Earth.png", 
    20, 20))

N.B. : Les syntaxes faisant appel aux pipes (tubes, %>% ici) sont de plus en plus fréquentes en R, en particulier avec l’utilisation de plus en plus courante de dplyr. Cet opérateur a été en premier lieu apporté par le package magrittr, et sa logique est de prendre en entrée l’argument qui le précède et de le fournir à la fonction qui suit. mean(scale(runif(n = 5))) est ainsi exprimé en runif(n = 5) %>% scale() %>% mean(), que l’on écrira le plus souvent sous la forme :

runif(n = 5) %>%
  scale() %>%
  mean()

Ce format rend la chaîne de traitement plus lisible, mieux documentable (on peut décrire sur chaque ligne ce qui est réalisé), et lui donne un mode de syntaxe fonctionnel.

Requête de récupération

Afin de faire appel à l’API d’OSRM, on utilisera le package osrm⁶, qui permet notamment de récupérer un objet directement interprétable par R plutôt qu’une longue de chaîne de caractères qu’il faudrait alors interpréter. La fonction osrmTrip() renvoie une liste des composantes connexes (segments non reliés) de l’itinéraire.

library(osrm)  # devtools::install_github('rCarto/osrm')

bestRoute <- osrmTrip(geogeogeoplaces)

Conversion en format R

Ce trajet ne contient qu’une composante connexe, on peut donc en récupérer le premier objet⁷, qui contient lui-même deux objets :

• trip, de type SpatialLinesDataFrame, contenant les segments géométriques de l’itinéraire.
• summary qui contient des informations générales sur l’itinéraire (durée et longueur totale).

geoRoute <- bestRoute[[1]]

summary(geoRoute$trip)

Object of class SpatialLinesDataFrame
Coordinates:
        min       max
x  2.332898  2.382103
y 48.826160 48.853168
Is projected: FALSE 
proj4string :
[+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
+towgs84=0,0,0]
Data attributes:
 start end      duration          distance     
 1:1   1:1   Min.   :0.02667   Min.   :0.0151  
 2:1   2:1   1st Qu.:2.60500   1st Qu.:1.5700  
 3:1   3:1   Median :2.68417   Median :1.6134  
 4:1   4:1   Mean   :3.19528   Mean   :2.1054  
 5:1   5:1   3rd Qu.:4.23083   3rd Qu.:2.7231  
 6:1   6:1   Max.   :6.45167   Max.   :4.7421  

geoRoute$summary

$duration
[1] 19.17167

$distance
[1] 12.6326

Cartographie

On peut alors cartographier l’itinéraire résultant, comme dans l’étape précédente.

plot(geoRoute$trip)

Et cartographier l’ensemble des éléments spatiaux récupérés jusque là, toujours avec les mêmes fonctions.

mapView(geoRoute$trip)

leaflet(data = geoRoute$trip) %>% 
addProviderTiles("CartoDB.Positron") %>% 
addPolylines(popup = ~sprintf("Distance : %s\nDurée :  %s", distance, duration)) %>% 
    
addMarkers(data = geogeogeoplaces, popup = ~as.character(Nom), icon = makeIcon("https://cdn2.iconfinder.com/data/icons/solar_system_png/512/Earth.png", 
    20, 20))

Segmenter le trajet

L’API renvoie l’élévation d’un point, alors que notre itinéraire est sous format linéaire. Il va donc falloir segmenter la polyligne afin d’obtenir les nœuds (points) qui la composent. Pour cela, on crée une fonction personnalisée qui va récupérer les coordonnées des nœuds composants les objets Lines, eux-mêmes contenus dans des lines⁸. On obtient en sortie un data.frame contenant les longitudes et latitudes de chacun des nœuds des polylignes.

library(magrittr)
coordsLists <- lapply(geoRoute$trip@lines, function(x) {
    x@Lines[[1]]@coords
})
coordsDF <- lapply(coordsLists, function(x) {
    as.data.frame(x)
}) %>% rbind_all() %>% set_colnames(c("lon", "lat"))

str(coordsDF)

Classes 'tbl_df', 'tbl' and 'data.frame':   50 obs. of  2 variables:
 $ lon: num  2.38 2.38 2.37 2.36 2.36 ...
 $ lat: num  48.8 48.8 48.8 48.8 48.8 ...

Échantilloner les points

La base de données du SRTM est principalement constituée d’un MNT dont la résolution est d’environ 90 mètres. Rien ne sert donc de chercher à récuperer l’altitude de points qui seraient plus proches que cette distance. En considérant qu’on dispose de environ 400 points pour un trajet de 12 km, soit un point tous les 30 mètres en moyenne, on va échantilloner tous les 4 points.

coordsDF$index <- as.numeric(row.names(coordsDF))
indexesToQuery <- seq(from = 1, to = nrow(coordsDF), by = 4)
altitudeDF <- coordsDF[seq(from = 1, to = nrow(coordsDF), by = 4), ]

Récuperer l’altitude

On peut alors lancer la récupération des élévations, grace à l’appel de la fonction GNsrtm3() du package geonames. Cette fonction prend en entrée la latitude et la longitude du point dont on souhaite l’altitude, il va donc falloir l’appliquer à chacun des points, par l’utilisation de la fonction apply. La fonction renvoie un data.frame constitué des coordonnées et de l’élévation (colonne srtm3, en mètres), que l’on va donc isoler et l’ajouter au data.frame utilisé pour l’appel.

library(geonames)
options(geonamesUsername = "parisgeo")  # L'utilisation d'un compte utilisateur est nécessaire pour utiliser cette API.

altitudePoints <- apply(altitudeDF, MARGIN = 1, FUN = function(x) {
    res <- GNsrtm3(lat = x["lat"], lng = x["lon"])
    res$srtm3
})

altitudeDF$alt <- altitudePoints

str(altitudeDF)

Classes 'tbl_df', 'tbl' and 'data.frame':   13 obs. of  4 variables:
 $ lon  : num  2.38 2.36 2.36 2.34 2.34 ...
 $ lat  : num  48.8 48.8 48.8 48.8 48.8 ...
 $ index: num  1 5 9 13 17 21 25 29 33 37 ...
 $ alt  : num  37 60 63 59 69 55 47 52 60 42 ...

Affichage du profil du trajet

On peut maintenant afficher le profil altimétrique via un graphique ggplot.

library(ggplot2)
ggplot(altitudeDF, aes(x = index * 30, y = alt)) + geom_line(group = 1) + coord_equal(ratio = 50) + 
    labs(x = "Distance", y = "Altitude")

L’ensemble du code de ce tutoriel est disponible sur le dépôt Git de cette séance du groupe ElementR : https://github.com/Groupe-ElementR/web-scraping