Archives de
Étiquette : Velov

Une carte des stations Vélo’V avec R (Partie 3 : Participer à OpenStreetMap)

Une carte des stations Vélo’V avec R (Partie 3 : Participer à OpenStreetMap)

OpenStreetMap est un formidable projet de cartographie collaborative en ligne. Dans deux précédents billets (ici et ) nous avons vu comment importer des données en provenance d’OSM dans R et un exemple (parmi une infinité) de traitement et de rendu cartographique. Mais les données que nous avons utilisées n’étaient pas vraiment complètes. Là où réside tout l’intérêt d’OSM c’est que tout le monde peut intervenir en complétant la base de donnée sur la base de ses propres connaissances, d’une sortie terrain ou encore avec l’imagerie satellite fournie par Bing.

Dans notre cas toutes les stations Vélo’V de la ville de Lyon ont déjà été localisé par les cartographes d’OSM mais pour un certain nombre d’entre elles les tags associés ne sont pas assez informatifs. On s’est par exemple intéressé à la capacité des stations (i.e. le nombre de vélos que l’on peut y garer) qui est renseignée dans la base de données avec le tag capacity. A l’heure où j’écris ces lignes – et cette précision est importante car les choses vont très vite sur OSM – la capacité est indisponible pour 87 stations. Il pourrait alors être intéressant de produire une carte de ces stations « incomplètes » pour optimiser le travail de terrain. Une carte statique comme on en a produit jusqu’à présent n’est pas idéale. L’idée ici est d’utiliser l’outil Google Maps qui permet d’une part d’accéder à une carte interactive mais aussi d’ajouter des points d’intérêts définis par l’utilisateur et accessoirement de profiter d’un module GPS.

Google Maps peut être très pratique pour s'orienter sur le terrain et localiser des éléments à cartographier

Dans ce billet, on réalisait l’importation des données à partir d’une extraction de la base téléchargée sur geofabrik mais dans le cas présent cette méthode se révèle lourde car il faudra re-télécharger le fichier à chaque fois qu’on voudra mettre à jour la carte. Ici il est plus intéressant de faire appel à une API pour extraire directement les données qui nous intéressent. J’utilise ici l’Overpass API, j’en ai essayé d’autre mais elles sont soit très peu réactives, soit HS, soit non à jour… Pour ce qui est du code on charge d’abord les données avec RCurl puis on parse le XML. Si les choses étaient simples on pourrait directement interpréter le XML avec la fonction as_osmar mais les API retournent souvent les objets en leur enlevant certains attributs ce qui pose problème à la fonction as_osmar. Une solution simple ici est de créer artificiellement les attributs manquants. Une fois l’objet converti au format sp, on exporte les données au format KML.

library(osmar)
library(sp)
library (XML)
library(rgdal)
osm.xml <- getURL("http://overpass-api.de/api/interpreter?data=node[amenity=bicycle_rental](45.64,4.57,45.86,5.12);out;", .encoding="UTF-8")
osm.xml <- xmlParse(osm.xml)
xmlSApply(xmlRoot(osm.xml), addAttributes,
"version"=1,
"timestamp"="2000-01-01T01:01:01Z",
"uid"=1,
"user"="a",
"changeset"=1)

stations <- as_osmar(osm.xml)
stations <-addtags(stations)
stations <- as_sp(stations, what="points")
stations$capacity <- as.numeric(as.character(stations@data$capacity))
stations.fix <- stations[which(stations$capacity==0 | is.na(stations$capacity)),]
writeOGR(stations.fix[,-5], dsn="stafix.kml", layer= "stations", driver="KML")

Il s’agit maintenant d’ouvrir le fichier KML avec Google Maps. Pour cela, il faut héberger son fichier KML soit sur votre propre serveur, soit à l’aide d’un service comme Google Sites. Ne reste plus qu’à ouvrir l’URL dans Google Maps et récupérer un lien abrégé pour retrouver sa carte sur n’importe quel appareil capable d’afficher Google Maps. Bien sur les données ne sont pas mises à jour de façon dynamique mais cette méthode a l’avantage d’être relativement simple et rapide à mettre en place.

Facebooktwitter
Une carte des stations Vélo’V avec R (Partie 2 : Density plot)

Une carte des stations Vélo’V avec R (Partie 2 : Density plot)

Je continue mon exploration des possibilités offertes par R pour la cartographie et l’analyse spatiale. Le nombre de packages et de fonctionnalités est assez impressionnant et on peut vite s’y perdre… A la suite de ce billet où je montre comment cartographier les stations Vélo’V de Lyon, j’ai voulu poursuivre en essayant générer un raster de la densité bidimensionnelle des stations dans la ville.

Matt Neale on Flickr (CC-BY)

Profiter de ggmap dans le système graphique traditionnel

La première étape était de pouvoir afficher un fond de carte dans le système graphique « classique » de R, le package ggmap ne permettant de travailler que dans le système « ggplot2 ». J’ai donc créé la fonction plotMAP qui permet d’afficher un objet ggmap créé avec get_map() dans le système traditionnel.

plotMAP <- function (x, bbox=NULL, crs="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0") {
  require(raster)
  lyr.bb <- as.numeric(attr(x, "bb"))
  lyr.m <- as.matrix(x)
  lyr.m.rgb <- apply(lyr.m, 2, col2rgb)
  lyr.a.rgb <- array(dim=c(dim(lyr.m),3))
    lyr.a.rgb[,,1] <- lyr.m.rgb[seq(from=1, to=dim(lyr.m.rgb)[1]-2, by=3), ]
    lyr.a.rgb[,,2] <- lyr.m.rgb[seq(from=2, to=dim(lyr.m.rgb)[1]-1, by=3), ]
    lyr.a.rgb[,,3] <- lyr.m.rgb[seq(from=3, to=dim(lyr.m.rgb)[1], by=3), ]
  lyr.brick <- brick(lyr.a.rgb, xmn=lyr.bb[2], xmx=lyr.bb[4], ymn=lyr.bb[1], ymx=lyr.bb[3], crs=crs)
  xsize <- lyr.brick@extent@xmax - lyr.brick@extent@xmin
  ysize <- lyr.brick@extent@ymax - lyr.brick@extent@ymin
  ratio <- max(xsize, ysize) / min(xsize, ysize)
  if (is.null(bbox)){
    plotRGB(lyr.brick, asp=ratio)
  } else {
    ext <- extent(bbox)
    plotRGB(lyr.brick, asp=ratio, ext=ext)
  }
}

La fonction s’appuie sur le package raster pour reconstituer l’image à partir d’un objet ggmap (x). On peut préciser une fenêtre avec un vecteur xmin, xmax, ymin, ymax et un système de projection (celui par défaut devrait convenir).

La carte de la densité spatiale des Vélo’V

Pour estimer ma fonction de densité spatiale j’ai décidé d’utiliser density.ppp() du package spatstat mais il existe beaucoup d’autres packages qui permettent de faire ça : MASS (kde2d), splancs (kernel2d)… Attention les yeux voilà le résultat !

La carte finale avec le fond Stamen, la densité des stations (heatmap et isolignes) et les stations (croix)

Je suis parti de l’objet stations (voir ce billet) et j’ai commencé par retirer les stations dont la capacité n’était pas renseignée. Ensuite petit passage par spatstat (conversion au format objet ppp obligatoire) pour estimer la fonction densité pondérée par la taille des stations et retour au format sp. Il ne reste plus qu’à afficher les différentes couches…

library(spatstat)
library(maptools)
library(ggmap)
library(sp)

stations$capacity <- as.numeric(as.character(stations@data$capacity))
stations.ok <- stations[which(stations$capacity>0 & !is.na(stations$capacity)),]

win <- owin(xrange=bbox(stations.ok)[c(1,3)], yrange=bbox(stations.ok)[c(2,4)])
stations.ppp <- ppp(coordinates(stations.ok)[,1], coordinates(stations.ok)[,2], window=win)
stations.ppp.den <- density(stations.ppp, sigma=0.005, weights=stations.ok$capacity)
stations.den <- as(stations.ppp.den, "SpatialGridDataFrame")

lyonraster <- get_map(bbox(stations.den)[c(1:4)], zoom=13, source="stamen", maptype="toner")
plotMAP(lyonraster)
image(stations.den, col = heat.colors(12, alpha=0.8), add=T)
contour(stations.den, col=2, lwd=2, add=T)
plot(stations.ok, add=T)

On distingue bien sur la carte deux zones de forte densité qui sont sans surprise la Presqu’ile et la Part Dieu. On observe aussi ce qui semble être un « effet campus » autour de la Doua et à proximité de l’ENS. On retrouve peut-être aussi cet effet vers la manufacture des tabacs (Lyon III). Ah mais ils ont pensé aux étudiants écolos et fauchés !  En revanche, on constate un « trou » vers Maisons Neuves. Une rapide inspection de la fonction G montre une tendance claire à la clusterisation. Pour trouver un Velo’V mieux vaut donc se trouver vers les Terraux. Mais c’est sans compter sur le subtil flux des utilisateurs. Ceux qui ont déjà erré vers l’Hôtel de Ville un samedi tard dans la nuit en quête du bicycle rouge salvateur ne me contrediront pas sur ce point…

Facebooktwitter
Une carte des stations Vélo’V avec R

Une carte des stations Vélo’V avec R

Je suis en train de lire (doucement) le livre de Bivand, Pebesma et Gómez-Rubio : Applied Spatial Data Analysis with R. Le livre est chouette, je le commenterai plus en détail dès que je l’aurai terminé. En attendant je me suis dit qu’il fallait que je pratique un peu le maniement de données spatiales. Et je suis parti je ne sais plus trop comment ni pourquoi sur les stations Vélo’v de Lyon. Ça me permet de jouer avec des données OpenStreetMap (mon péché mignon…). Ce billet explique comment extraire des données d’osm, les lire dans R et produire une carte.

La première étape consiste à récupérer les données d’osm. Je ne me suis pas encore penché sur l’extraction de données directement à partir de la base (est-ce seulement possible ?). Je pars donc d’une extraction de la région Rhône-Alpes qu’on peut trouver sur le site geofabrik au format XML .osm.

Une fois téléchargé, il n’est pas question d’ouvrir le fichier directement avec JOSM ou QGIS, il est bien trop lourd pour ça. On peut utiliser le logiciel Osmosis qui permet de faire les requêtes nécessaires en lignes de commandes. Ici en l’occurence je veux réduire ma bounding box :

osmosis --read-xml file="rhone-alpes.osm.bz2" --bounding-box top=45.8297562 left=4.7131348 bottom=45.6466883 right=5.0083923 --write-xml file="lyon_agglo.osm

et extraire les stations Vélo’v. On peux se contenter de sélectionner la valeur bicycle_rental sur la clé amenity car il n’y a que les Vélo’v qui propose ce service à Lyon.

osmosis --read-xml file="lyon_agglo.osm" --node-key-value keyValueList="amenity.bicycle_rental" --write-xml file="velov_stations.osm"

On a donc un fichier .osm qui contient les stations de Lyon. On peut le charger dans R avec le package osmar et le convertir en objet de la classe sp si on veut en faire une analyse plus poussée. Le truc qui me pose problème avec la fonction as_sp c’est qu’on perds toutes les métadonnées d’osm. Je fais donc tourner une petite fonction avant pour pouvoir retrouver ces informations dans le slot data de mon objet sp.

library(osmar)
library(sp)

stations <- get_osm(complete_file(), source = osmsource_file("velov_stations.osm"))

#Fonction qui ajoute les tags OSM au slot attrs d'un objet osmar pour conversion vers sp
 addtags <- function(obj){
 tag <- obj$nodes$tags
 tagf <- reshape(tag, idvar="id", timevar="k", direction="wide")
 colnames(tagf)[-1] <- substr(colnames(tagf)[-1], 3, 100000)
 obj$nodes$attrs <- merge(obj$nodes$attrs, tagf, by="id", all.x=TRUE)
 return(obj)
 }

stations <-addtags(stations)
stations <- as_sp(stations, what="points")

La classe sp a son lot de méthodes bien pratiques et permet de travailler avec beaucoup de packages. Un summary() nous donne un aperçu complet de l’objet et on peut afficher les positions relatives des points entre eux graphiquement avec la fonction plot().

Pour finir on peut afficher les points sur un fond de carte. J’ai utilisé le package ggmap qui fonctionne très bien et propose des tiles Google map, Osm, Stamen et CloudMade mais impose d’adopter la syntaxe ggplot2 et ne digère pas (encore?) les objets sp.

library(ggmap)
 lyonraster <- get_map("Lyon, France", zoom=13, source="osm")
 lyon <- ggmap(lyonraster, extent = "device")
 posist <- cbind(stations$lon,stations$lat, as.numeric(as.character(stations@data$capacity)))
 posist <- as.data.frame(posist)
 colnames(posist) <- c("lon","lat","capacity")
 lyon + geom_point(aes(x = lon, y = lat, size=capacity), data=posist, color=2)

Voilà ce que ça donne avec un rendu Stamen (toner). La taille des points est fonction de la taille de la station. Les stations de capacité nulle sont celles où l’info est manquante. Juste une question de temps 😉

Facebooktwitter