Histoire de continuer la série de billets sur la visualisation, et la manipulation de données ouvertes, je vais reprendre de codes de Tony, de la formation Data Science pour l’Actuariat, pour visualiser le transport dans Paris (et la région parisienne). Si j’ai le temps, dans les jours a venir, je ferais une analyse du réseau de métro, compare aux autres grandes villes européennes. Pour commencer, on va récupérer les données, fournies par le site d’open data du stif, le syndicat des transports d’Ile de France (https://opendata.stif.info). Les données sont découpées par semestre, ce qui rend le code un peu lourd… mais bon, ça n’est pas plus complique pour autant.
library(dplyr) library(stringr) library(ggplot2) library(xlsx) library(ggmap) |
On commence par lire tous les fichiers en ligne
nbvalid = list() download.file("https://opendata.stif.info/explore/dataset/emplacement-des-gares-idf-data-generalisee/download/?format=csv&timezone=Europe/Berlin&use_labels_for_header=true","Gares.csv") gares = read.csv("Gares.csv", sep=";", header = TRUE) distr_pers = list() download.file("https://opendata.stif.info/explore/dataset/validations-sur-le-reseau-ferre-profils-horaires-par-jour-type-1er-sem/download/?format=csv&timezone=Europe/Berlin&use_labels_for_header=true","horaires1.csv") distr_pers$S1 = read.csv("horaires1.csv", sep=";", header = TRUE) download.file("https://opendata.stif.info/explore/dataset/validations-sur-le-reseau-ferre-profils-horaires-par-jour-type-2e-sem/download/?format=csv&timezone=Europe/Berlin&use_labels_for_header=true","horaires2.csv") distr_pers$S2 = read.csv("horaires2.csv", sep=";", header = TRUE) download.file("https://opendata.stif.info/explore/dataset/validations-sur-le-reseau-ferre-nombre-de-validations-par-jour-1er-sem/download/?format=csv&timezone=Europe/Berlin&use_labels_for_header=true","validations1.csv") nbvalid$S1 = read.csv("validations1.csv", sep=";", header = TRUE) download.file("https://opendata.stif.info/explore/dataset/validations-sur-le-reseau-ferre-nombre-de-validations-par-jour-2e-sem/download/?format=csv&timezone=Europe/Berlin&use_labels_for_header=true","validations2.csv") nbvalid$S2 = read.csv("validations2.csv", sep=";", header = TRUE) download.file("https://freakonometrics.free.fr/Correspondance_NOM.csv","Correspondance_NOM.csv") Cooresp = read.csv("Correspondance_NOM.csv", sep=";", header = TRUE) |
On a ensuite besoin de définir les dates des vacances, pour 2017
Vacances = list() Vacances$Noel = append(seq(from = as.Date("01/01/2017", format="%d/%m/%Y"), to=as.Date("02/01/2017", format="%d/%m/%Y"), by=1),seq(from = as.Date("24/12/2017", format="%d/%m/%Y"), to=as.Date("31/12/2017", format="%d/%m/%Y"), by=1)) Vacances$Ski = seq(from = as.Date("04/02/2017", format="%d/%m/%Y"), to=as.Date("19/02/2017", format="%d/%m/%Y"), by=1) Vacances$Printemps = seq(from = as.Date("02/04/2017", format="%d/%m/%Y"), to=as.Date("17/04/2017", format="%d/%m/%Y"), by=1) Vacances$Ete = seq(from = as.Date("08/07/2017", format="%d/%m/%Y"), to=as.Date("03/09/2017", format="%d/%m/%Y"), by=1) Vacances$Toussaint = seq(from = as.Date("21/10/2017", format="%d/%m/%Y"), to=as.Date("05/11/2017", format="%d/%m/%Y"), by=1) Vacances$All=Reduce(append,Vacances) |
Après, un peu de nettoyage est nécessaire, avec des gares en double (par exemple quand passe a la fois le RER et le métro), et pour recuperer leur localisation spatiale (latitude et longitude)
gares$NOM_LONG = as.character(gares$NOM_LONG) DD = (gares$NOM_LONG[duplicated(gares$NOM_LONG)]) i = (gares$NOM_LONG %in% DD) & gares$MODE_=="Metro" gares$NOM_LONG[i] = paste(gares$NOM_LONG[i],"M", sep="-") i = (gares$NOM_LONG %in% DD) & gares$MODE_=="RER" gares$NOM_LONG[i] = paste(gares$NOM_LONG[i],"R", sep="-") gares$NOM_LONG=factor(gares$NOM_LONG) a=as.character(gares$Geo.Point) gares$Y=as.numeric(str_extract_all(a,"^[0-9]+.[0-9]+")) gares$X=as.numeric(str_extract_all(a,"[0-9]+.[0-9]+$")) |
On compte ensuite les nombres de validation de tickets, par gare
Manip_nbvalid = function(Data,DD,gares) { i=grep("^[a-zA-Z]+",as.character(Data$NB_VALD)) Data$NB_VALD[i]=as.integer(5) i=is.na(Data$NB_VALD) Data$NB_VALD[i]=as.integer(5) Data$LIBELLE_ARRET=as.character(Data$LIBELLE_ARRET) i=(Data$LIBELLE_ARRET %in% DD) & Data$CODE_STIF_TRNS=="100" Data$LIBELLE_ARRET[i]=paste(Data$LIBELLE_ARRET[i],"M", sep="-") i=(Data$LIBELLE_ARRET %in% DD) & Data$CODE_STIF_TRNS=="800" Data$LIBELLE_ARRET[i]=paste(Data$LIBELLE_ARRET[i],"R", sep="-") for (i in seq(1,nrow(Cooresp))) { Data$LIBELLE_ARRET=gsub(as.character(Cooresp$nbval[i]),as.character(Cooresp$gares[i]),Data$LIBELLE_ARRET) } gares$NOM_LONG=as.character(gares$NOM_LONG) Data=dplyr::left_join(Data,gares[,c("NOM_LONG","X","Y")],by=c("LIBELLE_ARRET"="NOM_LONG")) Data=Data[is.na(Data$CODE_STIF_ARRET)==FALSE,] Data=Data[Data$CODE_STIF_ARRET!="ND",] Data$NB_VALD=as.integer(as.character(Data$NB_VALD)) Data$JOUR=as.Date(Data$JOUR) Data$CODE_STIF_TRNS=factor(Data$CODE_STIF_TRNS) Data$CODE_STIF_RES=factor(Data$CODE_STIF_RES) Data$CODE_STIF_ARRET=factor(Data$CODE_STIF_ARRET) Data$LIBELLE_ARRET=factor(Data$LIBELLE_ARRET) Data$ID_REFA_LDA=factor(Data$ID_REFA_LDA) Data$CATEGORIE_TITRE=factor(Data$CATEGORIE_TITRE) Data$JOURSEM=weekdays(Data$JOUR) return(Data) } nbvalid=lapply(nbvalid, Manip_nbvalid,DD=DD,gares=gares) |
On a ainsi tous les comptages, pour toutes les gares. On fait ensuite un découpage par tranche horaire
Manip_dist_pers = function(DataFrame) { DataFrame=DataFrame[(DataFrame$TRNC_HORR_60)!="ND",] DataFrame$TRNC_HORR_60=factor(DataFrame$TRNC_HORR_60, levels = c("0H-1H", "1H-2H", "2H-3H", "3H-4H", "4H-5H", "5H-6H", "6H-7H", "7H-8H", "8H-9H", "9H-10H", "10H-11H", "11H-12H", "12H-13H", "13H-14H", "14H-15H", "15H-16H", "16H-17H", "17H-18H", "18H-19H", "19H-20H", "20H-21H", "21H-22H", "22H-23H", "23H-0H")) DataFrame=DataFrame[(DataFrame$CODE_STIF_ARRET)!="ND",] DataFrame$CODE_STIF_ARRET=factor(DataFrame$CODE_STIF_ARRET) DataFrame$TRANCHE=str_extract(as.character(DataFrame$TRNC_HORR_60),"^[0-9]{1,2}") return(DataFrame) } distr_pers=lapply(distr_pers, Manip_dist_pers) |
On peut ensuite recuperer la distribution des validation, par jour
distr_JOURV=list() distr_JOURV$S1 = nbvalid$S1 %>% group_by(JOUR, JOURSEM,CATEGORIE_TITRE) %>% summarise(NB_VALD=sum(NB_VALD)) distr_JOURV$S2 = nbvalid$S2 %>% group_by(JOUR, JOURSEM,CATEGORIE_TITRE) %>% summarise(NB_VALD=sum(NB_VALD)) distr_JOURV$Y=rbind(distr_JOURV$S1,distr_JOURV$S2) distr_JOUR=list() distr_JOUR$S1 = nbvalid$S1 %>% group_by(JOUR, JOURSEM) %>% summarise(NB_VALD=sum(NB_VALD)) distr_JOUR$S2 = nbvalid$S2 %>% group_by(JOUR, JOURSEM) %>% summarise(NB_VALD=sum(NB_VALD)) distr_JOUR$Y=rbind(distr_JOUR$S1,distr_JOUR$S2) distr_JOUR_Station=list() distr_JOUR_Station$S1 = nbvalid$S1 %>% group_by(JOUR, JOURSEM,CODE_STIF_ARRET,LIBELLE_ARRET) %>% summarise(NB_VALD=sum(NB_VALD), X=max(X), Y=max(Y)) distr_JOUR_Station$S2 = nbvalid$S2 %>% group_by(JOUR, JOURSEM,CODE_STIF_ARRET,LIBELLE_ARRET) %>% summarise(NB_VALD=sum(NB_VALD), X=max(X), Y=max(Y)) Manip_dist_Jour = function(DataFrame) { DataFrame$JOURSEM=factor(DataFrame$JOURSEM,levels = c("lundi","mardi","mercredi","jeudi","vendredi","samedi","dimanche")) DataFrame$TypeJ=character(nrow(DataFrame)) DataFrame$TypeJ[DataFrame$JOUR %in% Vacances$Ete]="Ete" DataFrame$TypeJ[DataFrame$JOUR %in% Vacances$Noel]="Noel" DataFrame$TypeJ[DataFrame$JOUR %in% Vacances$Ski]="Ski" DataFrame$TypeJ[DataFrame$JOUR %in% Vacances$Printemps]="Printemps" DataFrame$TypeJ[DataFrame$JOUR %in% Vacances$Toussaint]="Toussaint" DataFrame$TypeJ[DataFrame$JOUR %in% Vacances$All == FALSE]="HorsVacances" DataFrame$CAT_JOUR=character(nrow(DataFrame)) DFr=list() ii=(DataFrame$JOURSEM!="samedi" & DataFrame$JOURSEM!="dimanche") & DataFrame$TypeJ!="HorsVacances" DataFrame$CAT_JOUR[ii]="JOVS" DFr$JOVS$Data = DataFrame[ii,] DFr$JOVS$Nom="Jours ouvrés Vacances Scolaires" ii=(DataFrame$JOURSEM!="samedi" & DataFrame$JOURSEM!="dimanche") & DataFrame$TypeJ=="HorsVacances" DataFrame$CAT_JOUR[ii]="JOHV" DFr$JOHV$Data = DataFrame[ii,] DFr$JOHV$Nom="Jours ouvés Hors Vacances Scolaires" ii=DataFrame$JOURSEM=="samedi" & DataFrame$TypeJ!="HorsVacances" DataFrame$CAT_JOUR[ii]="SAVS" DFr$SAVS$Data = DataFrame[ii,] DFr$SAVS$Nom="Samedi VS" ii=DataFrame$JOURSEM=="samedi" & DataFrame$TypeJ=="HorsVacances" DataFrame$CAT_JOUR[ii]="SAHV" DFr$SAHV$Data = DataFrame[ii,] DFr$SAHV$Nom="Samedi HV" ii=DataFrame$JOURSEM=="dimanche" DataFrame$CAT_JOUR[ii]="DIJFP" DFr$DIJFP$Data = DataFrame[ii,] DFr$DIJFP$Nom="Dimanche" return(list("TypeJ"=DFr,"Distr"=DataFrame)) } res=Manip_dist_Jour(distr_JOUR$Y) distr_TypeJ=res$TypeJ distr_JOUR$Y=res$Distr res=Manip_dist_Jour(distr_JOURV$Y) distr_TypeJV=res$TypeJ distr_TypeJ_Station=list() res=Manip_dist_Jour(distr_JOUR_Station$S1) distr_TypeJ_Station$S1=res$TypeJ distr_JOUR_Station$S1=res$Distr res=Manip_dist_Jour(distr_JOUR_Station$S2) distr_TypeJ_Station$S2=res$TypeJ distr_JOUR_Station$S2=res$Distr rm(res) |
On peut alors tracer toutes sortes de graphiques, par exemple le nombre de validations, par jour, entre le 1er janvier et le 31 décembre, en fonction du jour de la semaine.
g0 = ggplot(distr_JOUR$Y, aes(x=JOUR, y=NB_VALD, color = JOURSEM)) + geom_point() g0 = g0 + labs(title="Nombres de validations chaque jours de 2017", x="Date", y="Nombre de validations") g0 |
On peut voir la très forte baisse les jours de semaine pendant les vacances d’été. Au lieu de regarder sur l’année, on peut regarder sur la journée
Fct_FqH = function(DataFrame,distr_pers) { DataFrame=dplyr::full_join(DataFrame,distr_pers[,c("CAT_JOUR","CODE_STIF_ARRET","pourc_validations","TRANCHE","TRNC_HORR_60")],by=c("CODE_STIF_ARRET"="CODE_STIF_ARRET","CAT_JOUR"="CAT_JOUR")) DataFrame$NB_VALD=DataFrame$NB_VALD*DataFrame$pourc_validations return(DataFrame) } distr_JOUR_Station$S1=Fct_FqH(distr_JOUR_Station$S1, distr_pers$S1) distr_JOUR_Station$S2=Fct_FqH(distr_JOUR_Station$S2, distr_pers$S2) distr_JOUR_Station$Y=rbind(distr_JOUR_Station$S1,distr_JOUR_Station$S2) distr_JOUR_Station$Y=distr_JOUR_Station$Y[is.na(distr_JOUR_Station$Y$NB_VALD)==FALSE,] |
On peut alors faire un graphique, en fonction de la tranche horaire, pour certaines périodes, par exemple en dehors de vacances scolaires, en semaine (par heure, on a ici un boxplot)
Graphique_HOR = function(DataFrame,TypeJ,NomJ) { # Graphique de la distribution de l'affluence par tranche horaire et type de jours g1 = ggplot(DataFrame[DataFrame$CAT_JOUR==TypeJ,], aes(x=TRNC_HORR_60, y=pourc_validations, color = TRNC_HORR_60,las=2)) + geom_boxplot() + ylim(c(0,100)) g1 = g1 + labs(title=paste(c("Distribution des validations par tranche horaire ",NomJ), sep="", collapse = ""), x="Jours", y="Nombre de validations") + theme(axis.text.x= element_text(size = 8, angle = 45)) g1 } Graphique_HOR(distr_JOUR_Station$Y,"JOHV","Jours ouvrés Hors Vacances Scolaires") |
ou bien le samedi
Graphique_HOR(distr_JOUR_Station$Y,"SAHV","Samedi Hors Vacances Scolaires") |
On peut tenter un peu de cartographie. Comme nombre de métros/bus, dans le monde, on a souvent uniquement accès aux nœuds d’entrée dans le réseau (et pas aux nœuds de sortie). Mais ça reste intéressant, et très informatif
get_Paris1 = get_map(c(2.3448688,48.8613029), zoom = 11) Paris1 = ggmap(get_Paris1) |
Par gare, et par heure, on peut regarder le nombre de validations de tickets
Median_Valid = distr_JOUR_Station$Y %>% group_by(CAT_JOUR, LIBELLE_ARRET, X, Y) %>% summarise(NB_VALD=median(NB_VALD)) Median_Valid_Station = distr_JOUR_Station$Y %>% group_by(CAT_JOUR, TRNC_HORR_60,LIBELLE_ARRET, X, Y) %>% summarise(NB_VALD=median(NB_VALD)) Carte_Densite = function(Nom,Carte,TypeJ,HOR,DataFrame) { if (HOR=="") { ii=DataFrame$CAT_JOUR==TypeJ NomSave=paste("Densité des validations",Nom,TypeJ) } else { ii=DataFrame$CAT_JOUR==TypeJ & DataFrame$TRNC_HORR_60==HOR NomSave=paste("Densité des validations",Nom,TypeJ,HOR) } U=DataFrame[ii,] n=round(log10(median(U$NB_VALD)))-1 n=max(1,10^n) Nb_Repete_Stations=ceiling(U$NB_VALD/n) U$Size_Stations=U$NB_VALD/max(U$NB_VALD) Z=U[rep(1:nrow(U),Nb_Repete_Stations),] Carte_A= Carte + geom_point(aes(x=X,y=Y),data=Z,col="coral", size=10*Z$Size_Stations) + geom_density2d(data = Z, aes(x=X,y=Y), size = 0.5) + stat_density2d(data = Z, aes(x=X,y=Y,fill = ..level.., alpha = ..level..),size = 0.01, bins = 16, geom = "polygon") + scale_fill_gradient(low = "chartreuse", high = "red",guide = FALSE) + scale_alpha(range = c(0, 0.3), guide = FALSE) + ggtitle(NomSave) + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x= element_blank(), axis.text.y = element_blank()) suppressWarnings(print(Carte_A)) } |
Par exemple, si on regarde les points de validations de tickets entre 5 et 6 heures du matin, on obtient
L=levels(Median_Valid_Station$TRNC_HORR_60) Carte_Densite("dans la petite ceinture",Paris1,"JOHV",L[6],Median_Valid_Station) |
avec beaucoup de ville dans la banlieue proche. Plus tard en journée, entre 11 heures et midi, les gares de validation sont davantage dans le cœur de Paris, avec la Défense a gauche et Saint-Denis au nord
Carte_Densite("dans la petite ceinture",Paris1,"JOHV",L[12],Median_Valid_Station) |
En fin de journée, c’est Paris et surtout la Défense qui ressortent
Carte_Densite("dans la petite ceinture",Paris1,"JOHV",L[19],Median_Valid_Station) |
Amusant, non ?
