Archives de
Étiquette : R

Visualizing the R packages galaxy

Visualizing the R packages galaxy

The idea of this post is to create a kind of map of the R ecosystem showing dependencies relationships between packages. This is actually quite simple to quickly generate a graph object like that, since all the data we need can be obtained with one call to the function available.packages(). There are a few additional steps to clean the data and arrange them, but with the following code everyone can generate a graph of dependencies in less than 5s.

First, we load the igraph package and we retrieve the data from our current repository.

library(igraph)
dat <- available.packages()

For each package, we produce a character string with all its dependencies (Imports and Depends fields) separated with commas.

dat.imports <- paste(dat[, "Imports"], dat[, "Depends"], sep = ", ")
dat.imports <- as.list(dat.imports)
dat.imports <- lapply(dat.imports, function(x) gsub("\\(([^\\)]+)\\)", "", x))
dat.imports <- lapply(dat.imports, function(x) gsub("\n", "", x))
dat.imports <- lapply(dat.imports, function(x) gsub(" ", "", x))

Next step, we split the strings and we use the stack function to get the complete list of edges of the graph.

dat.imports <- sapply(dat.imports, function(x) strsplit(x, split = ","))
dat.imports <- lapply(dat.imports, function(x) x[!x %in% c("NA", "R")])
names(dat.imports) <- rownames(dat)
dat.imports <- stack(dat.imports)
dat.imports <- as.matrix(dat.imports)
dat.imports <- dat.imports[-which(dat.imports[, 1]==""), ]

Finally we create the graph with the list of edges. Here, I select the largest connected component because there are many isolated vertices which will make the graph harder to represent.

g <- graph.edgelist(dat.imports)
g <- decompose.graph(g)
g <- g[[which(sapply(g, vcount) == max(sapply(g, vcount)))]]

Now that we have the graph we can compute many graph-theory related statistics but here I just want to visualize these data. Plotting a large graph like that with igraph is possible but slow and tedious. I prefer to export the graph and open it in Gephi, a free and open-source software for network visualization.
The figure below was created with Gephi. Nodes are sized according to their number of incoming edges. This is pretty useless and uninformative, but still I like it. This looks like a sky map and it is great to feel that we, users and developers are part of this huge R galaxy.

R packages network
Click on the picture to see the HD version
Facebooktwitter
Introducing rleafmap. An R package for interactive maps with Leaflet.

Introducing rleafmap. An R package for interactive maps with Leaflet.

Obviously, I am late…

I released rleafmap about 1 year ago and I am just writing this blog post today. During this time, I presented the package to the french R-users community at the 3eme Rencontres R in Montpellier and could get some good feedbacks. Now, I would like to communicate better on the project. My idea is to post news about the development and communicate on new features illustrated with examples on this blog. The documentation and tutorials will be published on the project website (http://www.francoiskeck.fr/rleafmap/) if I can save time for that.

Purpose and philosophy

rleafmap is an R package that can be used to generate interactive maps with your data. If you manipulate spatial data in the R environment, at some point you probably want to visualize them. The most common way to visualize spatial data is maps. Like other packages (googleVis, rMaps…) rleafmap is designed to produce maps with interactivity to bring a richer experience to the end user. This is made possible by the use of Leaflet, the amazing open-source javascript library created by Vladimir Agafonkin.

There are two things important to be aware for a good start with rleafmap.

  • First, the package use exclusively input data inheriting from the Spatial class of the sp package. These classes are central and allows to work with many packages. If you don’t know it, a good place to start is the vignette of the sp package on CRAN. If you prefer a good book have look to Bivand et al. 2013 [1].
  • The second point is about how the package works. Data layers are stored in independent R object with their own symbology. The map consists in a compilation of these objects. The idea is to stick to the philosophy of the traditional GIS software.

For a more complete view of the package I strongly recommend that you have a look to the website.

[1] Bivand R.S., Pebesma E.J. & Gómez-Rubio V. (2013) Applied Spatial Data Analysis with R, 2nd edn. Springer, New York.

Facebooktwitter
Convertir un tableau R vers Javascript

Convertir un tableau R vers Javascript

Pour un petit projet sur lequel je travaille, j’ai eu besoin de passer des données d’un dataframe de R vers Javascript. Les données en Javascript peuvent être chargées au format JSON que R exporte avec le package RJSON. Mais pour ce projet particulier, dans un souci de simplicité et parce que je débute en Javascript, je voulais importer mes données dans le format traditionnel des arrays de Javascript.

Un array en JS ça marche comme ça :

var montableau = [1, 2, 3, 4];

Ce qui équivaut grosso modo à un vecteur sous R :

montableau <- c(1, 2, 3, 4)

On peut aussi faire de la pseudo-2D :

var montableau = [[1, 2], [3, 4]];

Ce qui donnerai une matrice sous R :

montableau <- matrix(c(1, 2, 3, 4), ncol=2, byrow=TRUE)

Bref voici une petite fonction R pour passer un dataframe (ou un truc ressemblant) de R vers Javascipt :

toJSarray <- function(df){

  if (!is.data.frame(df)){
    df <- as.data.frame(df)
  }
  
  df.nbli <- dim(df)[1]
  df.nbco <- dim(df)[2]
  temp <- vector()
  a.1 <- vector()
  
  for (i in 1:df.nbli){
      for (j in 1:df.nbco){
          if (is.numeric(df[,j])){
            a.1[j] <- df[i,j]
          } else {
            a.1[j] <- paste("\"", df[i,j], "\"", sep="")   
          }
      }
      a.2 <- paste(a.1, collapse=", ")
      temp[i] <- paste ("[", a.2, "]", sep="")
  }
  
  temp.2 <- paste(temp, collapse=", ")
  if (df.nbli == 1){
    jsarr <- temp.2
  } else {
    jsarr <- paste ("[", temp.2, "]", sep="")
  }
  invisible(jsarr)
}

Après quoi on peut rapidement écrire un fichier JS contenant nos données.

write(paste("var montableau = ", toJSarray(mondf), sep=""), "monfichier.js")

Facebooktwitter
Google au hasard

Google au hasard

Il m’arrive parfois, dans de grand moments de solitude, de me trouver béat à l’adresse www.google.fr, perdu et sans trop savoir ce que je fais là. Je m’étonne alors de voir comment mon esprit pense à la manière de trouver une réponse avant même que je n’ai formulé ma question. Mais c’est une évolution de mon comportement qui est somme toute assez logique puisque quelque soit la question, 9 fois sur 10 c’est bien par google que j’obtiens la réponse.

C’est comme ça que j’explique que je puisse échouer sur google sans avoir rien à lui demander. Et parfois face au curseur qui clignote j’ai ce réflexe (simiesque ?) de taper des touches au hasard.

Mykl Roventine on Flickr (CC BY-NC-SA)

Même si quand on y réfléchit on est surement assez loin du hasard, je suis toujours impressionné par la capacité du moteur à renvoyer des résultats. Par exemple la recherche sdbxce renvoie 612000 résultats (65 si on recherche le terme exact). On constate néanmoins et c’est d’ailleurs assez évident, que le nombre de résultat diminue à chaque nouvelle lettre ajoutée.

La petite expérience du jour consiste donc à sonder un peu la monstruosité de l’index de Google en déterminant jusqu’à combien de lettre frappées au hasard le bestiau nous renvoie des résultats.
On peut bien sur faire ça à la main mais R a quelques avantages décisifs :

  • Il est plus fort que nous pour simuler le hasard
  • Il peut faire les recherches Google à notre place ce qui n’est pas sans intérêt sachant qu’on en a plusieurs centaines au programme…

Commençons donc par générer nos termes de recherche. Une matrice 100×15 : 100 réplications de mots de 1 à 15 caractères tirés au hasard.

nb.letters <- 15
nb.samp <- 100

req <- data.frame()
for (i in 1:nb.letters){
  for(j in 1:nb.samp){
    samp <- sample(letters, i, replace=T)
    req[j, i] <- paste(samp, collapse="")
  }
}

Maintenant voici une fonction qui lance une recherche Google à partir d’une chaîne de caractères et retourne le nombre de résultats renvoyés par le moteur. J’ai choisi de rechercher les termes exacts (recherche entre guillemets) car je trouvais ça plus parlant.

library (XML)
library(RCurl)

goonum <- function(q){
  Sys.sleep(runif(1))
  url <- paste("http://www.google.fr/search?aq=f&gcx=w&sourceid=chrome&ie=UTF-8&q=\"", q, "\"", sep="")
  gpage <- getURL(url)
  gpage <- htmlParse(gpage)
  nbres <- xpathApply(gpage, "//div[@id='resultStats']", xmlValue)
  nbres.num <- unlist(gregexpr("[0-9]", nbres))
  nbres.val <- vector()
  for(i in 1:length(nbres.num)) 
    {nbres.val[i] <- substr(nbres, nbres.num[i], nbres.num[i])}
  nbres.val <- as.numeric(paste(nbres.val, collapse=""))
  nbres.val <- ifelse(is.na(nbres.val), 0, nbres.val)
  return(nbres.val)
}

Il ne reste qu’à faire tourner la fonction en boucle sur notre matrice de requêtes. Et faire deux petit graphiques qui sont assez parlant.

g.ans <- apply(req, FUN=goonum, MARGIN=c(1,2))
plot(apply(log(g.ans+1), 2,mean), type="b", ylab="Nombre de résultats (log+1)", xlab="Nombre de lettres")
plot(apply(g.ans==0, 2, sum), type="b", ylab="% de résultats nuls", xlab="Nombre de lettres")

Hum, ça pour être parlant… Quand on vient du monde de la biologie (et pire encore, de l’écologie) on hallucine toujours un peu face à ce genre de données ! C’est doux comme de la soie !

Comme on pouvait s’y attendre, le nombre de réponses décroit de manière exponentielle avec l’augmentation du nombre de lettre… jusqu’à 7 lettres. A 8 lettres 90% des requêtes ne renvoient plus aucun résultat. A partir de 9, plus rien, en tout cas dans cet échantillon. M’enfin 7 lettres c’est déjà pas mal, des combinaisons à 7 lettres il y en a 26^7 soit plus de 8 milliards…

Facebooktwitter
Quelques réflexions sur le permis de conduire

Quelques réflexions sur le permis de conduire

J’ai eu mon permis de conduire il y a quelques semaines et je n’ai même pas fêté ça. Si j’avais eu quelques potes sous la main, nul doute que nous nous serions emparés du prétexte mais doit on se réjouir d’avoir englouti une énergie considérable et des flots d’euros pour avoir le droit de devenir un pollueur actif et de jouir d’une illusion de liberté que l’on paye cher à l’industrie pétrolière ? Comme je n’était pas un élève très doué la petite blague m’aura en tout coûté plus de 2000 euros. Voici quelques enseignements que j’ai pu tirer de cette expérience dont je me serais volontiers passé. Je les partage avec vous, souhaitant bonne chance à ceux qui devront en passer par là…

Martin Gommel on Flickr (CC BY-NC-ND)

Ne comptez par sur votre banquier (surtout à la Banque Postale)

Tout le monde le sait, le permis ça coûte un bras. Certains s’en tirent pas trop mal et puis il y a les boulets de mon espèce qui ratent l’examen et à qui il faut quinze heures de leçon pour arriver à passer les rapports… Dans tout les cas arrive le moment douloureux où il faut sortir le chéquier. Comme la plupart des candidats sont des jeunes souvent désargentés, l’État a mis en place le « Permis à 1 euro ». Ça ne veut pas dire qu’on va vous faire le permis à un euro faut pas déconner quand même ! Ça veut juste dire qu’une banque vous fait un prêt que vous remboursez 1 euro/jour et l’État paye les intérêts à la banque. Le principe est pas mal, voyons la pratique.

Je suis à la Banque Postale, je me pointe la gueule enfarinée au guichet. La fille ne comprends rien, elle n’a jamais vu le document que je lui tends. Finalement un responsable m’explique que je dois prendre rendez vous avec mon conseiller financier. Deux jour plus tard je me pointe au rendez-vous, ça dure plus de deux heures pendant lesquels le type me ré-explique toute la procédure puis fait le tour de mes comptes, puis demande des documents sur mon garant. A la fin il me dit que naturellement ce n’est pas lui qui va prendre la décision de m’accorder le prêt. Il téléphone à Paris, et me passe une fille qui à nouveau me repose 50 questions… Je sors du bureau, dubitatif. C’était comme si je leur demandais un prêt de 250000 euros. Au bout d’une semaine la grosse blague, je reçois un courrier me notifiant qu’après simulation et examen de mon dossier la banque me refuse ma demande de prêt. Mon conseiller m’explique qu’il en est désolé mais qu’il ne peut rien pour moi. Tout ça se décide à Paris dans des conditions pour le moins opaque.

C’est juste un scandale, ces abrutis de banquiers me refuse un prêt de 1200 euros alors qu’il font joujou avec des millions sur les marchés financiers. Précisons que j’étais leur client depuis plus de cinq ans et que durant tout ce temps je n’ai pas été dans le rouge plus de cinq fois. Mon garant était fonctionnaire de catégorie A avec un revenu fixe et largement suffisant pour un tel prêt. Autrement dit un paquet de gens doivent se pointer avec un moins bon dossier et se faire gentiment remballer. Ironique quand on sait que ces mesures ont été prises pour aider ceux qui ont le moins de moyens…

Jouer à GTA ou à Need for Speed ne vous servira à rien

C’est bien malheureux mais les heures que j’ai pu passer au volant dans les rues de Liberty City et les tournois endiablés à bord de monstrueuses caisses couvertes de néons n’ont pas fait de moi un as de la route. En fait ces jeux ont même plutôt tendance à donner de mauvais réflexes car écraser une mamie est éliminatoire le jour de l’examen. Mais que cela ne vous dissuade pas de vous amuser ! J’ai remarqué qu’après une leçon épuisante de conduite, une petite course sur NFS atteint un niveau de fun record.

Misez sur vos facultés intellectuelles à défaut de savoir conduire

Le permis est avant tout une question de maîtrise technique que l’on acquiert au fur et à mesure des leçons. Comme partout, certain sont plus doués et apprennent plus vite. D’autres (comme moi) galèrent beaucoup ! A ceux là je conseille de ne pas se décourager et d’envisager la conduite accompagnée (voir point suivant) mais aussi de compenser avec leur neurones. Car dans le permis aussi il y a de la théorie. L’examen du code par exemple est uniquement basé sur la connaissance du code de la route et l’interprétation de situations. Il y a des gens qui mettent des années à avoir leur code en enchaînant les échecs à l’examen et les auto-école en profitent à fond avec plein de petites lignes glissées dans le contrat, qui à la fin, coûtent beaucoup d’argent. La solution que j’ai adoptée est simple : travail intensif ! Je me suis mangé leur bouquin immonde et tout leurs cours vidéo en 2 semaines, puis pendant une semaine j’ai enchainé les tests (attention de veiller à ce que l’auto-école donne accès à une banque de tests en ligne). Y a que ça qui marche, au bout de 400 questions on voit venir tout les pièges. J’ai donc eu mon code en 3 semaines au tarif minimum.

Pour l’épreuve pratique de conduite les choses se corsent mais là aussi un peu de travail à la maison permet de gratter quelques points. Connaitre les commandes, effectuer les vérifications et s’installer correctement dans le véhicule, ça peut rapporter jusqu’à 9 points (sur 30). Au passage je vous conseille ce site bien foutu pour réviser.

Tablez sur la conduite accompagnée et restez vous-même (i.e. un gros geek)

La conduite accompagnée permet de prendre confiance au volant et de terroriser son accompagnateur. Il serait donc dommage de s’en priver. En plus pour les majeurs, il existe une formule light : la conduite supervisée.

En ouvrant mon livret d’apprentissage je constate que je dois relever chaque trajet avec la date et la distance parcourue. Hahaha je vais te faire ça avec Excel ça va être rigolo en plus d’être joli. Mais je me rends vite compte que ça me soule d’aller chercher les distances sur google map. R peut nous faire tout ça tout seul et ça va être très rigolo en plus d’être très joli. Maintenant vous savez comment j’en suis arrivé à développer les fonctions présentées dans mon dernier billet. Pour les exploiter dans le cadre de la conduite accompagnée, il suffit de remplir un tableau comme ça :

DATE FROMCITY FROMDEP TOCITY TODEP AR
24/03/2012 LYON 69 GRENOBLE 38 NO

On précise le n° de département pour lever les ambiguïtés lors de la résolution d’adresse. La colonne AR est utilisée pour préciser si le trajet est un aller-retour.

source("script_mapquest.R")
tab <- read.table("conduite.txt", h=T)
tab$FROMDEP <- sprintf("%02.0f", tab$FROMDEP)
tab$TODEP <- sprintf("%02.0f", tab$TODEP)

from.cont <- paste(tab$FROMCITY, tab$FROMDEP, sep="+")
to.cont <- paste(tab$TOCITY, tab$TODEP, sep="+")

DIST <- mapq.distance(from.cont, to.cont, print=FALSE)
DIST <- DIST*ifelse(tab$AR=="NO", 1, 2)
sum(DIST)

DATE <- strptime(tab$DATE, format="%d/%m/%Y")
plot(DATE,cumsum(DIST), type="b", pch=20, ylab="Distance Cumulée (Km)", xlab="Date")
abline(h=c(1000), lty="dashed")

La commande source sert à charger les fonctions que l’on trouve ici. On utilise sprintf pour formater les chiffres car R à la fâcheuse habitude de supprimer les 0 qui précèdent les nombres. La ligne 11 renvoie le nombre total de km parcourus et on peut faire un petit graphique de son expérience, qui va croissante.

Distance cumulée en fonction du temps. Ojectif 1000 km atteint :p

En guise de conclusion : réformez moi ce système à la con

Voilà ma conclusion. Le permis de conduire en France a grand besoin d’être réformé en profondeur. Et ça beaucoup de gens le disent sans que ça bouge pour autant. Pourtant je suis d’avis qu’il ne faut pas y aller avec le dos de la cuillère. Je rêve d’une nationalisation des auto-écoles que l’on intégrerait aux lycées. Mon immersion au sein de la mafia des auto-écoles et mes mésaventures à la banque ne m’ont pas transformé en communiste… mais presque ! L’apprentissage de la conduite pourrait être une option facultative (mais fortement conseillée), pour laquelle des créneaux serait réservés dans l’emploi du temps. Le prix de la formation serait fixe, indexé sur le QF et l’évaluation serait continue. Ce ne sont que quelques idées lancées en l’air et ça peut sonner un peu extrême mais si on maintient que le permis de conduire est un facteur important pour l’accès à l’emploi, il serait temps de remettre un peu d’égalité dans ce bazar.

Facebooktwitter
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 &lt;- function(query, format=&quot;xml&quot;, limit=1, details=1)
  {
  urlroot &lt;- &quot;http://open.mapquestapi.com/nominatim/v1/search?&quot;
  url &lt;- paste(urlroot, &quot;format=&quot;, format, &quot;&amp;q=&quot;, query, &quot;&amp;addressdetails=&quot;, details, &quot;&amp;limit=&quot;, limit, sep=&quot;&quot;)

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


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

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

Package maptools sous Linux, attention aux majuscules

Package maptools sous Linux, attention aux majuscules

Si cet article peut vous économiser une crise de nerfs, autant l’écrire. Hier, j’ai voulu me faire une jolie carte des dernières élections en utilisant le shapefile des communes de France mis à disposition par l’IGN pour un chouette résultat inspiré de ça.

Tout allait très bien jusqu’à ce que je cherche à charger le shapefile dans R avec la fonction readShapeSpatial du package maptools.

commune <- readShapeSpatial("COMMUNE")

Dans les heures qui ont suivies, R m’a renvoyé éternellement la même erreur : Erreur dans read.dbf(filen1) : impossible d’ouvrir le fichier DBF.

Après des heures de prise de tête voilà la solution au problème : la fonction est sensible à la casse de l’extension du fichier *.dbf et pas à celles des *.shp et *.shx. L’IGN fournissant ses fichiers avec l’extension en majuscule, il suffisait donc de passer en minuscules le *.DBF

Grosse galère pour un détail, comme souvent hélas ! A noter que ce problème se pose sous Linux et non sous Windows.


Facebooktwitter