{"word_count": 933, "character_count": 13375}
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 lignenbvalid = list()download.file("https://opendata.stif.info/explore/dataset/emplacement-des-gares-idf-data-generalisee/download/?format=csv&#38;timezone=Europe/Berlin&#38;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&#38;timezone=Europe/Berlin&#38;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&#38;timezone=Europe/Berlin&#38;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&#38;timezone=Europe/Berlin&#38;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&#38;timezone=Europe/Berlin&#38;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 2017Vacances = 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) &#38; gares$MODE_=="Metro"gares$NOM_LONG[i] = paste(gares$NOM_LONG[i],"M", sep="-")i = (gares$NOM_LONG %in% DD) &#38; 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 gareManip_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) &#38; Data$CODE_STIF_TRNS=="100"  Data$LIBELLE_ARRET[i]=paste(Data$LIBELLE_ARRET[i],"M", sep="-")  i=(Data$LIBELLE_ARRET %in% DD) &#38; 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 horaireManip_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 jourdistr_JOURV=list()distr_JOURV$S1 = nbvalid$S1 %&#62;% group_by(JOUR, JOURSEM,CATEGORIE_TITRE) %&#62;% summarise(NB_VALD=sum(NB_VALD))distr_JOURV$S2 = nbvalid$S2 %&#62;% group_by(JOUR, JOURSEM,CATEGORIE_TITRE) %&#62;% summarise(NB_VALD=sum(NB_VALD))distr_JOURV$Y=rbind(distr_JOURV$S1,distr_JOURV$S2)distr_JOUR=list()distr_JOUR$S1 = nbvalid$S1 %&#62;% group_by(JOUR, JOURSEM) %&#62;% summarise(NB_VALD=sum(NB_VALD))distr_JOUR$S2 = nbvalid$S2 %&#62;% group_by(JOUR, JOURSEM) %&#62;% 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 %&#62;% group_by(JOUR, JOURSEM,CODE_STIF_ARRET,LIBELLE_ARRET) %&#62;% summarise(NB_VALD=sum(NB_VALD), X=max(X), Y=max(Y))distr_JOUR_Station$S2 = nbvalid$S2 %&#62;% group_by(JOUR, JOURSEM,CODE_STIF_ARRET,LIBELLE_ARRET) %&#62;% 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" &#38; DataFrame$JOURSEM!="dimanche") &#38; 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" &#38; DataFrame$JOURSEM!="dimanche") &#38; 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" &#38; DataFrame$TypeJ!="HorsVacances"  DataFrame$CAT_JOUR[ii]="SAVS"  DFr$SAVS$Data = DataFrame[ii,]  DFr$SAVS$Nom="Samedi VS"  ii=DataFrame$JOURSEM=="samedi" &#38; 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$TypeJdistr_JOUR$Y=res$Distrres=Manip_dist_Jour(distr_JOURV$Y)distr_TypeJV=res$TypeJdistr_TypeJ_Station=list()res=Manip_dist_Jour(distr_JOUR_Station$S1)distr_TypeJ_Station$S1=res$TypeJdistr_JOUR_Station$S1=res$Distrres=Manip_dist_Jour(distr_JOUR_Station$S2)distr_TypeJ_Station$S2=res$TypeJdistr_JOUR_Station$S2=res$Distrrm(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")g0On 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éeFct_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 samediGraphique_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 informatifget_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 ticketsMedian_Valid = distr_JOUR_Station$Y %&#62;% group_by(CAT_JOUR, LIBELLE_ARRET, X, Y) %&#62;% summarise(NB_VALD=median(NB_VALD))Median_Valid_Station = distr_JOUR_Station$Y %&#62;% group_by(CAT_JOUR, TRNC_HORR_60,LIBELLE_ARRET, X, Y) %&#62;% 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 &#38; 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 obtientL=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 nordCarte_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 ressortentCarte_Densite("dans la petite ceinture",Paris1,"JOHV",L[19],Median_Valid_Station)Amusant, non ?