Archives de
Mois : septembre 2012

Distances routières avec R et l’API Mapquest

Distances routières avec R et l’API Mapquest

Voilà un petit bout de code que j’avais écrit pour obtenir rapidement la distance par la route entre deux lieux géographiques avec R. On utilise les librairies RCurl et XML pour interroger l’API Mapquest basée elle même sur les données d’OSM. Le code est constitué de trois fonctions :

  • mapq.loc() renvoie un certain nombre d’infos sur un lieu (y compris ses coordonnées géographiques)
  • mapq.dist() renvoie la distance entre deux points géographiques spécifiés par leurs coordonnées
  • mapq.distance() : renvoie la distance routière entre deux points géographiques spécifiés par leurs noms. Choix de l’unité possible (« k » pour km ou « m » pour miles)

Une petite application de ces fonction est prévue pour le prochain billet 😉
On peut noter qu’une fonction équivalente est désormais fournie avec le package ggmap, basée elle sur le service Google Maps.

library (RCurl)
library (XML)

mapq.loc <- function(query, format="xml", limit=1, details=1)
  {
  urlroot <- "http://open.mapquestapi.com/nominatim/v1/search?"
  url <- paste(urlroot, "format=", format, "&q=", query, "&addressdetails=", details, "&limit=", limit, sep="")

  xmllocinfo <- getURL(url)
  xmllocinfo <- htmlParse(xmllocinfo)
  xmllocinfo <- xmlRoot(xmllocinfo)
   
  verif <- function(express) tryCatch(express, error=function(e) NULL)
  locinfo <- list()
  locinfo$lat <- verif(as.numeric(unlist(xpathApply(xmllocinfo, "/html/body/searchresults/place/@lat"))))
  locinfo$lon <- verif(as.numeric(unlist(xpathApply(xmllocinfo, "/html/body/searchresults/place/@lon"))))
  locinfo$city <- verif(as.character(xpathApply(xmllocinfo, "/html/body/searchresults/place/city", xmlValue)))
  locinfo$county <- verif(as.character(xpathApply(xmllocinfo, "/html/body/searchresults/place/county", xmlValue)))
  locinfo$state <- verif(as.character(xpathApply(xmllocinfo, "/html/body/searchresults/place/state", xmlValue)))
  locinfo$country <- verif(as.character(xpathApply(xmllocinfo, "/html/body/searchresults/place/country", xmlValue)))
  locinfo$countrycode <- verif(as.character(xpathApply(xmllocinfo, "/html/body/searchresults/place/country_code", xmlValue)))
  return(locinfo)
  }


mapq.dist <- function(Alat, Alon, Blat, Blon, format="xml", unit="k")
  {
  urlroot <- "http://open.mapquestapi.com/directions/v0/route?"
  url <- paste(urlroot, "outFormat=", format, "&unit=", unit, "&from=", Alat, ",", Alon, "&to=", Blat, ",", Blon, sep="")
  
  xmldistinfo <- getURL(url)
  xmldistinfo <- htmlParse(xmldistinfo)
  xmldistinfo <- xmlRoot(xmldistinfo)
  
  verif <- function(express) tryCatch(express, error=function(e) NULL)
  distance <- verif(as.numeric(xpathApply(xmldistinfo, "/html/body/response/route/distance", xmlValue)))
  return(distance)
  }

mapq.distance <- function(from, to, unit="k", print=TRUE)
  {
  from.f <- chartr(" ","+", from)
  to.f <- chartr(" ","+", to)
  from.info <- mapq.loc(from.f)
  to.info <- mapq.loc(to.f)
  distance <- mapq.dist(from.info$lat, from.info$lon, to.info$lat, to.info$lon, unit=unit)
  if(unit=="k") unit.print <- "Km"
  if(unit=="m") unit.print <- "Miles"
  if (print==TRUE){
    cat("Distance between", "\n")
    cat(from, ",", from.info$state, ",", toupper(from.info$countrycode), "and", to, ",", to.info$state, ",", toupper(to.info$countrycode), "\n")
    cat(distance, unit.print, "\n")
  } else {
    invisible(distance)
  }
  }
Facebooktwitter
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