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

Laisser un commentaire

Votre adresse de messagerie ne sera pas publiée. Les champs obligatoires sont indiqués avec *