Robin Cura
ElementR, Vague 3, Séance 1 - 30/03/2016
Ce tutoriel vise à donner un exemple d'utilisation d'API^[Une API, ou Interface de Programmation est un service mis à disposition d'un public (complet ou restreint) lui permettant de récupérer du contenu (web ici) de manière automatique.] depuis R, autour d'un exercice complet qui consiste à préparer le parcours d'un événement scientifique inscrit dans un certain nombre de "lieux de la géographie" quantitative parisienne. La démarche est composée de 3 étapes :
-
Doter les lieux choisis de coordonnées géographiques, c'est le géocodage,
-
Rechercher un trajet entre ces lieux, c'est une opération de calcul d'itinéraire,
-
Établir le profil altimétrique du trajet résultant afin de contextualiser l'itinéraire.
Chacune de ces étapes fait appel à une API libre différente, choisie pour être aussi simple que possible à utiliser. Il existe naturellement de nombreuses autres manières de parvenir aux mêmes résultats, la démarche choisit ici est donc uniquement illustrative et pédagogique.
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^[Ensemble des packages créés par Hadley Wickham, un développeur R prolifique qui a notamment créé ggplot2, reshape2, (d)plyr. En voir une présentation complète]
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")
)On utilise pour cette étape le package photon^[En consulter la page de développement], 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)location name housenumber street postcode city state country osm_key osm_value lon lat msg
Institut de Géographie, 75005 Paris, France Institut de Géographie NA Rue Pierre et Marie Curie 75005 Paris Île-de-France France amenity university 2.342589 48.84469 NA
Rue du Four, 75006 Paris, France Rue du Four NA NA 75006 Paris Île-de-France France highway tertiary 2.333937 48.85270 NA
Rue Valette, 75005 Paris, France Rue Valette NA NA 75005 Paris Île-de-France France highway residential 2.346782 48.84749 NA
90 rue de Tolbiac, 75013 Paris, France NA 90 Rue de Tolbiac 75013 Paris Île-de-France France place house 2.364828 48.82695 NA
105 rue de Tolbiac, 75013 Paris, France NA 105 Rue de Tolbiac 75013 Paris Île-de-France France place house 2.364743 48.82658 NA
Rue Albert Einstein, 75013 Paris, France Rue Albert Einstein NA NA 75013 Paris Île-de-France France highway residential 2.381274 48.82741 NA
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^[Comme à peu près toutes les informations spatiales disponibles sur Internet.]).
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)Nom Adresse lat lon
Institut de Géographie Institut de Géographie, 75005 Paris, France 48.84469 2.342589 Géographie-cités Rue du Four, 75006 Paris, France 48.85270 2.333937 PRODIG Rue Valette, 75005 Paris, France 48.84749 2.346782 Centre PMF 90 rue de Tolbiac, 75013 Paris, France 48.82695 2.364828 Centre Montreal 105 rue de Tolbiac, 75013 Paris, France 48.82658 2.364743 Olympe de Gouges Rue Albert Einstein, 75013 Paris, France 48.82741 2.381274
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^[Dont le code EPSG est 4326].
library(sp)
geogeogeoplaces <- as.data.frame(geogeoplaces, stringsAsFactors = FALSE)
coordinates(geogeogeoplaces) <- ~lon + lat
proj4string(geogeogeoplaces) <- CRS("+init=epsg:4326")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 = "")
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))![https://cdn2.iconfinder.com/data/icons/solar_system_png/512/Earth.png","index":0},"iconRetinaUrl":{"data":20,"index":0},"iconWidth":20},null,null,{"clickable":true,"draggable":false,"keyboard":true,"title":"","alt":"","zIndexOffset":0,"opacity":1,"riseOnHover":false,"riseOffset":250},["Institut](https://cdn2.iconfinder.com/data/icons/solar_system_png/512/Earth.png","index":0},"iconRetinaUrl":{"data":20,"index":0},"iconWidth":20},null,null,{"clickable":true,"draggable":false,"keyboard":true,"title":"","alt":"","zIndexOffset":0,"opacity":1,"riseOnHover":false,"riseOffset":250},["Institut){kind=link}
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.
A partir de ce jeu de données spatialisées, nous allons maintenant chercher à trouver un itinéraire optimal entre tous ces points. Cette démarche fait appel au problème du voyageur de commerce, c'est-à-dire qu'on cherche l'itinéraire qui desserve chacun des points et revienne au point de départ tout en minimisant la distance parcourue. Il existe des outils permettant d'effectuer ce calcul depuis R, via des approximations successives concourant vers un optimum, mais on a ici choisi d'utiliser un service dédié, via l'API d'OSRM. Ce service de calcul d'itinéraire se base sur les données OpenStreetMap et permet d'en retirer des itinéraires précis susceptibles de répondre à différentes contraintes (itinéraires cyclistes, piétons etc.).
Afin de faire appel à l'API d'OSRM, on utilisera le package osrm^[Lui aussi développé par Timothée Giraud], 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)Ce trajet ne contient qu'une composante connexe, on peut donc en récupérer le premier objet^[Via l'utilisation des doubles crochets -- [[], qui contient lui-même deux objets :
trip, de typeSpatialLinesDataFrame, contenant les segments géométriques de l'itinéraire.summaryqui 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
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))Nous cherchons maintenant à enrichir cette information d'itinéraire, en prenant l'exemple de l'ajout d'un profil altimétrique, que l'on pourra afficher à côté de la carte.
Pour récupérer cette information sur l'altitude du trajet, on va utiliser l'API de Geonames, qui s'appuie notamment sur les données du SRTM afin de communiquer l'élévation d'un point donné. Cette API est accessible dans R via le package geonames.
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^[La manière dont sont structurés les objets spatiaux de sp n'est de toute évidence pas la plus évidente et simple à comprendre.].
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 ...
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), ]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 ...
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