Ce lundi, je participerais (en visio, je serai rentré sur Montréal) à la formation Data Science de l’Institut des Actuaires, pour parler d’équité. Les slides sont en ligne ici.
Tag Archives: DSA
Acheter un billet de train (pas trop cher)
Hier, je suis tombé sur article qui discutait des prix des billets de train, en France (et du prix très élevé, a certaines dates, genre pendant les vacances d’hiver). Tous ceux qui ont l’habitude de prendre le train savent que le prix que l’on paye dépend du moment ou on achète le billet (et de la souplesse que l’on peut avoir sur l’heure (voire la date) du trajet). Cet été, dans le cadre du projet de la formation Data Science pour l’Actuariat, pour mon cours, Pierre proposait de moissonner le site https://www.oui.sncf/ pour suivre un peu l’évolution du prix des billets.
Le principal soucis est que le site https://www.oui.sncf/ s’appuie sur du javascript pour les formulaires de saisie et pour l’affichage des résultats, ce qui empêche l’utilisation classique du package rvest, par exemple. J’avais évoqué dans un autre billet l’utilisation de wdman, pour scraper le site des incendies de forets. Ici, Pierre proposait de passer par casperjs, et je vais reprendre un peu ici sa stratégie:
- on va utiliser casperjs, un émulateur de navigateur écrit en javascript. Il permet d’émuler un véritable navigateur (même moteur que google chrome) et de résoudre le javascript intégrés à la page
- on va utiliser un petit code en bash pour lancer le code
Pour le code en bash, c’est juste que je suis sous mac et linux. Sous mac, on peut faire pas mal de choses en bash… Tout passe par des variables, que l’on peut définir, et afficher, par exemple
qui nous donne l’heure. Je peux aussi définir une variable, et l’incrémenter (pratique pour faire des boucles)
Plus intéressant, on peut planifier des taches. Pour ça, on tape
qui va ouvrir un éditeur,
Je demande ici de lancer un code R, tous les heures, a 13:50, 14:50, 15:50, etc. Pour demander tous les jours a 13:50, je tape
on va ensuite sauver l’instruction
On voit que la commande sera lancée, tous les jours : elle est dans la liste des taches a faire
notons qu’on peut lancer un code tout en passant des arguments : ici je lui dis quel objet manipuler (le premier argument) et le second sert pour créer un fichier (et pour le nommer).
Bref, c’est assez facile de lancer automatiquement des codes, pour scraper. Sous windows, on passe par le planificateur de taches. Un premier script permet d’extraire les trains ouverts dans la journée,
et le second, les trains ouverts depuis plus de 24 heures,
Ensuite, on va utiliser http://docs.casperjs.org/en/latest/ pour coder notre émulateur de navigateur internet (le code est ici en ligne).
On va ainsi créer plein de fichiers, contenant les informations que l’on veut ! Je vais passer un peu le retraitement, et juste présenter les informations qu’on peut en tirer. En particulier, Pierre avait stocke les données entre mars et juin dernier.
Ici, on est juste sur quelques trajets de grande ligne,
library(readr) library(rgdal) nomFichier = tempfile(fileext = ".zip") download.file("https://freakonometrics.free.fr/CarteFrance.zip", destfile = nomFichier, mode = "wb") unzip(zipfile = nomFichier, exdir = getwd()) download.file("https://freakonometrics.free.fr/LgTroncons.csv", destfile = "LgTroncons.csv", mode = "wb") download.file("https://freakonometrics.free.fr/CoordVilles.csv", destfile = "CoordVilles.csv", mode = "wb") fra0 = readOGR(dsn = paste(getwd(), "/CarteFrance", sep = ""), layer = "gadm36_FRA_0", verbose = F) LgTroncons = read_delim("LgTroncons.csv",";", escape_double = FALSE,locale = locale(decimal_mark = ","), trim_ws = TRUE) CoordVilles = read_delim("CoordVilles.csv",";", escape_double = FALSE,locale = locale(decimal_mark = ","), trim_ws = TRUE) NomsVilles = CoordVilles[CoordVilles$NOM_A_AFFICHER==1,] library(ggplot2) fr_df = fortify(fra0) ggp = ggplot() + geom_polygon(data=fr_df, aes(long, lat,group = group), fill = "#3A8EBA") ggp = ggp + geom_path(data = LgTroncons, aes(x = LONG, y = LAT, group = ID_TRONCON), colour= "#CC5500", lineend = "round", size=3) + geom_path(data = LgTroncons, aes(x=LONG, y=LAT, group = ID_TRONCON), colour="white", lineend = "round", size=1.75) ggp = ggp + geom_point(data = NomsVilles, aes(x=LONG, y=LAT), colour = "blue", fill = "white", shape=21, size = NomsVilles$PT_SIZE, stroke = NomsVilles$PT_STROKE) + theme_void() ggp = ggp + geom_text(data = NomsVilles, aes(x=LONG, y=LAT, label=NOM),hjust = NomsVilles$H_AJUST, vjust=NomsVilles$V_AJUST, colour = "white", fontface = "bold", size =3.25)+coord_fixed(1.47) ggp = ggp + ggtitle("Représentation des trajets étudiées") + theme(plot.title = element_text(hjust = 0.5, face="bold")) print(ggp) |
On va travailler sur les trajets suivants,
Comparons ici les billets pour un trajet Pars-Rennes, a partir des informations moissonnées pendant 3 mois, le vendredi soir, en particulier 2 vendredi de juin 2018 (les 15 et 22 juin). Pour ces deux jours, il y a eu 6 trains, entre 17 et 20 heures. Pour le 15 juin, les trois premiers ont commencé avec un prix de 45 €. Pour le 22 juin, le premier a commence a 45 €, mais les deux suivant on été lancés a 33 €. Assez rapidement, les prix sont monte a 45 €.
On peut regarder l’evolution du prix
Si on regarde plusieurs destination, on observe des comportement très différentes,
- pour Le Mans, les prix montent très vite, commençant a 15 €, montant a 18 € 10 heures après le lancement, 21 € le lendemdemain, 27 € au bout d’une semaine. En un mois, les prix ont presque doublé.
- pour Rennes, on observe une evolution similaire, passant de 20 € a 25 € en quelques heures, et 40 € deux semaines après !
- pour Toulouse en revanche, le prix initial est plus haut, 43 €, monte de 3 € en 10 heures, 6 € en 16 heures, pour rester a 49 €
Mais pour toutes les destinations, les prix sont croissants.
soit graphiquement
On peut aussi faire une carte. Si on regarde les prix a l’ouverture, Lille, Le Mans et Rennes sont peu chères.
Et les plus fortes variations sur 10 heures sont observées sur Nantes et Bordeaux.
Amusant non ?
Les transports publics parisiens
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 ?
Analyse des résultats au baccalauréat des séries générales
Pour continuer sur les manipulation de données publiques, je voulais m’inspirer du projet de Cédric, de la formation Data Science pour l’Actuariat sur les résultats au baccalauréat. Les données nécessaires à cette étude sont disponibles sur plusieurs sites,
- Ministère de l’éducation nationale
- Bases INSEE
- NosDonnees.fr
Il ne s’agit aucunement d’une analyse poussée des résultats, juste un peu de visualisation, sans aucune autre prétention ! Ah oui, même si on ne va pas faire de carte (je les trouve peu lisibles) on va quand même utiliser les données spatiales : les établissements scolaires sont géolocalises, et on peut obtenir des informations locales, sur le taux de chômage, ou le revenu médian. Et faire des graphiques.
Ce préambule passé, on peut commencer.
library(dplyr) library(readxl) library(sp) library(ggmap) library(raster) library(leaflet) library(DT) library(cowplot) library(gstat) library(tmap) |
On va commencer par récupérer par établissement, les résultats au bac.
url_resultat_etab = "https://data.education.gouv.fr/explore/dataset/fr-en-indicateurs-de-resultat-des-lycees-denseignement-general-et-technologique/download/?format=csv&timezone=Europe/Berlin&use_labels_for_header=true" download.file(url_resultat_etab,destfile = paste0(librairie,"import_resultat_etab.csv"), method="curl") df_resultat_etab = read.csv("import_resultat_etab.csv",header=TRUE, sep= ";", encoding="UTF-8") |
Comme bien souvent avec les données des administrations françaises, on a souvent des soucis de typographie. Pour simplifier, on va supprimer les accents, et uniformiser un peu les noms
MiseEnForme_Colonnes = function(text) { text <- gsub("è", "e", text) text <- gsub("é", "e", text) text <- gsub("_", ".", text) text <- gsub("serie.", "", text) text <- gsub("Effectif.Presents.", "Effectif.", text) text <- gsub("Taux.","Tx.",text) text <- gsub("Brut.de.Reussite.", "Admis.Etab.", text) text <- gsub("Reussite.Attendu.", "Admis.", text) text <- gsub("brut", "Etab", text) text <- gsub("attendu", "Academie", text) text <- gsub("toutes.", "TOTAL", text) text <- gsub("Total.", "TOTAL", text) text <- gsub("..Etablissement", ".Etab", text) text <- gsub("Pourcentage", "Tx", text) return(text) } for(i in 1:ncol(df_resultat_etab)){ colnames(df_resultat_etab)[i] <- MiseEnForme_Colonnes(names(df_resultat_etab)[i]) } |
On va ensuite supprimer les départements et régions d’outre-mer,
df_resultat_etab = df_resultat_etab[-which(toupper(df_resultat_etab$Departement) %in% c("GUADELOUPE","MAYOTTE","MARTINIQUE","REUNION","GUYANE")),] |
récupérer les noms des colonnes
Colonnes = colnames(df_resultat_etab) |
et comme on s’intéresse aux premières variables
Colonnes_Generiques = Colonnes[1:8] |
on les recupere, pour construire quelques statistiques pour colonnes relatives aux séries L, ES et S
Colonnes_Series = Colonnes[grepl("([a-zA-Z]*?.)*\\.S$|([a-zA-Z]*?.)*\\.ES$|([a-zA-Z]*?.)*\\.L$|([a-zA-Z]*?.)*\\.TOTAL$",Colonnes)] |
Et on finit avec les autres
Colonnes_Autres = Colonnes[grepl("(Tx.Bacheliers.*)|(Tx.acces.*)|(Effectif.de.*)|(libelle.region)|(code.region)|(element)",Colonnes)] df_resultat_etab = cbind(df_resultat_etab[Colonnes_Generiques],df_resultat_etab[Colonnes_Series],df_resultat_etab[Colonnes_Autres]) |
On peut aussi localiser les établissements
url_carto_etab <- "https://www.data.gouv.fr/s/resources/adresse-et-geolocalisation-des-etablissements-denseignement-du-premier-et-second-degres/20160526-143453/DEPP-etab-1D2D.csv" download.file(url_carto_etab,destfile=paste0(librairie,"import_carto_etab.csv")) df_carto_etab = read.csv2("import_carto_etab.csv",header=TRUE |
On récupère ici la géolocalisation de 66556 établissements nationaux ! On peut croiser avec des données socio-économiques des communes
nom_base_emploi = "base-cc-emploi-pop-act-2014" url_baseemploi_popactive = paste0("https://www.insee.fr/fr/statistiques/fichier/2862207/",nom_base_emploi,".zip") download.file(url_baseemploi_popactive,destfile=paste0(librairie,nom_base_emploi,".zip")) unzip(paste0(nom_base_emploi,".zip"),overwrite = TRUE) df_base_emploi_source = read_excel(paste0(nom_base_emploi,".xls"),sheet="COM_2014",skip=5) |
On va exclure les territoires d’outre-mer ici
df_base_emploi_source <- df_base_emploi_source[-which(df_base_emploi_source$DEP %in% c("971","972","973","974","975")),] df_base_emploi_colonnes = c("CODGEO","P14_POP1564","P14_H1564","P14_F1564","P14_ACT1564","P14_ACTOCC1564","P14_CHOM1564","P14_INACT1564", "P14_ETUD1564", "P14_RETR1564", "P14_AINACT1564", "P14_HCHOM1524", "P14_FCHOM1524", "C14_ACT1564","C14_ACT1564_CS1","C14_ACT1564_CS2","C14_ACT1564_CS3","C14_ACT1564_CS4","C14_ACT1564_CS5","C14_ACT1564_CS6","P14_POP15P") df_base_emploi = df_base_emploi_source[,names(df_base_emploi_source) %in% df_base_emploi_colonnes] |
et corriger les soucis classiques de la Corse,
MiseEnForme_CodeGeo = function(text) { text <- gsub("2A", "20", text) text <- gsub("2B", "20", text) } df_base_emploi$CODGEO = MiseEnForme_CodeGeo(df_base_emploi$CODGEO) |
On peut aussi utiliser des données de revenus, par communes
nom_base_revenus = "indic-struct-distrib-revenu-2014-COMMUNES" url_baserevenus = paste0("https://www.insee.fr/fr/statistiques/fichier/3126151/",nom_base_revenus,".zip") download.file(url_baserevenus,destfile=paste0(librairie,nom_base_revenus,".zip")) unzip(paste0(nom_base_revenus,".zip"),overwrite = TRUE) df_base_revenus = read_excel("FILO_DISP_COM.xls",sheet="ENSEMBLE",skip=5)[,c(1,4,7)] df_base_revenus$CODGEO = MiseEnForme_CodeGeo(df_base_revenus$CODGEO) |
On recupere des donnees spatiales relatives aux communes
url_geoloc_communes = "http://www.nosdonnees.fr/wiki/images/b/b5/EUCircos_Regions_departements_circonscriptions_communes_gps.csv.gz" download.file(url_geoloc_communes,destfile=paste0(librairie,"geoloc_communes.csv.gz")) df_geoloc_communes = read.csv2(gzfile("geoloc_communes.csv.gz"),header=TRUE, stringsAsFactors = FALSE,encoding="UTF-8") |
et comme toujours, un peu de corrections s’imposent
df_geoloc_communes = df_geoloc_communes[-which(df_geoloc_communes$numéro_département %in% c("971","972","973","974","975")),] df_geoloc_communes = df_geoloc_communes[,names(df_geoloc_communes) %in% c("code_insee","latitude","longitude","codes_postaux")] df_geoloc_communes_nb <- nrow(df_geoloc_communes) |
On va ensuite creer une fonction de remplacement des valeurs manquantes, et de correction des séparateurs décimaux
MiseEnForme_CoordonneesGeo = function(valeur){ pretraitement = ifelse(as.character(valeur)=="-","0",as.character(valeur)) traitement = as.numeric(ifelse(pretraitement==".","0",gsub(pattern=",",replacement=".",pretraitement))) return(traitement) } df_geoloc_communes$latitude = MiseEnForme_CoordonneesGeo(df_geoloc_communes$latitude) df_geoloc_communes$longitude = MiseEnForme_CoordonneesGeo(df_geoloc_communes$longitude) |
On passe ensuite a l’élimination des lignes en double
df_geoloc_communes = unique(df_geoloc_communes) |
On va ensuite changer les noms des colonnes pour harmoniser avec les autres bases
names(df_geoloc_communes) = c("Codes_Postaux","CODGEO","coordonnee_y","coordonnee_x") |
On peut ensuite rechercher les lignes en double sur les codes insee
liste_CODGEO2 = aggregate(x=df_geoloc_communes$Codes_Postaux,by=list(df_geoloc_communes$CODGEO),FUN="length") list_geoloc_communes_CODGEO2 = liste_CODGEO2[liste_CODGEO2$x>1,1] df_geoloc_communes_CODGEO2 = df_geoloc_communes[df_geoloc_communes$CODGEO %in% list_geoloc_communes_CODGEO2,1:2] |
Ici, un correction manuelle s’impose pour 4 configurations : les données propres à Lyon, Paris et Marseille ne sont pas géolocalisées, les données propres à la ville de Laguépie sont géolocalisées en doubles
df_geoloc_communes_propre = df_geoloc_communes[!df_geoloc_communes$CODGEO %in% list_geoloc_communes_CODGEO2,] df_geoloc_communes_corrige = data.frame(Codes_Postaux=c("13001","69001","75001","82250"), CODGEO=c("13055","69123","75056","82088"), coordonnee_y=c(43.3,45.75,48.85,44.15), coordonnee_x=c(5.4,4.85,2.31,1.97)) df_geoloc_communes = rbind(df_geoloc_communes_propre,df_geoloc_communes_corrige) |
On peut enfin fusionner les bases
df_etab = merge(df_resultat_etab,df_carto_etab,by="Cod.Etab") |
Certains établissements ne peuvent être géolocalisés pour certaines années
df_etab_total_nongeolocalises <- df_resultat_etab[!df_resultat_etab$Cod.Etab %in% df_carto_etab$Cod.Etab,] |
Comme l’étude ne porte que sur les seuls lycées d’enseignement général et technologique, le dataframe est réduit aux observations relatives d’une part aux lycées, d’autre part aux établissements d’enseignement polyvalent, général ou général et technologique.
df_etab = df_etab[grep("LYCÉE",toupper(df_etab$nature_uai_libe)),] df_etab = df_etab[grep("GÉNÉRAL|POLYVALENT",toupper(df_etab$nature_uai_libe)),] df_etab_nongeolocalises = df_etab[df_etab$Cod.Etab %in% df_etab_total_nongeolocalises$Cod.Etab,] df_etab_geolocalise = df_etab[!is.na(df_etab$coordonnee_x),] df_etab_geolocalise = df_etab_geolocalise[!is.na(df_etab_geolocalise$coordonnee_y),] |
Enfin, on va convertir les code géographiques (ici reconnus comme facteur) en chaines de caractères (de 5 caractères) pour pouvoir fusionner les tables
ConvertCODGEO = function(code) { if(is.character(code)) { code_character = ifelse(nchar(code)<5, paste0("0",code), code) return(code_character) } else if(is.factor(code)){ code_character = ifelse(code<10000, paste0("0",as.numeric(as.character(code))), as.numeric(as.character(code))) return(code_character) } else if(is.numeric(code)){ code_character = ifelse(code<10000, paste0("0",code), as.character(code)) return(code_character) } } df_etab_geolocalise$Code.commune = ConvertCODGEO(df_etab_geolocalise$Code.commune) df_etab_geolocalise$Secteur.Public.Prive = sapply(df_etab_geolocalise$Secteur.Public.Prive,function(nature) {ifelse(nature=="PU","Lycées Publics","Lycées Privés")}) |
On conservation alors les établissements dont la commune n’est pas manquante
df_etab_geolocalise = df_etab_geolocalise[!is.na(df_etab_geolocalise$Code.commune),] |
Pour finir, on va creer une base, pour ensuite faire une graphique
tbl_etab_nature_res_source = df_etab_geolocalise[,c(3,8,9,10,11,13,14,15)] for(i in c(6,7,8)){ temp = tbl_etab_nature_res_source[!is.na(tbl_etab_nature_res_source[i]),c(1,2,i-3,i)] temp$Serie = ifelse(i==6,"L",ifelse(i==7,"ES","S")) names(temp)[2:4] = c("Nature","Effectif","Tx.Admis") if(i==6){ tbl_etab_nature_result = temp } else{ tbl_etab_nature_result = rbind(tbl_etab_nature_result,temp) } } graph = ggplot(tbl_etab_nature_result,aes(x=Effectif,y=Tx.Admis,colour=factor(Annee))) graph = graph + geom_point(alpha=0.45) graph = graph + facet_grid(Serie~Nature) graph = graph + xlab("Effectifs de l'établissement en terminale (par série)") + ylab("Taux d'admission (%)") graph = graph + scale_color_discrete(name="Année des\nrésultats") graph = graph + theme(legend.title = element_text(size=9,face="bold"), legend.text = element_text(size=9), strip.background = element_rect(colour="black", fill="gray95"), panel.border = element_rect(linetype = "solid"), panel.grid.major = element_line(colour = "gray75",linetype = "dashed"), panel.grid.minor = element_line(colour = "gray95",linetype = "dashed"), axis.title.x = element_text(size=9, face="bold"), axis.text.x = element_text(size=8), axis.title.y = element_text(size=9, face="bold"), axis.text.y = element_text(size=8)) graph |
On a ici l’evolution des resultats en fonction de la taille des etablissements.
df_communes_CorrNaN = df_communes_Corr[which(!df_communes_Corr$TxChomage == "NaN" & !df_communes_Corr$TxCadres == "NaN" & !df_communes_Corr$TxOuvriers == "NaN" & !df_communes_Corr$NbPopulation == "NaN" & !df_communes_Corr$TxSenior == "NaN" & !df_communes_Corr$RevenusMedians == "NaN"),] df_communes_sp = SpatialPointsDataFrame(coords = df_communes_CorrNaN[, c("coordonnee_x", "coordonnee_y")], data = df_communes_CorrNaN) Grille = as.data.frame(makegrid(df_communes_sp, nsig=2, cellsize = 0.1)) names(Grille) = c("X", "Y") coordinates(Grille) = c("X", "Y") gridded(Grille) = TRUE fullgrid(Grille) = TRUE proj4string(Grille) = proj4string(df_communes_sp) |
On peut ensuite faire du krigeage, histoire de lisser un peu nos donnees de chomage et de revenu
df_communes_sp.TxChomage = krige(TxChomage ~ 1, df_communes_sp, Grille, nmax=1) df_communes_sp.RevenusMedians = krige(RevenusMedians ~ 1, df_communes_sp, Grille, nmax=1) sp_lycee_WGS84@data$TxChomage = extract(R.TxChomage,sp_lycee_WGS84) sp_lycee_WGS84@data$TxCadres = extract(R.TxCadres,sp_lycee_WGS84) sp_lycee_WGS84@data$TxOuvriers = extract(R.TxOuvriers,sp_lycee_WGS84) sp_lycee_WGS84@data$NbPopulation = extract(R.NbPopulation,sp_lycee_WGS84) sp_lycee_WGS84@data$TxSenior = extract(R.TxSenior,sp_lycee_WGS84) sp_lycee_WGS84@data$RevenusMedians = extract(R.RevenusMedians,sp_lycee_WGS84) |
On peut enfin conclure, en faisant une fonction generique de visualisation
Creation_Graphique = function(df, AnneeObs_Ouv, AnneeObs_Clo, Effectifs, Abscisses, Ordonnees, TitreAbs, TitreOrd, CouleurGraph, CouleurLiss, Serie) { df_temp = df[which(df$Annee>=AnneeObs_Ouv & df$Annee<=AnneeObs_Clo),] df_temp = df_temp[which(!is.na(df_temp[,Effectifs]) & !is.na(df_temp[,Abscisses]) & !is.na(df_temp[,Ordonnees])),] df_temp = df_temp[!df_temp[,Effectifs]==0,] df_temp = df_temp[,c(Abscisses,Ordonnees)] graphique = ggplot(df_temp,aes(x = df_temp[,Abscisses],y = df_temp[,Ordonnees])) graphique = graphique + geom_point(data = df_temp, aes(x = df_temp[,Abscisses],y = df_temp[,Ordonnees]),size=1, color=CouleurGraph,alpha=0.25) graphique = graphique + geom_density2d(aes(colour=..level..),show.legend=F) + scale_colour_gradient(low="gray55",high="gray25") graphique = graphique + scale_y_continuous(breaks= seq(80,100,by=2), limits = c(80,100)) graphique = graphique + xlab(TitreAbs) + ylab(TitreOrd) graphique = graphique + ggtitle(Serie) graphique = graphique + theme(plot.title = element_text(size=13,color=CouleurLiss, face="bold", hjust=0), axis.title.x = element_text(size=8, face="bold"), axis.text.x = element_text(size=8), axis.title.y = element_text(size=8, face="bold"), axis.text.y = element_text(size=8), panel.border = element_rect(linetype = "solid"), panel.grid.major = element_line(colour = "gray55",linetype = "dashed"), panel.grid.minor = element_line(colour = "gray75",linetype = "dashed")) graphique = graphique + stat_smooth(method = "loess",fill=CouleurLiss,color=CouleurLiss) return(graphique) } Production_Graphique_VI_1 = function(df, Titre_General, Axe_Abscisses, Titre_Abscisses, Annee_Observee_Ouv, Annee_Observee_Clo){ Graph_S = Creation_Graphique(df, Annee_Observee_Ouv, Annee_Observee_Clo, "Effectif.S", Axe_Abscisses, "Tx.Admis.Etab.S", Titre_Abscisses, "Taux d'admission (%)", "dodgerblue3","dodgerblue4","Série S") Graph_ES = Creation_Graphique(df, Annee_Observee_Ouv, Annee_Observee_Clo, "Effectif.ES", Axe_Abscisses, "Tx.Admis.Etab.ES", Titre_Abscisses, "Taux d'admission (%)","darkorange2","darkorange3","Série ES") Graph_L = Creation_Graphique(df, Annee_Observee_Ouv, Annee_Observee_Clo, "Effectif.L", Axe_Abscisses, "Tx.Admis.Etab.L", Titre_Abscisses, "Taux d'admission (%)", "chartreuse4","darkgreen","Série L") Graph_TS = Creation_Graphique(df, Annee_Observee_Ouv, Annee_Observee_Clo, "Effectif.Etab", Axe_Abscisses, "Tx.Admis.Etab", Titre_Abscisses, "Taux d'admission (%)", "indianred1","red4","Toutes séries") p = plot_grid(Graph_S, Graph_ES, Graph_L, Graph_TS, ncol = 2, nrow = 2,align = 'hv', scale = c(0.95, 0.95, 0.95, 0.95),vjust = 0.9, hjust=-0.5) titre <- ggdraw() + draw_label(Titre_General,fontface="bold", size=10) plot_grid(titre, p, ncol = 1, rel_heights=c(.25,5)) } |
On note ici
df_lycee <- sp_lycee_WGS84@data |
et on peut faire un premier graphique, avec le taux de chomage
Production_Graphique_VI_1(df = df_lycee, Titre_General = "Taux d'admission par série en fonction du taux de chômage \n dans la population active - Tous lycées confondus", Axe_Abscisses = "TxChomage", Titre_Abscisses = "Taux de chômage dans la population active (%)", Annee_Observee_Ouv = "2013", Annee_Observee_Clo = "2015") |
On ne va pas enfoncer les portes ouvertes de l’inference ecologique, en affirmant des choses aussi stupides que “on a moins de chances d’avoir le bac quand on est au chômage”. Mais on peut noter que dans les zones avec un fort taux de chômage, les résultats au bac sont moins bons.
On peut ensuite regarder en fonction du revenu de la commune du lycée
Production_Graphique_VI_1(df = df_lycee, Titre_General = "Taux d'admission par série en fonction du niveau des revenus disponibles médians - Tous lycées confondus", Axe_Abscisses = "RevenusMedians", Titre_Abscisses = "Quantile du niveau des revenus disponibles médians (%)", Annee_Observee_Ouv = "2013", Annee_Observee_Clo = "2015") |
Fascinant, non ? Mais c’est clairement juste une première approche… il faudrait aller plus loin ensuite !
Scraper, ou pas ?
Ce matin, je mettais en ligne un billet scraper la base d’incendies de forêts expliquant comment scraper la base Promethee, en remplaçant automatique une formulaire. Dans la soirée, Alek Vladimir Nevski sur Twitter me faisait remarquer que c’était intéressant, mais peut être inutile (sur cet exemple particulier en tous cas).
Sur Firefox, il existe un “Moniteur Réseau“, qui s’ouvre en cliquant sur le clavier Ctrl + Maj + E ( Command + Option + E quand on est est sur un Mac). Ça ouvre un espace en bas de la fenêtre
Si maintenant je lance a la main ma requête, on voit ce qui est fait
Une methode GET apparait. En cliquant dessus, on a toute l’information nécessaire a droite
Dans l’onglet ‘Reponse’, on voit l’information sur le fichier json créé,
on a plus d’information bien entendu, si on agrandit
et surtout si on va dans l’onglet ‘En-tetes’
on a l’URL complet de la requête ! On peut ensuite copier le lien, et l’ouvrir
On a ainsi a accès a presque toute la table. Quand on regarde l’URL, seuls 20 incendies sont renvoyés, a cause du &nbrLigne=20& mais on peut tenter de le changer. Alek Vladimir Nevski propose de mettre &nbrLigne=10000& et il affirmait que ça marchait. J’avoue avoir essayé, sans succès.
Mais je retiens le tuyau en tous cas. Parfois, on peut faire simple !
Scraper la base d’incendies de forêts
Aller, deux billets pour parler de scraping, ou harvesting serait peut être plus juste ? Disons qu’on va faire un code automatique pour récupérer des informations en ligne. Plus particulièrement, on va s’intéresser au cas ou les données sont obtenues par une requête manuelle, et le site n’a pas d’API pour générer un joli fichier csv. Pierre-Marie, de la formation Data Science pour l’Actuariat a propose un exemple dont je vais m’inspirer (avec une autre methode pour scraper car j’ai eu des soucis pour installer le logiciel).
En 1973 a été lancée en France Prométhée, une base de données sur les incendies de forêts de la région méditerranéenne (couvrant 15 départements du Sud-Est). Oui, Προμηθεύς est ce titan qui a voulu jouer avec le feu… On peut consulter la base sur http://promethee.com/default/incendies,
En remplissant les champs, on peut faire une extraction suivant tel ou tel critère. Par exemple, en 2016, en Corse, il y a eu 6 incendies de nature “accidentelle”,
On le voit, le lien de la page ne change pas… il faut donc trouver une methode pour remplir manuellement les champs, et récupérer les informations. En particulier, on voit qu’on peut avoir accès aux informations relatives aux différents sinistres en cliquant sur la loupe
Idéalement, on voudrait un code R, automatique, pour faire ça ! Une première solutions, proposée par Pierre-Marie, est de passer par RSelenium. Le package R (non maintenu sur le CRAN depuis quelques semaines) est base sur Selenium, “Selenium automates browsers” comme dit la publicité. Le code pour scraper une page comme par
require("devtools") devtools::install_github("ropensci/RSelenium") require(RSelenium) remDr = remoteDriver(remoteServerAddr = "localhost" , port = 4445L , browserName = "firefox") |
Le soucis est que ce code ne tourne pas sur mon mac. Il y a de nombreuses pages qui parlent des spécificités de mac (en particulier le soucis du port par défaut 4445, qu’il faut changer). J’ai tout essayé, sans succès.
Une autre piste est de passer par un autre package
require("devtools") devtools::install_github("ropensci/wdman") |
wdman, pour “Webdriver Manager“, propose aussi de le faire, en s’appuyant sur diverses technologies, selenium, chromedriver, phantomJS binary ou encore geckodriver. Moyennant l’installation de chrome, la seconde option est intéressante
library(wdman) cDrv = chrome() eCaps = list(chromeOptions = list( args = c('--headless', '--disable-gpu', '--window-size=1280,800') )) remDr = remoteDriver(browserName = "chrome", port = 4567, extraCapabilities = eCaps) remDr$open() |
qui marche ! On peut lui demander d’aller sur la base promethee, et de faire une copie d’écran
remDr$navigate("http://www.promethee.com/default/incendies") remDr$screenshot(display = TRUE) |
On va ensuite remplir les champs, pour notre requête, sur la page. Par exemple, pour la date de début, on peut remonter au début des années 80
webElem = remDr$findElement(using = 'name', 'dtAlerteDeb') webElem$clearElement() webElem$sendKeysToElement(list('01/01/1981')) |
et terminer fin 2017
webElem = remDr$findElement(using = 'name', 'dtAlerteFin') webElem$clearElement() webElem$sendKeysToElement(list('31/12/2017')) |
On se focalise seulement sur la Corse
webElem = remDr$findElement(using = 'name', 'codeg') webElem$sendKeysToElement(list('CORSE')) |
et on veut 50 observations par page
webElem = remDr$findElement(using = 'name', 'nbrLigne') webElem$sendKeysToElement(list('50')) |
Ah oui, et on va temporiser entre chaque requête
webElem = remDr$findElement(using = 'name', 'btnSubmit') webElem$submitElement() Sys.sleep(15) remDr$screenshot(TRUE) |
(c’est une règle assez classique quand on scrap: on évite d’aller trop vite, pour pas se faire repérer…). On va ensuite commencer a regarder ce que l’on demande
text_base = read_html(remDr$getPageSource()[[1]]) list_tr = html_nodes(text_base, "tr") a = html_nodes(list_tr[length(list_tr)], "a") page_max = as.numeric(html_text(a[10])) |
On a environ 620 pages ! avec au moins 10 secondes pour changer de pages, et en rajoutant 15 secondes (on en parlait avant), ca fait quelques heures pour scraper.
Que veut-on récupérer en grattant la page ? On peut récupérer les coordonnées de l’incendie (latitude/longitude), les caractéristiques de l’incendie (dates, lieu, surface), et retourner le tout dans une table
scrapping = function(feux, page) { list_a = html_nodes(read_html(feux), "a") if (length(list_a) >= 2) { att = html_attrs(list_a[2]) longitude = as.numeric(str_extract(att, "[0-9]\\.([0-9])*")) latitude = as.numeric(str_extract(att, "([0-9]){2}\\..([0-9])*")) } else { longitude = 0 latitude = 0 } list_span = html_nodes(read_html(feux), "span") date = html_text(list_span[2]) commune = html_text(list_span[4]) surface = as.numeric(html_text(list_span[5])) list_div = html_nodes(read_html(feux), "div") if (length(list_a) >= 2) { infos1 = html_text(list_div[3]) infos2 = html_text(list_div[5]) } else { infos1 = html_text(list_div[2]) infos2 = html_text(list_div[4]) } alerte = gsub("First alert : |First intervention", "", str_extract(infos1, "First alert : (.)*First intervention")) deb_inter = gsub("First intervention : |End of intervention", "", str_extract(infos1, "First intervention : (.)*End of intervention")) fin_inter = gsub("End of intervention : |Origins of the alert", "", str_extract(infos1, "End of intervention : (.)*Origins of the alert")) code_INSEE = gsub("INSEE code : |Locality", "", str_extract(infos2, "INSEE code : (.)*Locality")) lieudit = str_trim(gsub("Locality : ", "", str_extract(infos2, "Locality : (.)*"))) liste_feux = data.frame(commune, date, code_INSEE, longitude, latitude, lieudit, surface, alerte, deb_inter, fin_inter, page) return(liste_feux) } |
Reste a faire tourner sur les 620 pages. Pour chaque ligne, on récupère les informations !
df = data.frame() for (page in 1:page_max) { text_base = read_html(remDr$getPageSource()[[1]]) list_tr = html_nodes(text_base, "tr") if (length(html_nodes(read_html(as.character(list_tr[14])), "a")) == 8) { deb = 15 } else { deb = 14 } for (i in deb:(deb+49)) { feux = as.character(list_tr[i]) df = rbind(df, scrapping(feux, page)) } if (page == 1) { webElem = remDr$findElement(using = 'xpath', '//*[contains(concat( " ", @class, " " ), concat( " ", "page_img", " " ))]') } else { webElem = remDr$findElement(using = 'xpath', '//*[contains(concat( " ", @class, " " ), concat( " ", "page_img", " " ))][3]') } webElem$clickElement() heure = (floor((page_max-page)*15/60) - floor((page_max-page)*15/60) %% 60) / 60 minute = floor((page_max-page)*15/60) %% 60 pb = progress_bar$new(total = 100) for (i in 1:100) { if (i==1) {cat("page", page, "sur", page_max," - (reste", heure, "h", minute, "min)", "\n")} pb$tick() Sys.sleep(15 / 100) } remDr$screenshot(TRUE) } feux = df |
Ça y est, on a une base, avec presque 31 000 lignes ! on va ensuite ranger un peu, et clarifier
colnames(feux)[colnames(feux)=="longitude"] = "x" colnames(feux)[colnames(feux)=="latitude"] = "y" feux[, date := as.Date(date, "%d/%m/%Y")] feux[, date_alerte := as.Date(substr(alerte, 1, 10), "%d/%m/%Y")] feux$heure_alerte[nchar(feux$alerte) > 10] = substr(feux$alerte[nchar(feux$alerte) > 10], 12, 16) feux$heure_alerte[!(nchar(feux$alerte) > 10)] = "" feux[, date_deb_inter := substr(deb_inter, 1, 10)] feux[, date_deb_inter := as.Date(date_deb_inter, "%d/%m/%Y")] feux$heure_deb_inter[nchar(feux$deb_inter) > 10] = substr(feux$deb_inter[nchar(feux$deb_inter) > 10], 12, 16) feux$heure_deb_inter[!(nchar(feux$deb_inter) > 10)] = "" feux$heure_deb_inter[feux$heure_deb_inter == "00:00"] = "" feux[, date_fin_inter := substr(fin_inter, 1, 10)] feux[, date_fin_inter := as.Date(date_fin_inter, "%d/%m/%Y")] feux$heure_fin_inter[nchar(feux$fin_inter) > 10] = substr(feux$fin_inter[nchar(feux$fin_inter) > 10], 12, 16) feux$heure_fin_inter[!(nchar(feux$fin_inter) > 10)] = "" feux$heure_fin_inter[feux$heure_fin_inter == "00:00"] = "" feux$annee = year(feux$date_alerte) feux$month = month(feux$date_alerte) |
Pour voir un peu nos données, rien de tel que des statistiques descriptives. On peut commencer avec du barplot du nombre de feux de forêts par année
agg_annee = aggregate(alerte~annee, FUN=length, data=feux) colnames(agg_annee)[colnames(agg_annee)=="alerte"] = "Nombre" nb = ggplot(agg_annee) + geom_bar(aes(x=annee, y=Nombre, fill="Feux de forêts"), stat="identity", width = 0.8) + scale_colour_manual(name=NULL, values=c("orange")) + ggtitle("Nombre de feux par an") + theme(plot.title = element_text(hjust = 0.5, size=14), legend.position = "none", axis.text.x = element_text(angle=50, vjust=0.5, size=10), axis.text.y = element_text(size=12), axis.title.x = element_blank(), axis.title.y = element_blank()) + scale_fill_discrete("") + scale_x_discrete(limit=seq(1981, 2017, by = 2)) |
ou un barplot de la surface de feux de forêts par anee
agg_surface = aggregate(surface~annee, FUN=sum, data=feux) surf = ggplot(agg_surface) + geom_bar(aes(x=annee, y=surface, fill="Feux de forêts"), stat="identity", width = 0.8) + scale_colour_manual(name=NULL, values=c("orange")) + ggtitle("Surface totale incendiée par an (en ha)") + theme(plot.title = element_text(hjust = 0.5, size=14), legend.position = "none", axis.text.x = element_text(angle=50, vjust=0.5, size=10), axis.text.y = element_text(size=12), axis.title.x = element_blank(), axis.title.y = element_blank()) + scale_fill_discrete("") + scale_x_discrete(limit=seq(1981, 2017, by = 2)) + scale_y_discrete(limit=seq(0, 25000, by = 5000)) |
Et on fait le graphique
require(grid) grid.newpage() pushViewport(viewport(layout = grid.layout(1, 2))) define_region = function(row, col){ viewport(layout.pos.row = row, layout.pos.col = col) } print(nb, vp = define_region(1, 1)) print(surf, vp = define_region(1, 2)) |
Passons aux choses serieuses, avec des cartes (la partie la plus serieuse, on en conviendra, c’etait la phase initiale de scraping). On commence par des informations relatives a la Corse
url="https://raw.githubusercontent.com/data-PM/Promethee/master/communes_corse-shp.zip" temp="corse.zip" download.file(url, "corse.zip") unzip("corse.zip") communes_corse = sp::spTransform(shapefile("/Users/arthur/communes_corse-shp/communes_corse.shp"), CRS("+proj=longlat +datum=WGS84")) |
On se limite aux feux de plus de 1000m2
feux_min1000m2 = feux[feux$surface > 0.1, ] |
Maintenant, on peut faire la carte
feux_INSEE = aggregate(x = feux_min1000m2$code_INSEE, by = list(unique.values = feux_min1000m2$code_INSEE), FUN = length) colnames(feux_INSEE)[colnames(feux_INSEE)=="unique.values"] = "insee" communes_corse = merge(x = communes_corse, y = feux_INSEE, by = "insee", all.x = TRUE) communes_corse$x[is.na(communes_corse$x)] = 0 max = max(communes_corse$x) communes_corse$x2 = communes_corse$x/max * 100 corse2 = gSimplify(communes_corse, tol = 0.001, topologyPreserve = TRUE) corse2 = SpatialPolygonsDataFrame(corse2, communes_corse@data, match.ID = FALSE) opar = par(mar = c(0,0,1.2,0)) plot(corse2, col = "#DAE3E6", border = "#8A0641", lwd = 0.7, bg = "#B5D0D0") |
Le fond de carte (avec les communes) ressemble a ca
On peut maintenant passer a la carte choroplethe
choroLayer(spdf = corse2, df = corse2@data, var = "x", border = NA, col = carto.pal("wine.pal", 6), legend.pos = "topleft", add = TRUE, method = "fisher-jenks", nclass = 6, legend.title.txt = NULL) top_villes = corse2[corse2$nom %in% c("Ajaccio", "Bastia", "Corte", "Porto-Vecchio", "Propriano", "Bonifacio", "Calvi"), ] graphics::points(coordinates(top_villes), pch = 20, cex = 0.5) labelLayer(top_villes, top_villes@data, txt = "nom", cex = 1, pos = 2, font = 4, offset = 0.5) layoutLayer(title = "Nombre de feux de forêts depuis 1981 (minimum 1000 m2)", source = "Source : 1973-2018 Prométhée ©", author = NULL, scale = 10, frame = FALSE, col = "#cdd2d4", coltitle = "#8A5543") |
Comme souvent avec R, plusieurs formats sont envisageables. On peut passer a du ggplot2, mais on commence par aller chercher un shapefile global de la Corse
corse = shp_zip("https://raw.githubusercontent.com/data-PM/Promethee/master/corse-shp.zip", "corse") corse = spTransform(corse, CRS("+proj=longlat +datum=WGS84")) corse = fortify(corse) p=data %>% ggplot() + geom_polygon(data=corse, aes(x=long, y = lat, group = group), fill="grey", alpha=0.3) + geom_point(data=feux_min100ha, aes(x=x, y=y, size=surface, color=surface, text=mytext, alpha=surface) ) + labs(title = "\nPrincipaux feux de forêts depuis 1981 \n (plus de 100 Ha)") + xlim(8.5, 9.6) + ylim(41.3,43.1) + scale_size_continuous(range=c(1,5), ) + scale_color_viridis(option="inferno", trans="log", name="Surface", limits=c(99, 5700), breaks=c(100, 200, 500, 1000, 2500, 5000)) + scale_alpha_continuous(trans="log") + theme_void() + theme(legend.title = element_text(size=10, face="bold")) t=800 div(ggplotly(tooltip="text", width=t/1.7, height=t), align = "center") |
Joli, non ?
Carte météo dynamique
Pour aller un peu plus loin par rapport au précédent billet, on peut faire une carte dynamique, pour visualiser une tempête, sur une journée. Je vais partir du fait qu’on a encore les objets du précédant billet en mémoire.
On commence par garder uniquement les données relatives a la journée qui nous intéresse,
date_sel = "2016-10-13" date_sel_OK = paste0(substr(date_sel,9,10),"/",substr(date_sel,6,7),"/",substr(date_sel,1,4)) donnees_carte_0 = subset(data, data$date==as.POSIXct(date_sel)) horaires = sort(unique(donnees_carte_0$heure)) |
(mais on pourrait prendre une semaine). On recupere ensuite le contour de la France,
download.file("http://biogeo.ucdavis.edu/data/gadm2.8/rds/FRA_adm0.rds","FRA_adm0.rds") FR0=readRDS("FRA_adm0.rds") P1=FR0@polygons[[1]]@Polygons[[355]]@coords P2=FR0@polygons[[1]]@Polygons[[27]]@coords |
car seule la France métropolitaine nous intéressera (incluant la Corse, soit 2 gros polygones) dans l’esprit du précédant billet sur le zonier. On va ainsi definir un maillage sur lesquel on va lisser la pluviometrie,
grille = expand.grid(seq(min(donnees_carte_0$longitude),max(donnees_carte_0$longitude),length=101),seq(min(donnees_carte_0$latitude),max(donnees_carte_0$latitude),length=101)) paslong=(max(donnees_carte_0$longitude)-min(donnees_carte_0$longitude))/100 paslat=(max(donnees_carte_0$latitude)-min(donnees_carte_0$latitude))/100 f=function(i){ (point.in.polygon (grille[i, 1]+paslong/2 , grille[i, 2]+paslat/2 , P1[,1],P1[,2])>0)+(point.in.polygon (grille[i, 1]+paslong/2 , grille[i, 2]+paslat/2 , P2[,1],P2[,2])>0) } indic=unlist(lapply(1:nrow(grille),f)) grille=grille[which(indic==1),] |
Pour lisser, on utilise les k-plus proches voisins
knn=function(i,k=10){ d=distHaversine(grille[i,1:2],donnees_carte[, c("longitude","latitude")], r=6378.137) r=rank(d) ind=which(r<=k) weighted.mean(donnees_carte[ind,"rr3"],(1/d[ind])/sum(1/d[ind]))} |
On a ensuite la fonction suivante, pour faire une carte
carto_prec<-function(){ grille2<-grille grille2$rr=Vectorize(knn)(i=1:nrow(grille2)) bk=seq(0,50,length=21) grille2$cuty=cut(grille2$rr,breaks=bk,labels=1:20) cols = rev(carto.pal(pal1 = "blue.pal", n1=20, pal2 = "white.pal", n2=1)) plot(FR0,border=NA) polygon(P1) polygon(P2) points(grille2[,1]+paslong/2,grille2[,2]+paslat/2,col=cols[grille2$cuty],pch=19) points(donnees_carte$longitude,donnees_carte$latitude, col="black",pch=19,cex=.5) title(main = paste0(date_sel_OK," à ",heure_sel,"H"), sub = "Précipitations des 3 dernières heures", cex.main = 1.5, font.main= 4, cex.sub = 1, font.sub = 3) legend(8.2, 50, legend=seq(0,50,length=5), title='en mm', fill=cols[seq(1,20,length=5)], cex=0.8) } |
On va ensuite utiliser 8 fichiers – car on regarde toutes les 3 heures
for (hh in horaires) { heure_sel<-hh donnees_carte <- subset(donnees_carte_0, donnees_carte_0$heure==heure_sel) png(paste0("CartePrec",date_sel,"_",hh,".png")) carto_prec() dev.off() } |
(mais ça peut se raffiner, avec un lissage spatio-temporel). Et on va ensuite créer notre animation en concaténant les images
frames <- image_morph( c(image_scale(image_read("CartePrec2016-10-13_1.png")), image_scale(image_read("CartePrec2016-10-13_3.png")), image_scale(image_read("CartePrec2016-10-13_6.png")), image_scale(image_read("CartePrec2016-10-13_9.png")), image_scale(image_read("CartePrec2016-10-13_12.png")), image_scale(image_read("CartePrec2016-10-13_15.png")), image_scale(image_read("CartePrec2016-10-13_18.png")), image_scale(image_read("CartePrec2016-10-13_21.png"))), frames = 16) frames.anim <- image_animate(frames) image_write(frames.anim, paste0("CartePrecDay",date_sel,".gif")) |
Il faut avouer que ça a de gueule… On voit clairement la tempête arriver…
Faire ses cartes météo
Allez, un peu de scraping aujourd’hui (promis, il y aura d’autres billets d’ici la fin de la semaine sur le sujet). Je vais m’inspirer des codes de Romain de la formation Data Science pour l’Actuariat. Le but est de récupérer des données de températures, de précipitations, de vitesse de vent, journalières, en France. On commence par charger quelques librairies,
library(plyr) library(stringr) library(OpenStreetMap) library(leaflet) library(shiny) library(rsconnect) library(mapview) library(png) library(magick) library(yaml) |
On va récupérer les deux dernières années, mais comme toujours, le code peut facilement s’adapter.
annee = c(2016:2017) mois = c('01','02','03','04','05','06','07','08','09','10','11','12') aaaamm = sort(sub(pattern=' ',replacement = '', x=outer(annee,mois,paste))) myCols = c(NA,NA,"NULL","NULL","NULL",NA,NA,"factor","NULL","NULL","NULL",NA,"NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL",NA,NA,NA,NA,NA,"NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL") for (mois in aaaamm) { fichier = paste0("synop_",mois,".csv.gz") if (!file.exists(fichier)) { page = paste0('https://donneespubliques.meteofrance.fr/donnees_libres/Txt/Synop/Archive/synop.',mois,'.csv.gz') download.file(page, fichier, quiet = TRUE, cacheOK = FALSE) } } for (mois in aaaamm) { fichier = paste0("synop_",mois,".csv.gz") if (mois=='201601') { data = read.table(gzfile(fichier),header=TRUE, sep=";",quote='"',fileEncoding="UTF-8",colClasses = myCols) } else { data = rbind(data,read.table(gzfile(fichier),header=TRUE, sep=";",quote='"',fileEncoding="UTF-8",colClasses = myCols)) } } |
On supprime une des colonnes qui ne sert a rien
data[,"X"] = NULL |
et on bricole un peu (sinon on a des soucis avec minuit, qui est décalé de 24 heures)
data$temp = as.numeric(substr(data$date,1,14)) decalage_minuit = function(dateheure){ r = dateheure if (dateheure-trunc(dateheure/1000000)*1000000==0){ r = dateheure + 10000 } return(r) } data$temp = sapply(data$temp,decalage_minuit) data$jourheure_obs = strptime(data$temp, format = '%Y%m%d%H%M%S', 'UTC') data$date = as.POSIXct((strptime(data$jourheure_obs,format = "%Y-%m-%d"))) data$heure = as.numeric(substr(data$jourheure_obs, start=12, stop=13)) data[,"temp"] = NULL data$t = as.numeric(as.character(data$t))-273.15 |
oui, on doit convertir en degrés Celsius. On peut aussi travailler sur les vitesses de vent
data$dd = as.numeric(as.character(data$dd)) data$ff = as.numeric(as.character(data$ff))/1000*3600 data$ff = round(data$ff/5)*5 |
que l’on va convertir en km/h, et enfin les précipitations
data$rr1 = pmax(0,as.numeric(as.character(data$rr1))) data$rr3 = pmax(0,as.numeric(as.character(data$rr3))) data$rr6 = pmax(0,as.numeric(as.character(data$rr6))) data$rr12 = pmax(0,as.numeric(as.character(data$rr12))) data$rr24 = pmax(0,as.numeric(as.character(data$rr24))) |
On peut d’ailleurs commencer par une simple visualisation de séries de températures
data_temperature = subset(data, !is.na(data$t)) temp_min = aggregate(x=data_temperature$t,by=list(date=data_temperature$date),FUN=min) names(temp_min)[2] = "temp_min" temp_max = aggregate(x=data_temperature$t,by=list(date=data_temperature$date),FUN=max) names(temp_max)[2] = "temp_max" temp_moy = aggregate(x=data_temperature$t,by=list(date=data_temperature$date),FUN=mean) names(temp_moy)[2] = "temp_moy" temp_stat = cbind(temp_min,temp_max[],temp_moy) rm(temp_min,temp_max,temp_moy) temp_stat[,3] = NULL temp_stat[,4] = NULL |
avec le maximum observé, par jour, sur toutes les stations, le minimum, et la moyenne des stations
plot(x=temp_stat$date,y=temp_stat$temp_moy,'l',ylim=c(-15,40),xlab='Date',ylab='Température (°C)') title('Températures françaises métropolitaines de 01/2016 à 12/2017') points(x=temp_stat$date,y=temp_stat$temp_min,col="blue",'l') points(x=temp_stat$date,y=temp_stat$temp_max,col="red",'l') |
Pour faire des jolies cartes, il nous faut des informations sur les stations météo.
stations = read.table('https://donneespubliques.meteofrance.fr/donnees_libres/Txt/Synop/postesSynop.csv',header=TRUE, sep=";",quote='"',fileEncoding="UTF-8") stations = rename(stations, replace=c("ID"="numer_sta")) stations = rename(stations, replace=c("Nom"="nom")) stations = rename(stations, replace=c("Latitude"="latitude")) stations = rename(stations, replace=c("Longitude"="longitude")) stations = rename(stations, replace=c("Altitude"="altitude")) data = merge(x=data,y=stations) |
Pour la carte, on va se restreindre aux données en France métropolitaine
data = subset(data, data$latitude >= 39) data = subset(data, data$latitude <= 53) data = subset(data, data$longitude >= -5.3) data = subset(data, data$longitude <= 9.8) coordonnees_stations = unique(data[, c("nom","latitude","longitude")]) |
Commencons par visualiser le maximum observee le 19 juillet 2016
date_sel = "2016-07-19" date_sel_OK = paste0(substr(date_sel,9,10),"/",substr(date_sel,6,7),"/",substr(date_sel,1,4)) donnees_carte_0 = subset(data, data$date==as.POSIXct(date_sel)) temp = aggregate(x=donnees_carte_0$t,by=list(numer_sta=donnees_carte_0$numer_sta),FUN="max") donnees_carte_0 = merge(x=donnees_carte_0,y=temp) donnees_carte_0 = rename(donnees_carte_0, replace=c("x"="TEMP")) donnees_carte_0$TEMP = round(donnees_carte_0$TEMP) donnees_carte = unique(donnees_carte_0[,c('numer_sta','nom','longitude','latitude','TEMP')]) get_france = get_map(c(lon=2.25,lat=46), zoom=5, col='bw') colfunc = colorRampPalette(c("darkblue","royalblue","cyan","lightblue","orange","red","darkred")) liste_couleurs = data.frame(TEMP=as.numeric(seq(-49,50,1)),couleur=as.character(colfunc(100))) donnees_carte = merge(x=donnees_carte,y=liste_couleurs) france_tmax = ggmap(get_france) + ggtitle(paste0('Températures maximales du ', date_sel_OK)) + scale_x_continuous(limits = c(-5, 10), expand = c(0, 0)) + scale_y_continuous(limits = c(41.3, 51.1), expand = c(0, 0)) france_tmax = france_tmax + geom_label(data=donnees_carte,aes(x=longitude,y=latitude,label=TEMP),col=as.character(donnees_carte$couleur)) france_tmax |
On peut tenter aussi une visualisation des précipitations. On se place le 13 octobre 2016 (oui, il y a eu des inondations dans la region de Montpelier ce jour la)
date_sel = "2016-10-13" colfunc = colorRampPalette(c('white','darkblue')) data$precip_en_mm = data$rr3 data_precipitation = subset(data, !is.na(data$precip_en_mm)) precipitations = aggregate(x=data_precipitation$precip_en_mm,by=list(date=data_precipitation$date,station=data_precipitation$nom),FUN=sum) names(precipitations)[3] = "pp" precipitations$pp = pmax(precipitations$pp,0) donnees_carte = subset(precipitations, precipitations$date==date_sel) qt_prec = quantile(round(donnees_carte$pp),seq(0, 1, 0.05)) donnees_carte$Couleur_prec = colfunc(21)[(findInterval(round(donnees_carte$pp), qt_prec, all.inside=TRUE))] donnees_carte = rename(donnees_carte, replace=c("station"="nom")) donnees_carte = merge(x=donnees_carte,y=coordonnees_stations) |
On peut tenter un leaflet ici (ça sera plus joli)
france = leaflet() %>% addTiles() %>% fitBounds(lng1=-4.412, lat1=41.92, lng2=9.485, lat2=50.57) france = france %>% addCircles(lng=donnees_carte$longitude, lat=donnees_carte$latitude, color=donnees_carte$Couleur_prec, opacity = 1, fillColor=donnees_carte$Couleur_prec, fillOpacity = 1, radius=donnees_carte$pp*500) france |
Comme souvent, j’ai des soucis pour integrer du leaflet en wordpress, alors on se contentera dans le billet d’une copie d’ecran (mais un copier/coller du code suffit).
Et pour finir, la meme chose avec les vitesses de vent. On se place le 6 mars 2017 (la encore, la date n’est pas choisie au hasard, c’était la tempete Zeus)
date_sel = "2017-03-06" date_sel_OK = paste0(substr(date_sel,9,10),"/",substr(date_sel,6,7),"/",substr(date_sel,1,4)) heure_sel = 18 donnees_carte_0 = subset(data, data$date==as.POSIXct(date_sel)) donnees_carte_0 = subset(donnees_carte_0, donnees_carte_0$heure==heure_sel) donnees_carte = donnees_carte_0 |
Pour les precipitations, des cercles plus ou moins grands suffisaient. La, on va mettre des fleches
logo = image_trim(image_read("https://image.freepik.com/icones-gratuites/fleche-noire-vers-le-haut_318-30934.jpg")) if (file.exists('flecheblanche90.png')==FALSE){ for (angle in seq(10,360, by=10)) { temp = image_rotate(logo,angle) image_write(temp,paste0('flecheblanche',angle,'.png')) temp = image_colorize(image=temp,opacity=40,color='green') image_write(temp,paste0('flecheverte',angle,'.png')) temp = image_colorize(image=temp,opacity=40,color='yellow') image_write(temp,paste0('flechejaune',angle,'.png')) temp = image_colorize(image=temp,opacity=40,color='orange') image_write(temp,paste0('flecheorange',angle,'.png')) temp = image_colorize(image=temp,opacity=40,color='red') image_write(temp,paste0('flecherouge',angle,'.png')) temp = image_colorize(image=temp,opacity=40,color='black') image_write(temp,paste0('flechenoire',angle,'.png')) } } arrowIcons = icons( iconUrl = ifelse(donnees_carte$dd == 0, ifelse(donnees_carte$ff<39, 'flecheverte360.png', ifelse(donnees_carte$ff<79, 'flechejaune360.png', ifelse(donnees_carte$ff<119, 'flecheorange360.png', ifelse(donnees_carte$ff<159, 'flecherouge360.png', ifelse(donnees_carte$ff>159, 'flechenoire360.png', 'flecheblanche360.png'))))), ifelse(donnees_carte$ff<39, paste0('flecheverte',donnees_carte$dd,'.png'), ifelse(donnees_carte$ff<79, paste0('flechejaune',donnees_carte$dd,'.png'), ifelse(donnees_carte$ff<119, paste0('flecheorange',donnees_carte$dd,'.png'), ifelse(donnees_carte$ff<159, paste0('flecherouge',donnees_carte$dd,'.png'), ifelse(donnees_carte$ff>159, paste0('flechenoire',donnees_carte$dd,'.png'), paste0('flecheblanche',donnees_carte$dd,'.png'))))))), iconWidth = 30, iconHeight = 30) |
Cette fois, on est bon ! On peut faire la carte
france = leaflet(options = leafletOptions()) %>% addTiles() %>% fitBounds(lng1=-4.412, lat1=41.92, lng2=9.485, lat2=50.57) france = addMarkers(map = france,lng=donnees_carte$longitude,lat=donnees_carte$latitude,icon=arrowIcons,popup=as.character(donnees_carte$ff)) france |
Visualiser la localisation des dinosaures, grâce aux quaternions…
Allez, un petit billet que je rêve de faire depuis des années, mais je viens enfin de trouver un prétexte, grâce à François de la formation Data Science pour l’Actuariat (à qui je vais emprunter l’application qui va suivre). Il y a plusieurs années, j’avais évoqué dans un billet les difficultés de travailler sur les données spatiales, surtout quand on est au pole. Le pole est un point de singularité dans la représentation classique (latitude, longitude) de la sphère. A l’époque, j’avais propose un bricolage simpliste (mettre les points de singularité dans un coin), mais il y a plus propre, à l’aide des quaternions d’Hamilton. Pour les présenter rapidement, et le lien avec les données spatiales, sur le globe terrestre, je vais reprendre la page wikipedia qui explique de manière incroyablement claire l’idée générale, et plus particulièrement les liens des quaternions avec les rotations (qui servent de base dans la représentation de la sphère qu’est le globe terrestre).
Chaque rotation en dimension trois consiste à tourner d’un certain angle \alpha autour d’un certain axe \vec{v}. Pour un angle petit mais non nul (s’il est nul, on parle de rotation identité), l’ensemble des rotations possibles est une petite sphère entourant la rotation identité, où chaque point de la sphère représente un axe pointant dans une direction particulière. Des rotations d’angles de plus en plus grands s’éloignent progressivement de la rotation identité. Aussi, au voisinage de la rotation identité, l’espace abstrait des rotations ressemble à l’espace ordinaire en trois dimensions (qui peut également être vu comme un point central entouré de sphères de différents rayons).
On peut assimiler les différentes directions à partir du pôle (c’est-à-dire les différents méridiens) aux différents axes de rotations et les différentes distances au pôle Nord aux différents angles : on a ainsi une analogie de l’espace des rotations. Mais la surface de la sphère est en deux dimensions alors que les axes de rotation utilisent déjà trois dimensions. L’espace des rotations est donc modélisé par une sphère de dimension 3 dans un espace à 4 dimensions (une hypersphère). Il faut alors penser la sphère ordinaire comme à une section de l’hypersphère (de la même façon qu’un cercle est une section de sphère). On peut prendre la section pour représenter, par exemple, uniquement les rotations d’axes dans le plan xy. Et on peut légitimement penser maintenant aux rotations comme à des points de la sphère en dimension 4.
Bon, rentrons dans les détails. On paramètre la surface d’une sphère à l’aide de deux coordonnées, comme la latitude et la longitude. Mais cela pose des soucis importants aux pôles (comme je l’avais note dans un ancien billet). Le théorème de la boule chevelue montre en fait qu’il n’existe aucun système de coordonnées à deux paramètres évitant cette dégénérescence. On va donc plonger la sphère dans l’espace à trois dimensions, en la paramétrant a l’aide de trois coordonnées cartésiennes (ici w, x et y). Par convention, on place le pôle Nord à (w,y,z) = (1, 0, 0), le pôle Sud à (w,y,z) = (-1, 0, 0) et l’équateur sera le cercle d’équations w = 0 et x^2+y^2=1. Un point (w,x,y) de la sphère représente une rotation de l’espace ordinaire autour de l’axe horizontal dirigé par le vecteur {\displaystyle {\vec {v}}={\begin{pmatrix}x\\y\\0\end{pmatrix}}} et d’angle {\displaystyle \alpha =2\cos ^{-1}w=2\sin ^{-1}{\sqrt {x^{2}+y^{2}}}}
C’est l’idée générale.
Pour parler un peu des quaternions, pour rappel, un plan en dimension 2 peut être paramétré en utilisant les nombres complexes, en introduisant un symbole abstrait \mathbf {i} qui vérifie la règle \mathbf {i}^2=-1. On peut faire la meme chose en dimension 4, en introduisant des symboles abstraits \mathbf {i}, \mathbf {j} et \mathbf {k}. La partie imaginaire {\displaystyle b\mathbf {i} +c\mathbf {j} +d\mathbf {k} } d’un quaternion se comporte comme un vecteur {\displaystyle {\vec {v}}={\begin{pmatrix}b\\c\\d\end{pmatrix}}} d’un espace vectoriel à trois dimensions.
Définissons le quaternion {\displaystyle \mathbf {q}=w+x\mathbf {i} +y\mathbf {j} +z\mathbf {k} =\cos(\alpha /2)+\frac{\vec {v}}{\|v\|}\sin(\alpha /2)} ou bien {\displaystyle \mathbf {q}=w+x\mathbf {i} +y\mathbf {j} +z\mathbf {k} =\cos(\alpha /2)+{\vec {u}}\sin(\alpha /2)} où {\displaystyle {\vec {u}}} est un vecteur unitaire. Soit également {\displaystyle {\vec {v}}} un vecteur ordinaire de l’espace en 3 dimensions, considéré comme un quaternion avec une coordonnée réelle nulle. On pourrait que le produit de quaternions{\displaystyle \mathbf{q}{\vec {\nu}}\mathbf{q}^{-1}} renverrait le vecteur {\displaystyle {\vec {\nu}}} tourné d’un angle {\displaystyle \alpha } autour de l’axe dirigé par {\displaystyle {\vec {u}}} . Et c’est effectivement ce qui se passe. Cette opération est connue comme la conjugaison par \mathbf {q}.
Aussi, la multiplication de quaternions correspond à la composition de rotations, car si \mathbf {p} et \mathbf {q} sont des quaternions représentant des rotations, alors la rotation (conjugaison) par \mathbf {pq} est{\displaystyle \mathbf{pq}q{\vec {\nu}}(\mathbf{pq})^{-1}=\mathbf{pq}{\vec {\nu}}\mathbf{q}^{-1}\mathbf{p}^{-1}=\mathbf{p}(\mathbf{q}{\vec {\nu}}\mathbf{q}^{-1})\mathbf{p}^{-1}} ce qui revient à tourner (conjuguer) par \mathbf {q}, puis par \mathbf {p}.
Le quaternion inverse d’une rotation correspond à la rotation inverse, car {\displaystyle q^{-1}(q{\vec {v}}q^{-1})q={\vec {v}}} . Et assez naturellement, le carré d’un quaternion – noté \mathbf {q}^2 – correspond à la rotation de deux fois le même angle autour du même axe. Plus généralement, \mathbf {q}^{ {n}} correspond à une rotation de {n} fois l’angle autour du même axe que \mathbf {q}. Par convention, on peut considérer un réel arbitraire {r}, ce qui permet de calculer des rotations intermédiaires de façon fluide entre des rotations de l’espace.
Un petit exemple. Considérons la rotation f autour de l’axe dirigé par {\displaystyle {\vec {v}}=\mathbf {i} +\mathbf {j} +\mathbf {k} } et d’angle 120°, soit 2\pi/3.
La norme de {\displaystyle {\vec {v}}} est \sqrt{3} , le demi-angle est \pi/3 (ou 60°), le cosinus de ce demi-angle est 1/2, et le sinus est \sqrt{3}/2. Nous devons donc conjuguer avec le quaternion unitaire {\displaystyle \mathbf{q}=\cos {\frac {\pi }{3}}+\sin {\frac {\pi }{3}}\cdot {\frac {1}{\sqrt {3}}}{\vec {v}}} soit \mathbf{q}={\frac {1}{2}}+{\frac {\sqrt {3}}{2}}\cdot {\frac {1}{\sqrt {3}}}{\vec {v}}={\frac {1}{2}}+{\frac {\sqrt {3}}{2}}\cdot {\frac {\mathbf {i} +\mathbf {j} +\mathbf {k} }{\sqrt {3}}} qui peut finalement s’écrire simplement{\frac {1+\mathbf {i} +\mathbf {j} +\mathbf {k} }{2}}
Si f est la fonction de rotation,{\displaystyle f(a\mathbf {i} +b\mathbf {j} +c\mathbf {k} )=\mathbf{q}(a\mathbf {i} +b\mathbf {j} +c\mathbf {k} )\mathbf{q}^{-1}}
On peut prouver que l’on obtient l’inverse d’un quaternion unitaire simplement en changeant le signe de ses coordonnées imaginaires. Autrement dit{\displaystyle \mathbf{q}^{-1}={\frac {1-\mathbf {i} -\mathbf {j} -\mathbf {k} }{2}}} et donc f(a\mathbf {i} +b\mathbf {j} +c\mathbf {k} ) s’ecrit {\displaystyle {\frac {1+\mathbf {i} +\mathbf {j} +\mathbf {k} }{2}}(a\mathbf {i} +b\mathbf {j} +c\mathbf {k} ){\frac {1-\mathbf {i} -\mathbf {j} -\mathbf {k} }{2}}} En appliquant les règles ordinaires de calcul avec les quaternions, on obtient {\displaystyle f(a\mathbf {i} +b\mathbf {j} +c\mathbf {k} )=c\mathbf {i} +a\mathbf {j} +b\mathbf {k} } Ah oui, et de la même manière qu’on peut associer une matrice 2\times2 à un nombre complexez=a+b\mathbf{i} ~\rightarrow~\begin{pmatrix}a&-b\\b&a\end{pmatrix} on peut associer une matrice 4\times4 à un quaternion\mathbf{q}=a+b\mathbf{i}+c\mathbf{j}+d\mathbf{k} ~\rightarrow~\begin{pmatrix}\quad a&\quad -b&\quad -c&\quad -d\\\quad b&\quad a&\quad -d&\quad c\\\quad c&\quad d&\quad a&\quad -b\\\quad d&\quad -c&\quad b&\quad a\end{pmatrix} (on peut aussi passer par unE matrice 2\times2 à coefficients complexes, mais ça ne servirait qu’à compliquer, ici).
On a vu qu’on pouvait associer une rotation à un quaternion, et un quaternion à une matrice. Si la rotation est d’axe \vec {OM}, où O est le centre de la terre, et M un point sur la surface terrestre, décrit par sa latitude et sa longitude, avec pour angle \alpha (comme décrit dans un billet sur stackoverflow), on a la fonction R suivante
quat = function(lat,long,ang=NA){ n = length(lat) lat = lat/180*pi long = long/180*pi x = cos(lat) * cos(long) y = cos(lat) * sin(long) z = sin(lat) if (is.na(ang)){ Q = matrix(c(x,y,z,rep(0,n)), ncol = 4) } else { Q = matrix(c(sin(ang/2*pi/180) * c(x,y,z), cos(ang/2*pi/180)), ncol =4) } return(Q) } |
avec la réciproque, permettant de passer d’un quaternion à une rotation, avec une description de l’axe comme auparavant (un point sur la sphère – sur le globe terrestre) et un angle (coordonnées polaires)
polaire = function(Q, digits=2) { Q = Q/norme(Q) ang = round(acos(Q[4])*2*180/pi,digits) n = norme(Q[1:3]) x = Q[1]/n y = Q[2]/n z = Q[3]/n lat = asin(z) * 180/pi if (z**2 == 1){ long = 0 } else { phi = (x+1i*y)/sqrt(1-z**2) long = Im(log(phi)) * 180/pi } c(lat,long,ang) } |
à condition de définir au préalable la norme du quaternion
norme = function(Q){ sqrt(sum(Q**2)) } |
On peut aussi définir le produits de quaternion (toutes les opérations sont décrites dans la page wikipedia), qui sera noté \otimes par la suite
pdt_quat = function(Q1,Q2){ Q=rep(0,4) Q[1:3] = Q1[4]*Q2[1:3]+Q2[4]*Q1[1:3]+pdt_vect(Q1[1:3],Q2[1:3]) Q[4] = Q1[4]*Q2[4] - pdt_scal(Q1[1:3],Q2[1:3]) return(Q) } |
mais aussi un produit scalaire
pdt_scal = function(M,N){ return(sum(M*N)) } |
un produit vectoriel
pdt_vect = function(M,N){ return(c(M[2]*N[3]-M[3]*N[2], M[3]*N[1]-M[1]*N[3], M[1]*N[2]-M[2]*N[1])) } |
et finalement l’inverse du quaternion \mathbf{q}^{-1}
inv_quat = function(Q){ (c(0,0,0,2*Q[4])-Q)/norme(Q) } |
On peut aussi demander les coordonnées d’un point M' obtenu comme transformation d’un point M par la rotation \mathbf{q}
rotation = function(M, Q){ polaire(pdt_quat(pdt_quat(Q,M),inv_quat(Q)))[1:2] } |
Maintenant, on va pouvoir passer aux choses sérieuses….
Je l’ai évoqué en introduction, les quaternions peuvent permettre de contourner certains problèmes, comme manipuler des objets (comme la calotte glaciaire) qui sont situés autour du pole (qui est un point de singularité dans la représentation par coordonnées polaires). Une autre application, présentée par François, est celle du déplacement des plaques tectoniques. En particulier, on peut utiliser un fichier de rotations des plaques tectoniques entre aujourd’hui et une certaine date dans le passé (cette idée se retrouve dans le projet gplates programmé avec des librairies python). Ou plus généralement entre deux dates, t_1 et t_2. Le fichier a notre disposition contient ainsi des quaternions \mathbf{q}^{P_0}_{t,P} pour une date t et une plaque P (défini comme un polygone décrit par une collection de latitudes et de longitudes), ou le deplacement de la plaque est décrit relativement a la plaque P_0. Pour des soucis de calculs, on va supposer qu’on peut interpoler linéairement les quaternions,\mathbf{q}^{P_0}_{t,P}=(1-\lambda)\mathbf{q}^{P_0}_{t_1,P}+\lambda\mathbf{q}^{P_0}_{t_2,P}avec\lambda=\frac{t-t_1}{t_2-t_1}et que récursivement, on peut composer les rotations, au sens ou\mathbf{q}^{P_3}_{t,P_1}\mathbf{q}^{P_3}_{t,P_2}\otimes\mathbf{q}^{P_2}_{t,P_1}Aussi, a partir de notre fichier de rotations, on peut creer une fonction qui calcule l’ensemble des quaternions de rotation, pour les plaques plates données. A priori un quaternion \mathbf{q}^{P_0}_{t,P} ne sera calculé qu’une fois
projecteur = function(t,plates,rot){ ll = length(plates) Q0 = matrix(rep(0,4),ll,4) for (i in seq(ll)){ cur_plate = i Q = c(0,0,0,1) while(cur_plate > 0){ df = rot[rot$Start = t & rot$plate == plates[cur_plate],] if (dim(df)[1] > 0) { Q1 = Q0[cur_plate,] cur_plate = 0 if (norme(Q1)==0){ Q_St = quat(df$lat_St, df$lon_St, df$ang_St) Q_End = quat(df$lat_End, df$lon_End, df$ang_End) periode = df$End - df$Start pct = 0 if (periode > 0) pct = min(max(0,(t - df$Start)/periode),1) Q1 = Q_St + pct * (Q_End - Q_St) cur_plate = which(plates==df$anchor) if (length(cur_plate) == 0) cur_plate = 0 } Q = pdt_quat(Q1,Q) } else {cur_plate = 0} } Q0[i,] = Q } return(Q0) } |
On peut maintenant appliquer ces outils a des données. Ici, trois bases issues du site http://paleobiodb.org seront exploitées, pour visualiser ou les dinosaures vivaient :
– une base des collections recensant les sites de fouilles (et en particulier leur géolocalisation)
– une base d’occurrences recensant les spécimens trouvés, par collection
– une base des spécimens décrivant les spécimens de dinosaures.
site="http://paleobiodb.org/data1.2/" req=".txt?datainfo&rowcount&max_ma=999&min_ma=0" limit = "" #"&limit=100" names = c("colls/list", "occs/taxa", "occs/list") destfile = rep('',length(names)) for (i in 1:length(names)) { destfile[i] = paste0("data/",sub("/","_",names[i]),".csv") download.file(paste0(site,names[i],req,limit),destfile=destfile[i]) } |
Pour les données de collection
collection = read.csv(destfile[1],skip=17,sep = ",", header = TRUE) coll = select(collection, c(collection_no, lng, lat, max_ma, min_ma)) remove(collection) names(coll) = c('no', 'lng', 'lat', 'max_ma', 'min_ma') |
et pour les données d’occurrence
occurence = read.csv(destfile[3],skip=17, sep = ",", header = TRUE) occ = select(occurence, c(occurrence_no, collection_no, accepted_no)) remove(occurence) names(occ) = c('no', 'coll_no', 'taxo_no') |
Les espèces que nous étudierons ici sont les dinosaures de l’ordre des Ornithischia et des Saurischia (qui incluent les grandes familles classiques de dinosaures – a ce que j’ai pu comprendre)
taxonomie = read.csv(destfile[2],sep = ",", skip=20, header = TRUE) taxo = select(taxonomie,c(orig_no, accepted_rank, accepted_name, parent_no, container_no)) remove(taxonomie) names(taxo) = c('N0', 'rang0', 'nom0', 'parent', 'container') taxo = taxo[!is.na(taxo$N0),] taxo["N1"]=as.character(taxo$container) taxo[taxo$container=='',]$N1 = as.character(taxo[taxo$container=='',]$parent) taxo$container = as.factor(taxo$container) taxo$parent = NULL taxo$container = NULL t0 = taxo i = 0 while(dim(taxo[!is.na(taxo[paste0("N",i)]),])[1] > 0){ i = i+1 colnames(t0) = c(paste0('N',i),paste0('rang',i), paste0('nom',i), paste0('N',i+1)) taxo = merge(taxo,t0,all.x=TRUE) if (i==10) break } t1 = select(taxo,N0) for (niveau in levels(taxo$rang0)){ if (niveau != ""){ t1[niveau]="" for (j in seq(0,i)){ test = which(taxo[paste0("rang",j)]==niveau) t1[test,niveau] = as.character(taxo[test,paste0("nom",j)]) } } } taxo=t1[t1["unranked clade"%in%c("Ornithischia","Saurischia")],] head(taxo,5) |
On peut alors fusionner nos bases
M1 = merge(occ,taxo, by.x=c('taxo_no'),by.y=c('N0'), all=FALSE) M2 = merge(M1,coll,by.x=c('coll_no'),by.y=c('no'), all.x=TRUE) paste(dim(M2)[1], "specimens étudiés") |
ce qui donne 9771 spécimens
head(M2,5) coll_no taxo_no no class family genus 1 5195 55999 373398 Nodosauridae Pawpawsaurus 2 10755 55580 130209 Chaoyangsauridae Chaoyangsaurus 3 10760 38561 144305 Dromaeosauridae 4 10764 66066 130295 Caudipterygidae Caudipteryx 5 10764 66068 130294 Protarchaeopteryx infraclass kingdom order phylum species 1 Chordata Pawpawsaurus campbelli 2 Chordata Chaoyangsaurus youngi 3 Avetheropoda Chordata 4 Avetheropoda Chordata Caudipteryx zoui 5 Avetheropoda Chordata Protarchaeopteryx robusta subclass subfamily subgenus suborder subspecies superclass superfamily 1 2 3 4 5 superphylum tribe unranked clade lng lat max_ma min_ma 1 Ornithischia -97.3000 32.86667 105.3 99.60 2 Ornithischia 123.9667 42.93330 150.8 132.90 3 Saurischia 21.0500 46.11667 70.6 66.00 4 Saurischia 120.7333 41.80000 130.0 122.46 5 Saurischia 120.7333 41.80000 130.0 122.46 |
Voila pour les dinosaures. On peut maintenant chercher des informations sur les plaques tectoniques,
chemin = "data/Shapefile" download.file('https://www.earthbyte.org/webdav/ftp/earthbyte/GPlates/SampleData_GPlates2.0/Individual/FeatureCollections/Coastlines.zip', 'data/coastlines.zip') coast_file = 'Matthews_etal_GPC_2016_Coastlines' unzip(zipfile='data/coastlines.zip', exdir= chemin, junkpaths = TRUE) continents = readOGR(dsn=chemin,layer=coast_file,verbose=TRUE) |
Matthews et al. (2016) a mis en ligne un fichier de rotations simulant la dérive des plaques
download.file('https://www.earthbyte.org/webdav/ftp/earthbyte/GPlates/SampleData_GPlates2.0/Individual/FeatureCollections/Rotations.zip', 'data/rot.zip') rot_file = 'Matthews_etal_GPC_2016_410-0Ma_GK07.rot' unzip(zipfile='data/rot.zip', files = c(paste0('Rotations/',rot_file)), exdir= 'data', junkpaths = TRUE) rot = paste0("data/",rot_file) |
On va corriger quelques anomalies
x = readLines(rot) y = gsub( "!101 !", "!", x ) cat(y, file=rot, sep="\n") remove(x,y) |
et on charge les données en mémoire (pour faire ensuite notre visualisation)
rot_file2 = read.csv(file=rot,header=FALSE,sep='', comment.char = '!') ll = dim(rot_file2)[1] rot_file3 = cbind(rot_file2[1:ll-1,],rot_file2[2:ll,]) names(rot_file3) = c('plate','Start','lat_St', 'lon_St', 'ang_St', 'anchor','plate2','End','lat_End', 'lon_End', 'ang_End', 'anchor2') rot_file = rot_file3[rot_file3$plate==rot_file3$plate2,] rot_file$plate2 = NULL rot_file$anchor2 = NULL PLATES = sort(unique(rot_file$plate)) |
On ne va garder que les dinosaures qui peuvent etre rattaches à une plaque tectonique
X = M2%>%select(lng, lat) Y = SpatialPoints(X,proj4string = continents@proj4string) plaques = over(Y,continents)$PLATEID1 filtre = which(!is.na(plaques)) X0 = X[filtre,] NBX = dim(X0)[1] print(paste0(NBX, " spécimens retenus")) |
ce qui laisse quand même 9617 spécimens
PERIOD = M2%>%select(max_ma, min_ma) PERIOD = PERIOD[filtre,] plaques = plaques[filtre] plaque_id = rep(0,NBX) for (j in seq(1,NBX)){ plaque_id[j]=which(PLATES==plaques[j]) } dataX = M2[filtre,] dataX["plaques"] = plaques |
On y est presque… on va maintenant remonter de -250 millions d’annees a aujourd’hui, en faisant des bonds de 10 millions d’annees
TMAX = 250 TMIN = 0 PAS = 10 |
Pour toutes ces dates, on calcule les quaternions
ROT = array(rep(0,4),c(TMAX/PAS,length(PLATES),4)) for (t in seq(1,TMAX/PAS)){ ROT[t,,] = projecteur(t*PAS, PLATES, rot_file) } QM = list() plaque=list() xy = list() ll = length(continents@polygons) for (i in seq(1,ll)){ M = continents@polygons[[i]]@Polygons[[1]]@coords xy[[i]] = M QM[[i]] = quat(M[,2],M[,1]) plaque[[i]] = which(PLATES==continents$PLATEID1[i]) } QX = quat(X0[,2],X0[,1]) |
On va ensuite projeter
projete=list() Xt = list() X = X0 cpt = 0 for (TIME in seq(TMAX,TMIN,-PAS)){ # setTxtProgressBar(pb, -TIME) cpt = cpt+1 projete[[cpt]] = continents Xt[[cpt]] = X0 if (TIME > 0 ) { for (i in seq(1,ll)){ M = xy[[i]] for (j in seq(1,dim(M)[1])){ M[j,] = rev(rotation(QM[[i]][j,],ROT[TIME/PAS,plaque[[i]],])) } inf_180 = which(M[,1] < -90); inf_180_ = length(inf_180) sup_180 = which(M[,1] > 90); sup_180_ = length(sup_180) if (inf_180_ > 0 & sup_180_ > 0) { if(sup_180_>inf_180_) { M[inf_180,] = t(t(M[inf_180,]) + c(360,0)) } else { M[sup_180,] = t(t(M[sup_180,]) - c(360,0))} } projete[[cpt]]@polygons[[i]]@Polygons[[1]]@coords = M } # setTxtProgressBar(pb, -TIME + PAS/2) for (j in seq(1,NBX)){ if (PERIOD$max_ma[j] > TIME & PERIOD$min_ma[j] <= TIME){ X[j,] = rev(rotation(QX[j,],ROT[TIME/PAS,plaque_id[j],])) } else{ X[j,] = c(NA, NA) } } filtre = which(!is.na(X$lng)) if (length(filtre) > 0){ Xt[[cpt]]=SpatialPointsDataFrame(coords = X[filtre,], data = dataX[filtre,]) } else { Xt[[cpt]]=X } } } |
On peut faire en première carte, 70 millions d’années avant notre ere
t = 1 + (TMAX-70)/PAS leaflet(options = leafletOptions(minZoom = 1)) %>% addPolygons(data=projete[[t]], weight=2) %>% addMarkers(data = Xt[[t]], popup = ~paste(sep = " ", paste("espèce :", Xt[[t]]$species), paste("genre :", Xt[[t]]$genus), paste("famille :", Xt[[t]]$family), paste("ordre : ", Xt[[t]]["unranked clade",]) ), icon = dinoIcon) |
Je mets ici une copie d’écran du leaflet ainsi créé
(l’idée est qu’on peut zoomer, ce qui rend l’analyse plus interactive)
Mais on peut aussi aller 200 millions d’années avant notre ere
t = 1 + (TMAX-200)/PAS leaflet(options = leafletOptions(minZoom = 1)) %>% addPolygons(data=projete[[t]], weight=2) %>% addMarkers(data = Xt[[t]], popup = ~paste(sep = " ", paste("espèce :", Xt[[t]]$species), paste("genre :", Xt[[t]]$genus), paste("famille :", Xt[[t]]$family), paste("ordre : ", Xt[[t]]["unranked clade",]) ), icon = dinoIcon) |
Amusant, non ? en tout cas, merci François pour cette jolie application des quaternions ! Et merci d’avoir suggéré d’utiliser autre chose que des points rouges sur une carte !
dinoIcon = makeIcon(iconUrl = "https://www.ludeek.com/wp-content/uploads/2015/03/uploadfsdfsdf1426350179.1426350368774.png", iconWidth = 30, iconHeight = 50, iconAnchorX = 15, iconAnchorY = 25) |
Le sport en France
Je voulais profiter de la rentree pour mettre en ligne quelques billets sur la data science (comme on dit), en particulier en me basant sur des projets R de la formation en Data Science pour l’Actuariat. L’an passe, j’avais déjà mis en ligne un billet sur le sport (“le sport, une activité de riches“). Cette fois, en m’inspirant de ce qu’a proposé Benoit, on va regarder qui sont les licenciés des différentes fédérations sportives, et ou ils vivent. Comme toujours en R, on charge les librairies qu’on va utiliser…
library(rgdal) library(sp) library(reshape2) library(data.table) library(ggplot2) library(gridExtra) library(ggmap) library(RColorBrewer) library(classInt) library(backports) library(OpenStreetMap) |
J’ouvre une parenthèse rapide, mais en pratique on sait rarement ce qui va servir… ex-post, on va les ramener ce chargement de librairies au début. Je pense que ça serait mieux de les charger juste quand on les utilise. Bon, ensuite, il faut les donnees
Url_Licences = "https://www.data.gouv.fr/s/resources/recensement-des-licences-et-clubs-aupres-des-federations-sportives-agreees-par-le-ministere-charge-d/20180131-163516/Licences_2015.csv" Licences_2015 = read.csv(file=Url_licences, header=TRUE, sep=",",stringsAsFactors = FALSE) Url_Federation = "http://freakonometrics.free.fr/Projet_R/Code_federation.csv" Code_Fede = read.csv(Url_Federation, sep=";",header=FALSE, skip=3) colnames(Code_Fede) = c("Code_Federation","Libelle_Federation") |
On change ici le nom des variables, ça sera plus simple ensuite, et on retient juste quelques lignes interessantes
Code_Fede = Code_Fede[c(1:31,33:92),c(1:2)] |
Il faut ensuite les coordonnées des villes pour faire une carte
Commune = read.csv(file="https://www.data.gouv.fr/fr/datasets/r/554590ab-ae62-40ac-8353-ee75162c05ee", sep=";", header=TRUE) |
En fait, juste la latitude de la longitude nous interesse
Geocod = colsplit(Commune$coordonnees_gps, ",", c("Latitude", "Longitude")) Commune = data.frame(Commune,Geocod) |
Un peu de menage ne fera pas de mal
Commune$Ligne_5 = NULL Commune$coordonnees_gps = NULL doublons = which(duplicated(Commune$Code_commune_INSEE)) #détecte les lignes où il y a doublon Commune_Indiv = Commune[-doublons,] |
On rajoute maintenant un libelle pour chaque sport
Licences_2015 = merge(x=Licences_2015, y=Code_Fede, by.x="fed_2014", by.y="Code_Federation", all.y=TRUE) |
Et on supprime également les lignes ou les codes commune ne sont pas renseignés (car les données ne seront pas exploitables)
Licences_2015 = Licences_2015[!is.na(Licences_2015$newcog2),] |
On a besoin de faire un peu attention a Paris et Marseille, car on a des données par arrondissement,
for (i in 1:nrow(Licences_2015)){ if (Licences_2015[i,c("newcog2")]=="75056") { (Licences_2015[i,c("newcog2")] = "75101")} if (Licences_2015[i,c("newcog2")]=="13055") { (Licences_2015[i,c("newcog2")] = "13101")}} Licences_2015 = merge(x=Licences_2015, y=Commune_Indiv, by.x="newcog2", by.y="Code_commune_INSEE", all.x=TRUE) |
On y est presque. On va créer la variable taux de licenciés (nombre de licences rapporté a la population) pour chaque commune
Licences_2015$Taux_Licencies = ifelse(Licences_2015$pop_2014 != 0,Licences_2015$l_2015/Licences_2015$pop_2014,0) |
Maintenant, on peut jouer ! Ou presque… reste a faire quelques regroupements en fonction de ce qu’on veut représenter.
df_Nb_Lic_Agg_Fed = aggregate(data.frame( Nb_Licence = Licences_2015$l_2015, Nb_hommes = Licences_2015$l_h_2015, Nb_femmes = Licences_2015$l_f_2015, NbLicences_0_4_Ans=Licences_2015$l_0_4_2015, NbLicences_5_9_Ans=Licences_2015$l_5_9_2015, NbLicences_10_14_Ans=Licences_2015$l_10_14_2015, NbLicences_15_19_Ans=Licences_2015$l_15_19_2015, NbLicences_20_29_Ans=Licences_2015$l_20_29_2015, NbLicences_30_44_Ans=Licences_2015$l_30_44_2015, NbLicences_45_59_Ans=Licences_2015$l_45_59_2015, NbLicences_60_74_Ans=Licences_2015$l_60_74_2015, NbLicences_75_Ans=Licences_2015$l_75_2015, Nb_0_4_Ans=Licences_2015$pop_0_4_2014, Nb_5_9_Ans=Licences_2015$pop_5_9_2014, Nb_10_14_Ans=Licences_2015$pop_10_14_2014, Nb_15_19_Ans=Licences_2015$pop_15_19_2014, Nb_20_29_Ans=Licences_2015$pop_20_29_2014, Nb_30_44_Ans=Licences_2015$pop_30_44_2014, Nb_45_59_Ans=Licences_2015$pop_45_59_2014, Nb_60_74_Ans=Licences_2015$pop_60_74_2014, Nb_75_Ans=Licences_2015$pop_75_2014, Pop_femmes=Licences_2015$popf_2014, Pop_hommes=Licences_2015$poph_2014, Pop_Totale=Licences_2015$pop_2014), by = list(Federation = Licences_2015$Libelle_Federation), sum, na.rm = TRUE) |
On peut ainsi calculer le “taux de féminisation” de chaque sport
df_Nb_Lic_Agg_Fed$tx_femmes = ifelse(df_Nb_Lic_Agg_Fed$Nb_Licence!=0,df_Nb_Lic_Agg_Fed$Nb_femmes/df_Nb_Lic_Agg_Fed$Nb_Licence,0) |
ou la répartition par classe d’âge du nombre de licenciés par fédération
df_Nb_Lic_Agg_Fed$Nb_Licence_Norme = df_Nb_Lic_Agg_Fed$NbLicences_0_4_Ans+ df_Nb_Lic_Agg_Fed$NbLicences_5_9_Ans+ df_Nb_Lic_Agg_Fed$NbLicences_10_14_Ans+ df_Nb_Lic_Agg_Fed$NbLicences_15_19_Ans+ df_Nb_Lic_Agg_Fed$NbLicences_20_29_Ans+ df_Nb_Lic_Agg_Fed$NbLicences_30_44_Ans+ df_Nb_Lic_Agg_Fed$NbLicences_45_59_Ans+ df_Nb_Lic_Agg_Fed$NbLicences_60_74_Ans+ df_Nb_Lic_Agg_Fed$NbLicences_75_Ans |
Pour la classe d’age 0-14 ans, on pose alors
df_Nb_Lic_Agg_Fed$Tx_Licences_0_14_Ans = ifelse(df_Nb_Lic_Agg_Fed$Nb_Licence_Norme != 0, (df_Nb_Lic_Agg_Fed$NbLicences_0_4_Ans+df_Nb_Lic_Agg_Fed$NbLicences_5_9_Ans+df_Nb_Lic_Agg_Fed$NbLicences_10_14_Ans)/df_Nb_Lic_Agg_Fed$Nb_Licence_Norme,0) |
et pour la classe d’age 15-29 ans
df_Nb_Lic_Agg_Fed$Tx_Licences_15_29_Ans = ifelse(df_Nb_Lic_Agg_Fed$Nb_Licence_Norme != 0, (df_Nb_Lic_Agg_Fed$NbLicences_15_19_Ans+ df_Nb_Lic_Agg_Fed$NbLicences_20_29_Ans)/df_Nb_Lic_Agg_Fed$Nb_Licence_Norme,0) |
pour la classe d’age 30-44 ans
df_Nb_Lic_Agg_Fed$Tx_Licences_30_44_Ans = ifelse(df_Nb_Lic_Agg_Fed$Nb_Licence_Norme != 0,(df_Nb_Lic_Agg_Fed$NbLicences_30_44_Ans)/df_Nb_Lic_Agg_Fed$Nb_Licence_Norme,0) |
pour la classe d’age 45-59 ans
df_Nb_Lic_Agg_Fed$Tx_Licences_45_59_Ans = ifelse(df_Nb_Lic_Agg_Fed$Nb_Licence_Norme != 0, (df_Nb_Lic_Agg_Fed$NbLicences_45_59_Ans)/df_Nb_Lic_Agg_Fed$Nb_Licence_Norme,0) |
pour la classe d’age 60 ans et plus (on a compris le truc)
df_Nb_Lic_Agg_Fed$Tx_Licences_60_Ans = ifelse(df_Nb_Lic_Agg_Fed$Nb_Licence_Norme != 0, (df_Nb_Lic_Agg_Fed$NbLicences_60_74_Ans+ df_Nb_Lic_Agg_Fed$NbLicences_75_Ans)/df_Nb_Lic_Agg_Fed$Nb_Licence_Norme,0) |
On passe a la détermination des 25 premières fédérations en nombre de licenciés
dt_Nb_Lic_Agg_Fed = data.table(df_Nb_Lic_Agg_Fed) setorder(dt_Nb_Lic_Agg_Fed,-Nb_Licence,na.last=TRUE) dt_Nb_Lic_Agg_Main_Fed = dt_Nb_Lic_Agg_Fed[1:25,] graph1 = ggplot(data=dt_Nb_Lic_Agg_Main_Fed, aes(x=reorder(Federation,Nb_Licence), y=Nb_Licence)) + geom_bar(stat="Identity",fill = "blue")+ geom_text(aes(label=Nb_Licence),check_overlap = TRUE, vjust=0.5, hjust=0, color="blue")+ ggtitle("TOP 25 des fédérations sportives en termes de licenciés")+ ylim(0, 2500000)+ xlab("Fédérations") + ylab("Nombre de licences") graph1+coord_flip() |
On ordonne ensuite par taux de femmes,
setorder(dt_Nb_Lic_Agg_Main_Fed,-tx_femmes,na.last=TRUE) graph2 = ggplot(data=dt_Nb_Lic_Agg_Main_Fed) + aes(x =reorder(Federation,tx_femmes), y = tx_femmes) + geom_bar(stat="Identity",fill = "pink")+ geom_text(aes(label=paste(round(100*tx_femmes, 0), "%", sep="")),check_overlap = TRUE, vjust=0.5, hjust=0.5, color="black")+ xlab("Fédération") + ylab("part des licenciées femmes")+ ggtitle("la pratique sportive féminine par fédération") graph2+coord_flip() |
Et finalement on va regarder par classe d’age
df_Nb_Lic_Agg_Main_Fed = data.frame(dt_Nb_Lic_Agg_Main_Fed) Licence_Age = melt(df_Nb_Lic_Agg_Main_Fed, id=c("Federation"), measured=c("Tx_Licences_0_14_Ans","Tx_Licences_15_29_Ans", "Tx_Licences_30_44_Ans", "Tx_Licences_45_59_Ans","Tx_Licences_60_Ans")) Licence_Age_Clean = Licence_Age[(Licence_Age$variable=="Tx_Licences_0_14_Ans" | Licence_Age$variable=="Tx_Licences_15_29_Ans" | Licence_Age$variable=="Tx_Licences_30_44_Ans" | Licence_Age$variable=="Tx_Licences_45_59_Ans" | Licence_Age$variable=="Tx_Licences_60_Ans"),] dt_Licence_Age_Clean = data.table(Licence_Age_Clean) setorder(dt_Licence_Age_Clean,-variable,na.last=TRUE) setorder(Licence_Age_Clean,variable,na.last=TRUE) graph3 = ggplot(data=Licence_Age_Clean, aes(x=Federation, y=value, fill=variable)) + geom_bar(stat="identity")+ xlab("Fédération") + ylab("répartition par classe d'âge")+ ggtitle("Répartition des licenciés par classe d'âge") graph3+coord_flip()+scale_fill_brewer(palette="Paired") |
A la lecture du graphique ci-dessus, les sports pourraient être classés en 3 catégories :
- les “sports de jeunes” : ceux-ci ont plus de la motié de leurs licenciés âgés de moins de 15 ans : il s’agit de la gymnastique, du judo, du handball, de la natation, ou encore de la voile.
- les “sports de vieux” : on retrouve ici sans surprise la randonnée, le cyclotourisme, le golf, la pétanque, le tir ou encore les sports sous-marins. Ceux-ci voient leurs licenciés avoir plus de 45 ans pour tois quart d’entre eux.
- les “sports pour tous” qui correspondent à ceux qui n’ont pas encore été cités et pour lesquels classes d’âge apparaissent plus équilibrés
Finallement, on peut regarder quelques sports, sur une carte
map.France = get_map(location = c(lon=1.75, lat=46.70), zoom = 6) |
Rugby_2015 = Licence_Max_2015[Licence_Max_2015$fed_2014=="133",] Voile_2015 = Licence_Max_2015[Licence_Max_2015$fed_2014=="128",] Ski_2015 = Licence_Max_2015[Licence_Max_2015$fed_2014=="121",] PetanQ_2015 = Licence_Max_2015[Licence_Max_2015$fed_2014=="242",] Rugby = ggmap(map.France, extent = "normal") + geom_point(aes(x = Longitude, y = Latitude), data = Rugby_2015, colour="red", alpha = 0.5, size=2.0, na.rm=TRUE)+ theme_nothing(legend = TRUE) + theme(legend.position = "bottom")+ ggtitle("Rugby")+ theme(plot.title = element_text(size = 10, face = "bold", hjust=0.5, color="red")) Voile = ggmap(map.France, extent = "normal") + geom_point(aes(x = Longitude, y = Latitude), data = Voile_2015, colour="blue", alpha = 0.5, size=2.0, na.rm=TRUE)+ theme_nothing(legend = TRUE) + theme(legend.position = "bottom")+ ggtitle("Voile")+ theme(plot.title = element_text(size = 10, face = "bold", hjust=0.5, color="blue")) Ski = ggmap(map.France, extent = "normal") + geom_point(aes(x = Longitude, y = Latitude), data = Ski_2015, colour="grey", alpha = 0.5, size=2.0, na.rm=TRUE)+ theme_nothing(legend = TRUE) + theme(legend.position = "bottom")+ ggtitle("Ski")+ theme(plot.title = element_text(size = 10, face = "bold", hjust=0.5, color="grey")) Petanque = ggmap(map.France, extent = "normal") + geom_point(aes(x = Longitude, y = Latitude), data = PetanQ_2015, colour="chocolate3", alpha = 0.5, size=2.0, na.rm=TRUE)+ theme_nothing(legend = TRUE) + theme(legend.position = "bottom")+ ggtitle("pétanque et jeu provençal")+ theme(plot.title = element_text(size = 10, face = "bold", hjust=0.5, color="chocolate3")) grid.arrange(Rugby,Voile,Ski,Petanque, ncol=2, nrow = 2,top="visualisation géographique de sports \n à fort ancrage régional") |
Amusant, non?
(Advanced) R Crash Course, for Actuaries
The fourth year of the Data Science for Actuaries program started this morning. I will be there for the introduction to R. The slides are available online (created with slidify, the .Rmd file is also available)
A (standard) markdown is also available (as well as the .Rmd file). I have to thank Ewen for his help on slidify (especially for the online quizz, and the integration of leaflet maps or the rgl animated graph….)
Scraper pour avoir des infos sur les médecins sur Paris
Allez, quelques codes R pour s’amuser un peu… Je me suis inspiré du projet de Marine pour la formation Data Science pour l’Actuariat. L’idée était de scraper un site qui permet de prendre un rendez-vous chez un médecin, pour voir où se trouvent les médecins (sur Paris) et de voir leurs prix. Histoire de donner un avant goût, voilà les graphiques qu’on cherche à produire
Continue reading Scraper pour avoir des infos sur les médecins sur Paris
Les votes à l’assemblée nationale
Un dernier petit billet basé sur les projets de R que j’avais donnés pour la formation Data Science pour l’Actuariat. Aujourd’hui, je reviens sur des codes tirés du projet de Raphaël qui a scrapé les données de l’Assemblée Nationale. On commence par charger les librairies qui nous seront utiles,
require(xml2)
require(downloader)
require(stringr)
require(classInt)
require(plotrix)
require(FactoMineR)
require(sp)
La première partie de l’importation des données vise à importer l’ensemble des députés, c’est à dire l’ensemble des personnes (compris entre 1 et 248) dont le mandat est à l’ASSEMBLEE. Pour chacun de ces députés, on récupére le numéro de circonscription et de département, afin de pouvoir construire le code IDEN, identifiant unique de chaque circonscription (DDCC, deux chiffres pour le département, deux pour la circonscription), la référence acteur (acteurRef) et la référence mandat (uid) afin de pouvoir construire un identifiant pour chaque député du type (PA???XXXPMXXXXXX, chaque X representant un chiffre, les ? representant un chiffre ou l’absence de caractères – les PA n’ont pas tous le même nombre de caractères). On va ensuite créer des variables qui serviront par la suite : le nombre de vote total, le nombre de vote POUR, CONTRE, ABSTENTION et le nom de l’organe (groupe auquel appartient le deputé à l’assemblée).
path="http://data.assemblee-nationale.fr/static/openData/repository/AMO/deputes_senateurs_ministres_legislature/AMO20_dep_sen_min_tous_mandats_et_organes_XIV.csv.zip"
dest="deputes.zip"
download(path,destfile=dest,mode="wb")
loc=paste(getwd(),"/",dest,sep="" )
unzip(loc)
dest="acteurs.csv"
loc=paste(getwd(),"/",dest,sep="" )
t2=read.csv(loc,sep=";")
res=NULL
for (i in 1:248){
test=str_c("mandats.1..mandat.",i,"..typeOrgane.1.")
j=which(colnames(t2) == test)
#Je ne récupère que les personnes dont le i ème mandat est un mandat de député i.e siegeant à l'Assemblée Nationale
t3=subset(t2,t2[,j]=="ASSEMBLEE")
if(nrow(t3)!=0){
#Nom des différentes variables en concatenant les chaines
circo=str_c("mandats.1..mandat.",i,"..election.1..lieu.1..numCirco.1.")
dept=str_c("mandats.1..mandat.",i,"..election.1..lieu.1..numDepartement.1.")
acteur=str_c("mandats.1..mandat.",i,"..acteurRef.1.")
mandat=str_c("mandats.1..mandat.",i,"..uid.1.")
# Récupération du numero de la colonne du data frame dont le nom correspond a chacune des variable voulues
k=which(colnames(t3) == circo)
l=which(colnames(t3) == dept)
m=which(colnames(t3) == acteur)
n=which(colnames(t3) == mandat)
#Je cree un data frame correspondant "au numero de mandat"
t4=data.frame(as.numeric(as.character(t3[,k])),as.numeric(as.character(t3[,l])),str_c(t3[,m],t3[,n]))
colnames(t4)=c("circo","dept","identifiant")
#Et j'ajoute l'ensemble des députés
res=rbind(res,t4)
}
}
res$iden=str_c(res$dept,res$circo)
res$nbVote=0
res$oui=0
res$non=0
res$abst=0
res$nomOrgane=0
Ensuite, pour récupérer les votes de chacun des députés à chaque scrutin, on récupére le fichier XML du site data.assemblee-nationale.fr.
path="http://data.assemblee-nationale.fr/static/openData/repository/LOI/scrutins/Scrutins_XIV.xml.zip"
dest="Scrutins_XIV.xml.zip"
download(path,destfile=dest,mode="wb")
loc=paste(getwd(),"/",dest,sep="" )
unzip(loc)
dest="Scrutins_XIV.xml"
loc=paste(getwd(),"/",dest,sep="" )
En parcourant la hiérarchie du fichier XML (à l’aide de la fonction xml_children) on obtient les différents fils d’un noeud.
t=read_xml(loc)
liste=xml_children(t)
Le premier niveau de noeud correspond à l’ensemble des scrutins. En bouclant sur chacun d’eux, on récupéree dans les fils le résultat du vote de chaque scrutin pour chaque député. On va se limiter aux scrutins dont le mode de publication scrutin (14ème fils) est “Decompte nominatif” : en effet, pour les scrutins par “Decompte Dissident”, il n’y avait que les votes des députés qui n’avaient pas voté dans le sens la majorité, ce qui introduirait un biais ensuite.
Le 16ème fils correspond à la variable ventilation du vote. Ce noeud recense toutes les données concernant le résultat du vote.
Pour chacun des groupes représentés à l’Assemblée présent dans le fils ventilation vote, on récupère le vote de chaque député du groupe dont les modalités sont : “Non Votant”, “Pour”, “Contre” et “Abstention”.
En bouclant sur l’ensemble des scrutins on comptabilise également le nombre de vote réalisé pour chaque député et le nom de son groupe (organeRef).
for (a in liste){
if(xml_text(xml_children(a)[14])=="DecompteNominatif"){
numero=xml_text(xml_children(a)[2])
ventil=xml_children(a)[16]
groupe=xml_children(xml_children(xml_children(ventil)))
tempDataFrame=data.frame(res$identifiant,NA)
colnames(tempDataFrame)=c("identifiant",str_c("scrutin",numero))
for (b in groupe){
intermediaire=xml_children(xml_children(b))
nomGroupe=xml_text(xml_children(b)[1])
for (i in 1:4){
df3=data.frame(xml_text(xml_children(xml_children(intermediaire[3])[i])))
if(nrow(df3)!=0){
for (j in df3[,]){
if (i==1){
j=strsplit(strsplit(strsplit(j,"MG")[[1]],"PSE")[[1]],"PAN")[[1]]
}
res[res$identifiant==j,]$nomOrgane=as.character(nomGroupe)
if(i!=1)
res[res$identifiant==j,]$nbVote=res[res$identifiant==j,]$nbVote+1
if (i==1)
tempDataFrame[tempDataFrame$identifiant==j,2]="NV"
if (i==2){
res[res$identifiant==j,]$oui=res[res$identifiant==j,]$oui+1
tempDataFrame[tempDataFrame$identifiant==j,2]="POUR"
}
if (i==3){
res[res$identifiant==j,]$non=res[res$identifiant==j,]$non+1
tempDataFrame[tempDataFrame$identifiant==j,2]="CTRE"
}
if (i==4){
res[res$identifiant==j,]$abst=res[res$identifiant==j,]$abst+1
tempDataFrame[tempDataFrame$identifiant==j,2]="ABST"
}}}}}
res=data.frame(res, tempDataFrame[,2])
}}
res2=subset(res,res$dept<96)
res2$circo=ifelse(nchar(res2$circo)==1, str_c("0",res2$circo) , str_c("",res2$circo))
res2$dept=ifelse(nchar(res2$dept)==1, str_c("0",res2$dept) , str_c("",res2$dept))
res2$iden=str_c(res2$dept,res2$circo)
res2=subset(res2,!is.na(iden))
res2=res2[order(res2$iden),]
parIden=aggregate(res2$nbVote, by=list(res2$iden), sum)
Par soucis de visualisation, on va se limiter ici à la France métropolitaine. Il a fallu uniformiser le code IDEN de chaque département pour qu’il corresponde au format utilisé par la carte de France utilisé (DDCC, deux chiffres pour le département, deux pour la circonscription). Le fond de carte est le suivant
path="http://www.laspic.eu/data/documents/circosshp_v3.zip"
dest="circosshp_v3.zip"
download(path,destfile=dest,mode="wb")
loc=paste(getwd(),"/",dest,sep="" )
unzip(loc)
dest="circosSHP_v3.RData"
loc=paste(getwd(),"/",dest,sep="" )
load(loc)
On agrége enfin par par député afin d’obtenir un nombre de vote par circonscription. Notons que dans certaines circonscriptions, le député démissionait ou mourrait et était donc remplacé par un autre député. J’ai donc choisi de sommer le vote de l’ensemble des députés ayant représenté une circonscription pour le représenter graphiquement : on raisonne par circonscription, et par par représentant.
nuancier <- findColours(classIntervals(parIden$x, 6, style = "quantile"), smoothColors("white",98,"#0C3269"))
plot(fdc, col=nuancier)
leg <- findColours(classIntervals( round(parIden$x), 6, style="quantile"), smoothColors("white",98,"#0C3269"), under="moins de", over="plus de", between="–",cutlabels=FALSE)
legend("bottomleft",fill=attr(leg, "palette"), legend=names(attr(leg,"table")),title = "Nombre de Votes",bty="n")
title( main="Nombre de votes par circonscription",cex.main=1.5)
parIdenScrutin=aggregate(res2[,9:ncol(res2)], by=list(res2$iden),na.omit)
for (i in 1:nrow(parIdenScrutin))
parIdenScrutin$nomOrgane[i]=ifelse(parIdenScrutin$nomOrgane[i][[1]][1]!="0",as.character(parIdenScrutin$nomOrgane[i][[1]][1]),as.character(parIdenScrutin$nomOrgane[i][[1]][2]))
for (i in 3:ncol(parIdenScrutin))
parIdenScrutin[,i]=as.factor(levels(parIdenScrutin[,i][[1]])[as.numeric(parIdenScrutin[,i])])
parIdenScrutin=replace(parIdenScrutin,is.na(parIdenScrutin),"NV")
for (i in 1:nrow(parIdenScrutin)) {
if (parIdenScrutin$nomOrgane[i]=="PO656014" || parIdenScrutin$nomOrgane[i]=="PO713077" || parIdenScrutin$nomOrgane[i]=="PO656002")
parIdenScrutin$nomOrgane[i] <- "SER"
if (parIdenScrutin$nomOrgane[i] == "PO656006" || parIdenScrutin$nomOrgane[i] =="PO707869")
parIdenScrutin$nomOrgane[i] <- "LR"
if (parIdenScrutin$nomOrgane[i] == "PO656022")
parIdenScrutin$nomOrgane[i] <- "RRDP"
if (parIdenScrutin$nomOrgane[i] == "PO656010")
parIdenScrutin$nomOrgane[i] <- "UDI"
if (parIdenScrutin$nomOrgane[i] == "PO656018")
parIdenScrutin$nomOrgane[i] <- "GDR"
if (parIdenScrutin$nomOrgane[i] == "PO645633")
parIdenScrutin$nomOrgane[i] <- "NI"
}
parIdenScrutin$Group.1 <- as.factor(parIdenScrutin$Group.1)
parIdenScrutin$nomOrgane <- as.factor(as.character(parIdenScrutin$nomOrgane))
On peut finir par une petite analyse des correspondances : les députés ont voté a 644 scrutins (en colonne). Pour chaque scrutin, un député peut avoir voté pour (POUR), contre (CTRE), s’être abstenu (ABST) ou ne pas s’être présenté. Si on ne garde que les deux premiers axes pricinpaux, on obtient
acm <- MCA(parIdenScrutin, quali.sup=2, graph=FALSE)
head(acm$eig,5)
plot(acm$ind$coord[, 1:2], type= "n", xlab=paste0("Axe 1 (" , round(acm$eig[1,2], 1), " %)"), ylab=paste0("Axe 2 (", round(acm$eig[2,2], 1), " %) "), main= "Nuage des individus selon les partis", cex.main=1, cex.axis=1, cex.lab=1, font.lab=3)
abline(h=0, v=0, col= "grey", lty=3, lwd=1)
points(acm$ind$coord[,1:2], col=as.numeric(parIdenScrutin$nomOrgane), pch=19, cex=0.5)
legend("topleft", legend=levels(parIdenScrutin$nomOrgane), bty= "o", text.col=1:10, col=1:10, pch=18, cex=0.8)
text(acm$ind$coord[,1:2], labels=rownames(acm$ind$coord), col=as.numeric(parIdenScrutin$nomOrgane), cex=0.7, pos=4)
On a ici le nuage des individus. Ces derniers sont placés en fonction de leur tendance de votes. Les individus les plus proches sont ceux qui votent sensiblement de la même façon tandis que les individus éloignés votent différemment à chaque scrutin. Les couleurs représente le parti auquel appartient le député. Elles sont présentent à titre indicatif pour valider ou non l’hypothèse de “vote en groupes“. Sur ce graphique, parmi les groupes représentés, deux sont très nets à gauche et à droite. Un 3ème se remarque dans la partie centre haute du graphique et enfin un 4e groupe dans la supérieure droite du graphique.
Si on se limite aux couleurs (les partis), les deux groupes très nets qui s’opposent sont Les Républicains (LR – Droite, à gauche sur la projection) et le groupe Socialiste, Ecologiste et Républicain (SER – Gauche, à droite sur la projection). Le 3ème groupe qui se distingue représente le groupe d’extrème Gauche GDR (Gauche démocrate et républicaine). Ces trois groupes sont disctinct car ils votent de façon opposé. Ensuite en regardant toujours les partis, on remarque que l’Union des Démocrates et Indépendants (UDI – Droite) ont des votes proches de ceux du groupe LR. Le groupe Radical, Républicain, démocrate et progressiste ont une façon de voter similaire à SER. Les votes des députés non inscrits varient énormément entre les scrutins, ils ne se rapporchent d’aucun parti.
Accidents de voiture et de deux roues à Paris
Il y a quelques mois, j’avais commencé à reprendre des codes que j’avais vu passer dans les projets de R que j’avais donnés pour la formation Data Science pour l’Actuariat. J’avais commencé avec un billet sur sur cartographier le vote pour le Brexit (qui avait été repris – et amélioré – sur le site des voisins de rgeomatic). Aujourd’hui, en rangeant un peu mon ordinateur, je suis tombé sur un code de Quentin sur les acccidents de la route à Paris.
Comme toujours on charge quelques packages
require(downloader)
require(ggmap)
require(stringr)
options(encoding = "UTF-8")
et on récupère les données, en ligne sur https://opendata.paris.fr
path="https://opendata.paris.fr/explore/dataset/accidentologie/download/?format=csv&timezone=Europe/Berlin&use_labels_for_header=true"
dest="test"
download(path,destfile=dest,mode="wb")
loc=paste(getwd(),"/",dest,sep="" )
obj <- read.csv(file = loc, stringsAsFactor = FALSE,sep=";")
file.remove(dest)
On peut commencer par regarder quand on eu lieu les accidents, quelle jour et à quelle heure. On va alors agréger les données
t=data.frame(obj,counter=1)
parJour=aggregate(t$counter, by=list(format(as.Date(t$Date),"%A")), sum)
et créer une base qui contient toutes les informations
vec_dat=data.frame();
vec_dat[ (1:3),1]=parJour$Group.1[3:5]
vec_dat[ (1:3),2]=parJour$x[3:5]
vec_dat[ 4,1]=parJour$Group.1[2]
vec_dat[ 4,2]=parJour$x[2]
vec_dat[ 5,1]=parJour$Group.1[7]
vec_dat[ 5,2]=parJour$x[7]
vec_dat[ 6,1]=parJour$Group.1[6]
vec_dat[ 6,2]=parJour$x[6]
vec_dat[ 7,1]=parJour$Group.1[1]
vec_dat[ 7,2]=parJour$x[1]
on va rajouter un nom, pour que ce soit plus clair
names(vec_dat)[1]="Jour";names(vec_dat)[2]="Nbsinistre"
On choisit éventuellement une palette de couleurs
Palette<-colorRampPalette(c("lightblue","darkblue"))
et on fait un barplot, pour compter le nombre d’accidents, par jour de la semaine, sur Paris
barplot(vec_dat$Nbsinistre,names.arg=vec_dat$Jour,main="Répartition du nombre d'accidents par jour de la semaine ", col=Palette(10), xlab="jour de la semaine", ylab="nombre d'accidents")
On peut ensuite regarder par heure de la journée (en données brutes, sans lissage)
parHeure=aggregate(t$counter, by=list(format(strptime(t$Heure,format="%H:%M"),"%H")), sum)
plot(parHeure,type="l", col=c("blue"), main=" Répartition du nombre d'accidents par heure de la journée", xlab="heure de la journee", ylab="nombre d'accidents")
Passons maintenant à de la géolocalisation
t2 = subset(t, Coordonnées!= "")
coord=data.frame(t(matrix(as.numeric(unlist(strsplit(t2$Coordonnées,"[,]"))),nrow=2)))
t2$lat=coord[,1]
t2$lon=coord[,2]
t2 = subset(t2, ((48.81 <= lat) & (lat<=48.90) & (2.25<=lon) &(lon<=2.43) ) )
On peut regarder sur une carte où ont eu lieu les accidents
map<-get_map(location = "Paris", zoom=12, maptype="roadmap", color="bw")
vis <- ggmap(map) + geom_point(data = t2, aes(x=lon, y=lat, color=1)) + theme( legend.position="none") + labs(title="Cartographie des accidents à Paris")
vis
On peut tenter une densité (brute, non corrigée d’une densité de routes, ou de population)
overlay<-stat_density2d(data=t2, aes(x=lon, y=lat, fill=..level.., alpha=..level..), contour=T, n=100, geom="polygon")
densi<- ggmap(map) + overlay + scale_fill_gradient("Accident Density") +scale_alpha(range=c(0.4,0.75),guide=FALSE) + guides(fill=guide_colorbar(barwidth=1.5, barheight = 10))
densi = densi + labs(title="Densite des accidents à Paris")
densi
et non-corrigée ici d’effets de bords (ce qui donne l’illusion qu’il y a plus d’accidents intra-muros que sur le boulevard périphérique par exemple). On peut également regarder en fonction de variables explicatives, comme le type de véhicules impliqués, ou la gravité
densi_byGravUs<- densi + facet_wrap(~ Usager.1.Grav, nrow=2)
densi_byGravUs = densi_byGravUs + labs(title="Densité des accidents à Paris par gravité d'accident")
densi_byGravUs
Bref, on peut faire facilement plein de chose à partir de ces données.
Install R Packages on the Ubuntu Virtual Machine
For the (Advanced) R Crash Course of the Data Science for Actuaries program, we will use the Ubuntu virtual machine. There might be some issues when installing some packages… One trick can be to open a terminal
and then to use the sudo command, to install some packages,
(after entering the password). Just type (or copy/paste)
sudo apt-get install libcurl4-openssl-dev
sudo apt-get install libxml2-dev
sudo apt-get install openjdk-8-*
update-alternatives --config java
sudo apt-get install aptitude
sudo aptitude install libgdal-dev
sudo aptitude install libproj-dev
sudo apt-get install build-essential libcurl4-gnutls-dev libxml2-dev libssl-dev
sudo apt-get install curl
sudo apt-get build-dep r-cran-rgl
sudo apt-get install r-cran-plyr r-cran-xml r-cran-reshape r-cran-reshape2 r-cran-rmysql
sudo apt-get install r-cran-rjava
sudo apt-get install r-cran-glmnet
sudo apt-get build-dep r-cran-boot
sudo apt-get build-dep r-cran-class
sudo apt-get build-dep r-cran-cluster
sudo apt-get build-dep r-cran-codetools
sudo apt-get build-dep r-cran-foreign
sudo apt-get build-dep r-cran-kernsmooth
sudo apt-get build-dep r-cran-lattice
sudo apt-get build-dep r-cran-mass
sudo apt-get build-dep r-cran-matrix
sudo apt-get build-dep r-cran-mgcv
sudo apt-get build-dep r-cran-nlme
sudo apt-get build-dep r-cran-nnet
sudo apt-get build-dep r-cran-rpart
sudo apt-get build-dep r-cran-spatial
sudo apt-get build-dep r-cran-survival
sudo apt-get build-dep r-cran-rodbc
sudo apt-get build-dep
Then, in RStudio, enter
install.packages("RCurl")
install.packages("xml")
install.packages("rJava")
install.packages("rgdal")
install.packages("xlsx")
install.packages("devtools")
It should be fine…