Category Archives: Computer

Tous des (potentiels) terroristes ?

En fin de semaine dernière, je commencais surveiller et prévenir: L’ère de la pénalité prédictive de Nicolas Bourgoin, et au début, il évoque des mesures de sécurité renforcées dans un rayon de 20km des gares, des aéroports et des ports, qui pourraient être prises dans le cadre de la réforme de lois sur la sécurité intérieure.  Cette information était reprise récemment sur rtl,

Le périmètre est aussi agrandi : les contrôles pourront avoir lieu aux abords des gares internationales (et non plus seulement à l’intérieur) ainsi que dans un rayon de 20 kilomètres autour des aéroports et des ports.

ou sur le site de l’obs

ces opérations de contrôles seront mises en place “aux abords” de 373 gares, ports et aéroports, dans un rayon de 20 kilomètres. Une extension considérable puisque jusqu’à présent, ces contrôles restaient cantonnées à l’intérieur de ces espaces accessibles au public.

En lisant le livre, je trouvais incroyable cette histoire de 20km, car tout le monde habite à moins de 20km d’une gare (ne sachant trop quelles pouvaient être ces 373 gares, je pars du fait que toutes les gares pourraient être concernées). C’est un peu ce que note la cimade, en écrivant

Le projet de loi prévoit ainsi de permettre les contrôles d’identité aux frontières pour une durée de 12 heures (contre 6 aujourd’hui), de les élargir « aux abords » de 373 gares, ports et aéroports, ainsi que dans un rayon de 20 km des 118 points de passages frontaliers. Bien au-delà des simples frontières de l’Hexagone, c’est presque tout le territoire qui est couvert. Le dispositif porte ainsi une atteinte disproportionnée aux droits et libertés des personnes. Des villes entières, comme Paris et toute la région Île-de-France, Lyon, Nantes, Rennes, Bordeaux, Montpellier, Toulouse ou Marseille seraient soumises à un régime de légalisation du contrôle au faciès. Des personnes assimilées par la police comme étant étrangères, quelle que soit leur situation en France, risquent ainsi d’être les victimes de ces contrôles d’identité.

Etant d’un naturel (très) sceptique, j’ai voulu vérifier par moi même, non seulement en terme de surface (ce qui est dit ici) mais surtout en terme de population. La première étape est de récupérer la liste des gares géolocalisées, sur https://ressources.data.sncf.com/. On peut aussi récupérer la liste des aéroports sur https://www.data.gouv.fr/ si on veut aller plus loin. Mais contentons nous des gares pour l’instant. La seconde étape est de récupérer les contours des communes et leur population. En fait, plus que la superficie, ce qui m’intéresse le plus, c’est le nombre de personnes qui y habitent. On peut trouver un tel fichier sur https://www.data.gouv.fr/ là encore. Mais commencons par charger tous nos packages de cartographie,

library(maptools)
library(rgeos)
library(rgdal)
library(ggplot2)
library(plyr)
library(maptools)

On peut ensuite récupérer les données des gares

loc = "/home/arthur/referentiel-gares-voyageurs.shp"
gare = readOGR(loc)

Pour regarder où sont ces gares, on récupère un fond de carte

Proj = CRS("+proj=longlat +datum=WGS84")
France = readShapePoly("departements-20140306-100m.shp", verbose=TRUE, proj4string=Proj)
metropolitaine = as.character(1:95)
metropolitaine[1:9] = paste("0",1:9,sep="")
France = France[France$code_insee%in%metropolitaine,]

On utilise aussi une base avec les contours des communes,

loc2 = "/home/arthur/communes-20150101-5m.shp"
communes_lim = readOGR(loc2)
communes_lim = spTransform(communes_lim, CRS("+proj=longlat +datum=WGS84"))

et une base pour la population des communes

 base = read.csv(
"http://freakonometrics.free.fr/popfr19752010.csv",
header=TRUE)
base$insee = base$dep*1000+base$com

Cette fois on est prêt. Considérons un rayon de 20km

r=20

On va constituer tous les polygônes correspondant à une cercle de rayon 20km, centré sur les gares francaises

PL = data.frame(i=1:(3000*120),id=rep(1:3000,each=120),lon=NA,lat=NA)
 for(i in 1:nrow(gare)){
   x=as.numeric(as.character(gare$longitude_w[i]))
   y=as.numeric(as.character(gare$latitude_wg[i]))
   vx=c(x+u*r/111,x+rev(u)*r/111)
   vy=c(y+v*r/111,y-rev(v)*r/111)
   polygon(vx,vy,border=NA,col="blue",pch=19)
   PL[PL$id==i,2:4]=data.frame(id=i,lon=vx,lat=vy) }

(comme 1 degré fait environ 111km). On doit alors bricoler un peu pour constituer une collection de polygones, que l’on pourra ensuite manipuler à notre guise

PL=PL[!is.na(PL$lat),]
 PLdf=PL[,c(3,4,2)]
 PLdf[,3]=as.factor(PLdf[,3])
 PL_list <- split(PLdf, PL$id)
 PL_list <- lapply(PL_list, 
  function(x) { x["id"] <- NULL; x })
 PPL <- sapply(PL_list, Polygon)
 PPL1 <- lapply(seq_along(PPL), 
    function(i) Polygons(list(PPL[[i]]),
        ID = names(PL_list)[i]  ))
PPLs <- SpatialPolygons(PPL1, 
proj4string = CRS("+proj=longlat +datum=WGS84") )

Si on visualise ces périmètres de sécurité, on obtient

G=gUnionCascaded(PPLs)
F=gUnionCascaded(France)
FG=gIntersection(F,G)
plot(F)
plot(G,add=TRUE,col="blue")

soit, si on regarde l’intersection entre les deux,

plot(F)
plot(FG,add=TRUE)

Pour trouver la population dans cette vaste région, on va supposer la population uniformément répartie sur une commune. Et on regarde la proportion de la surface de la commune qui est à moins de 20km d’une gare (c’est en gros la technique qu’on utilisait dans Kernel density estimation based on Ripley’s correction pour corriger d’effet de bords dans du lissage spatial, sur des cartes, avec Ewen Gallic). Bref, pour une commune donnée, on a le code suivant, qui permet de récuperer la proportion de la surface située à moins de 20km d’une gare, et sa population

f=function(i){
insee=as.numeric(as.character(communes_lim@data$insee[i]))
  POPULATION=base[base$insee==insee,"pop_2010"]
  B_list=list()
  for(j in 1:length(communes_lim@polygons[[1]]@Polygons)){ B_list[[j]]=data.frame(communes_lim@polygons[[i]]@Polygons[[j]]@coords,id=j)}
  B_list <- lapply(B_list,function(x) { x["id"] <- NULL; x })
  BL <- sapply(B_list, Polygon)
  BL1 <- lapply(seq_along(BL), function(i) Polygons(list(BL[[i]]),ID = names(PL_list)[i]  ))
  BLs <- SpatialPolygons(BL1, proj4string = CRS("+proj=longlat +datum=WGS84") ) 
  t=try(FGB<-gIntersection(BLs,FG),silent=TRUE)  t1=try(l<-length(BLs@pointobj@coords),silent=TRUE)
  if((!inherits(t1, "try-error"))){ 
    a_list=list()
    for(j in 1:length(BLs@pointobj@coords)){
      a_list[[j]]=BLs@polyobj@polygons[[j]]@area}
    a1=sum(unlist(a_list))}
  if((inherits(t1, "try-error"))){ 
    a1=BLs@polygons[[1]]@area}
  a2=0
  if(!is.null(FGB)){
t2=try(l<-length(FGB@pointobj@coords),silent=TRUE) if((!inherits(t2, "try-error"))){ a_list=list() for(j in 1:length(FGB@pointobj@coords)){ a_list[[j]]=FGB@polyobj@polygons[[j]]@area} a2=sum(unlist(a_list))} if((inherits(t2, "try-error"))){ a2=FGB@polygons[[1]]@area}} p=c(1,NA,0) if((!inherits(t, "try-error"))&(!is.null(t))&(length(POPULATION)==0)) p=c(a2/a1,a1,0) if((!inherits(t, "try-error"))&(!is.null(t))&(length(POPULATION)>0) ) p=c(a2/a1,a1,POPULATION) 
cat(i,insee,sum(unlist(lapply(B_list,nrow))),p,"\n")
  return(p)}

On n’a plus qu’à faire tourner sur notre 36000 communes,

F = lapply(unique(communes_lim_df$id),f)

On obtient alors

> F0=F[!is.na(F[,2]),]
> sum(F0[,1]*F0[,2])/sum(F0[,2])
[1] 0.9031736
> F=F1
> sum(F[,1]*F[,3])/sum(F[,3])
[1] 0.9661656

Autrement dit, dans un rayon de 20km des gares (seulement), on a 90% du territoire, et plus de 95% de la population. Si on passe à 10km, on note que l’on couvre environ 75% du territoire, et plus de 90% de la population, et en passant à 5km, on couvre moins de 50% de la surface du territoire, et environ 75% de la population.

On retrouve des résultats du même ordre que ceux obtenus dans un précédant billet, qui établissait que 80% de la population francaise était à moins de 3km d’une agence bancaire. On pourrait bien sûr rajouter les ports, et les aéroports, et surprimer quelques gares… mais je doute qu’on ait une conclusion très différente… Après, les codes sont disponibles, il suffit de les adapter sur une autre base de centres de cercles…

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

I Got The Feelin’

Last week, I’ve been going through my CD collection, trying to find records I haven’t been listing for a while. And I got the feeling that music I listen to nowadays is slower than the one I was listening to in my 20’s. I was wondering if that was an age issue, or it was simply the fact that music in the 90s was “faster” than the one released in 2015. So I tried to scrap the BPM database to get a more appropriate answer about this “feeling” I have. I did extract two information: the BPM (beat per minute) and the year (of release).

Here is a function to extract information from the website,

> library(XML)
> extractbpm = function(VBP,P){
+ url=paste("https://www.bpmdatabase.com/music/search/?artist=&title=&bpm=",VBP,"&genre=&page=",P,sep="")
+ download.file(url,destfile = "page.html")
+ tables=readHTMLTable("page.html")
+ return(tables)}

For instance

> extractbpm(115,13)
$`track-table`
Artist Title
1 Eros Ramazzotti y Claudio Guidetti Dimelo A Mi
2 Everclear Volvo Driving Soccer Mom
3 Evils Toy Dear God
4 Expose In Walked Love
5 Fabolous ft. 2 Chainz When I Feel Like It
6 Fabolous ft. 2 Chainz When I Feel Like It
7 Fabolous ft. 2 Chainz When I Feel Like It
8 Fanny Lu Fanfarron
9 Featurecast Ain't My Style
10 Fem 2 Fem Obsession
11 Fernando Villalona Mi Delito
12 Fever Ray Triangle Walks
13 Firstlove Freaky
14 Fito Blanko Pegadito Suavecito
15 Flechazo Del Norte Mariposa Traicionera
16 Fluke Switch/Twitch
17 Flyleaf Something Better
18 FM Static The Next Big Thing
19 Fonseca Eres Mi Sueno
20 Fonseca ft. Maffio & Nayer Eres Mi Sueno
21 Francesca Battistelli Have Yourself A Merry Little Christmas
22 Frankie Ballard Young & Crazy
23 Frankie J. More Than Words
24 Frank Sinatra The Hucklebuck
25 Franz Ferdinand The Dark Of The Matinée
Mix BPM Genre Label Year
1 — 115 — Sony 2009
2 — 115 — Capitol Records 2003
3 — 115 — — —
4 — 115 — Arista Records 1994
5 Explicit 115 Urban Def Jam/Island Def Jam 2013
6 — 115 Urban Def Jam/Island Def Jam 2013
7 Radio Edit 115 Urban Def Jam/Island Def Jam 2013
8 — 115 Latin Pop Universal Latino 2011
9 Psychemagik Dub 115 — Jalapeno 2012
10 — 115 — Critique Records 1993
11 — 115 — Mt&vi Records/caminante Records 2001
12 Rex The Dog Remix 115 — Little Idiot/Mute 2012
13 — 115 — Jwp Music 2000
14 — 115 Merengue Mambo Crown Loyalty 2012
15 — 115 — Hacienda 2010
16 Album Version 115 — One Little Indian Records 2004
17 — 115 Alternative A&M/Octone 2013
18 — 115 — Tooth & Nail Records 2007
19 — 115 Merengue Mambo 10 2012
20 Urban Version 115 — 10 2012
21 — 115 — Word/Fervent/Warner Bros 2009
22 — 115 Country Warner Bros 2015
23 Mynt Rocks Radio Edit 115 — Columbia 2005
24 — 115 Jazz Columbia 1950
25 — 115 New Wave — 2004

We have here one of the few old songs, a 1950 tune by Frank Sinatra. If we scrap the website, with a simple loop (where the bpm is from 40 to 200). Start with

BASE=NULL
> vbp=40
> p=1

and then, a loop based on

> while(vbp<=2017){
+ F=extractbmp(vbp,p)
+ if(length(F)==1){
+ BASE=rbind(BASE,F[[1]][,c("Artist","Title","BPM","Year")])
+ p=p+1}
+ if(length(F)==0){
+ vbp=vbp+1
+ p=1}}

Then we should clean the dataset

BASE=BASE[-duplicated(BASE),]
BASE=BASE[-which(BASE$Year=="—"),]
BASE$y=as.numeric(as.character(BASE$Year))
BASE$bpm=as.numeric(as.character(BASE$BPM))
BASE=BASE[BASE$y>=1940,]

and we end up with almost 50,000 tunes.

boxplot(BASE$bpm~as.factor(BASE$y),
col="light blue")

Over the past 20 years, it looks like speed of tunes has declined (let us forget tunes of 2017, clearly we have a problem here…)

library(mgcv)
plot(BASE$y,BASE$bpm)
reg=gam(bpm~s(y),data=BASE)
B=data.frame(y=1950:2017)
p=predict(reg,newdata=B)
lines(B$y,p,lwd=3,col="red")

which is confirmed with a (smoothed) regression

p2=predict(reg,newdata=B,se.fit=TRUE)
plot(B$y,p2$fit,lwd=3,col="red",type="l",ylim=c(90,140))
lines(B$y,p2$fit+p2$se.fit)
lines(B$y,p2$fit-p2$se.fit)

even when incorporating the confidence band. Bumps are probably related to smoothing parameters, but indeed, it looks like the average speed of music tune has decreased, from 110-115 in the 90’s to less than 100 nowadays. Now to be honest, I would love to have access to personal information from itunes, deezer or spotify, to get a better understanding (eg when in the week, in the day, do we like to listen to faster music for instance). But so far, I could not have access to such data. Too bad…

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.

Visualiser les résultats du premier tour

Juste quelques lignes de code, pour visualiser les résultats du premier tour en France. L’idée est de faire une carte assez minimaliste, avec des cercles centrés sur les centroïdes des départements. On commence par récupérer les données pour le fond de carte, un fichier 7z sur le site de l’ign.

download.file("https://wxs-telechargement.ign.fr/oikr5jryiph0iwhw36053ptm/telechargement/inspire/GEOFLA_THEME-DEPARTEMENTS_2016$GEOFLA_2-2_DEPARTEMENT_SHP_LAMB93_FXX_2016-06-28/file/GEOFLA_2-2_DEPARTEMENT_SHP_LAMB93_FXX_2016-06-28.7z",destfile = "dpt.7z")

On a dans ce fichier des informations sur les centroïdes

library(maptools)
library(maps)
departements<-readShapeSpatial("DEPARTEMENT.SHP")
plot(departements)
points(departements@data$X_CENTROID,departements@data$Y_CENTROID,pch=19,col="red")

Comme ça ne marche pas très bien, on va le refaire à la main, par exemple pour l’Ille-et-Vilaine,

pos=which(departements@data[,"CODE_DEPT"]==35)
Poly_35=departements[pos,]
plot(departements)
plot(Poly_35,col="yellow",add=TRUE)
departements@data[pos,c("X_CENTROID","Y_CENTROID")]
points(departements@data[pos,c("X_CENTROID","Y_CENTROID")],pch=19,col="red")
library(rgeos)
(ctd=gCentroid(Poly_35,byid=TRUE))
points(ctd,pch=19,col="blue")

Comme ça marche mieux, on va utiliser ces centroïdes.

ctd=as.data.frame(gCentroid(departements,byid=TRUE))
plot(departements)
points(ctd,pch=19,col="blue")

Maintenant, il nous faut les résultats des élections, par département. On peut aller scraper le site du ministère d’ intérieur. On a une page par département, donc c’est facile de parcourir. Par contre, dans l’adresse url, il faut le code région. Comme je suis un peu fainéant, au lieu de faire une table de correspondance, on teste tout, jusqu’à ce que ça marche. L’idée est d’aller cherche le nombre de voix obtenues par un des candidats.

candidat="M. Emmanuel MACRON"
library(XML)
voix=function(no=35){
testurl=FALSE
i=1
vect_reg=c("084","027","053","024","094","044","032","028","075","076","052",
"093","011","001","002","003","004")
region=NA
while((testurl==FALSE)&(i<=20)){
reg=vect_reg[i]
nodpt=paste("0",no,sep="")
# if(!is.na(as.numeric(no))){if(as.numeric(no)<10) nodpt=paste("00",no,sep="")}
url=paste("http://elections.interieur.gouv.fr/presidentielle-2017/",reg,"/",nodpt,"/index.html",sep="")
test=try(htmlParse(url),silent =TRUE)
if(!inherits(test, "try-error")){testurl=TRUE
region=reg}
i=i+1
}
tabs tab=tabs[[2]]
nb=tab[tab[,1]==candidat,"Voix"]
a<-unlist(strsplit(as.character(nb)," "))
as.numeric(paste(a,collapse=""))}

On peut alors tester

> voix(35)
[1] 84648

Comme ça semble marcher, on le fait pour tous les départements

liste_dpt=departements@data$CODE_DEPT
nbvoix=Vectorize(voix)(liste_dpt)

On peut alors visualiser sur une carte.

plot(departements,border=NA)
points(ctd,pch=19,col=rgb(1,0,0,.25),cex=nbvoix/50000)

Et on peut aussi tenter pour une autre candidate,

candidat="Mme Marine LE PEN"

et on obtient la carte suivante

plot(departements,border=NA)
points(ctd,pch=19,col=rgb(0,0,1,.25),cex=nbvoix/50000)

Reading text automatically

It is now very easy to read (automatically) some text that can be found in a pdf file. For instance, consider the program of the conference we had yesterday – and today – in Rennes

> library(pdftools)
> scan_pdf <- pdf_text("http://crem.univ-rennes1.fr/Documents/Docs_sem_divers/2017_03_10-11_JJD/JDD_prog.pdf")
> cat(scan_pdf)
Journées Jeunes Docteurs
Programme du jeudi 9 mars 2017
Faculty of Economics - Rennes - Amphi Henri Krier
9h- 9h30 - Accueil
9h30-10h15 :      Présentation du CREM, de la faculté et des activités de recherche liées du ou laboratoire
10h15-10h50 :     Emmanuel LORENZON (Université de Bordeaux, GREThA)
Collusion with a rent seeking agency in sponsored search auctions
10h50-11h25 :     Julien BERTHOUMIEU (Université de Bordeaux, GREThA)
The Impact of “At-the-Border” and “Behind-the-Border” Policies on Cost-Reducing Research
and Development
Co-écrit avec Antoine Bouët

(etc). As you can see, it is working well, even in French, where we have those weird letters (with accents). Here, it is working well because the pdf is vectorized, i.e. it was generated properly, by open office.

But sometimes, we can have only a scanned version of a letter

or just a picture with some typed text. I will not mention hand-writing because it is much more complex.

The other day, my friend Fleur did show me a picture, and some very simple lines of code,

> library('tesseract')
> pic1="https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2017/03/pic1.png"
> text_fr <- ocr(pic1, engine = tesseract("fra"))
> cat(text_fr)
Près de 14.400 décès

Si [épidémie de grippe est un phénomène récurrent. celle de
2016—2017 présente plusieurs spécificités. outre sa virulence :
une survenue plus précoce que d‘habitude. une activité
modérée en médecine ambulatoire. mais un impact fort en
milieu hospitalier.

It looks like we’ve be able to extract typed text from a picture ! I want to check. I have to admit, first of all that installation on a linux machine is tricky: one has to install first leptonica, and then follow some guidelines to install tesseract (see also Artem‘s advices). It took me some time, but I’ve been able to install the package.

The first important step, it to train the algorithm with some texts in French (because it is in French in my picture)

> library('tesseract')
> tesseract_download("fra")

Then, I did try with the picture that Fleur did send me (the picture was inserted in the core of the message)

> pic2="https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2017/03/pic2.png"
> text_fr <- ocr(pic2, engine = tesseract("fra"))
> cat(text_fr)
Près de 14400 decès

s. mm….agw«… ………«…m ……
a……u…u Dhs—ur; ;pmum…. ;: u……
… W»: »… w…q…na… … ……
…… ………u…_ …… mm……
…… nwm/u

… … mm…—mg…»— sa…… su a.…….…
: :mmræwesdæ ; ; m…decnflwtülflws WWW…
un»… M on m…… . … … .. m...… wma:
.,… … V, …… … …;………yg…gn…
…… pe- le…samemeuuœwpwv m…
mum

Clearly, something went wrong here. When I got that output, I thought that I did not train properly the function. But it was not the answer. As described in that post (in French) it is necessary to have a clean picture, to read it properly

And actually, if we zoom in our picture – the first one, used by Fleur, to show me that package – we have

while for the second one – with a lower resolution – we have

It is necessary to have a scan of a typed text with high resolution… And you have to admit that it is awesome….

The good thing is that I have to work with a judge, in France, to assess quality of experts. And since most of the reports are typed, and then scanned, I am glad to have such a function. I just have to make sure that the resolution is high enough…

Evolution du nombre de demandeurs d’emploi

Allez, on continue avec les projets rendus dans le cadre du projet R de la formation en Data Science pour l’Actuariat, Meryem Schalck a proposé une jolie animation des statistiques des demandeurs d’emploi. Sur http://dares.travail-emploi.gouv.fr/, on récupère les séries d’évolution du nombre de demandeurs d’emploi, par département (entre autres), sous forme de fichier Excel. Histoire de simplifier, je mets en ligne une copie en csv,

ufc<- read.csv("http://freakonometrics.free.fr/DemandeV2.csv", header=TRUE, sep = ";")

Pour charger les packages, Meryem propose une jolie petite fonction

Chargepackages <- function(x){
x <- as.character(x)
if(!require(x, character.only = TRUE)){
install.packages(x, dependencies = TRUE)
require(x, character.only = TRUE)
}
}
packages <- list("sp", "rgdal", "animation", "rgeos", "mapproj", "maptools", "dplyr", "ggplot2", "RColorBrewer", "classInt", "PBSmapping", "ggmap", "splancs", "osmar", "maps", "scales", "osmar", "geosphere", "RCurl", "bitops", "XML", "splancs", "PBSmapping", "htmltools"
)
lapply(packages, Chargepackages)

On va retravailler un peu les données pour commencer

x<-ufc$nombre
y<-ufc$demande
dt<-as.Date(ufc$date,"%d/%m/%Y")
n<-length(x)

et ensuite on fait notre animation,

colfunc <- colorRampPalette(c("green","yellow", "red"))
YlOrBr <- data.frame(col=colfunc(length(unique(y))),num=unique(y)[order(unique(y))])
ufc %>%
left_join(YlOrBr, by=c("demande"="num")) -> YlOrBr
oopt = ani.options(interval = 0.5)
library(scales)
YlOrBr$Date<-as.Date(YlOrBr$date,"%d/%m/%Y")
ani.record(reset = TRUE) # clear history before recording
for (i in 1:length(unique(dt))) {
YlOrBr[i>=x,] %>%
ggplot() +  geom_point(aes(x=Date,y=demande),colour=as.character(YlOrBr$col[i>=x]))+
scale_x_date(breaks=date_breaks("6 month"),limits = c(min(dt),max(dt))) +
ylim(c(1000,max(y))) +
ggtitle ("Nombre de demandeurs d'emploi par mois entre 01/2000 et 04/2016") +
xlab("Date") +
ylab("Nombre de demandeurs d'emploi cat A ensemble") +
theme(axis.text.x  = element_text(angle=45, vjust=0.5, size=5)) -> p
print(p)
ani.record()
}
oopts = ani.options(interval = 0.2)
ani.record()
saveHTML(ani.replay(),img.name = "record_plot")

et voilà

Localisation des accidents corporels

Encore un billet dans le cadre du projet R de la formation en Data Science pour l’Actuariat, avec quelques lignes de codes pour visualiser les accidents corporels, en France. Il y aura d’autres billets sur ces données dans les jours à venir. Le premier s’inspire du projet de Guillaume Gerber. L’idée est de normaliser par la population d’un département, quand on va visualiser le nombre d’accidents.

On va commencer par récupérer notre fond de carte

library(rgdal)
library(sp)
library(data.table)
library(dplyr)
library(plyr)
download.file(url = "http://professionnels.ign.fr/sites/default/files/GEOFLA_1-1_SHP_LAMB93_FR-ED111.tar.gz",
destfile="GEOFLA.tar.gz")
untar("GEOFLA.tar.gz")

Ces données contiennent la population. Mais par commune. Donc on va aggréger pour l’avoir par département

com <- readOGR(dsn = "./GEOFLA_1-1_SHP_LAMB93_FR-ED111/COMMUNES/", layer="COMMUNE", stringsAsFactors=FALSE,encoding = "UTF-8")
com <- readOGR(dsn = "/home/arthur/Téléchargements/projet-R-DSA-2016/GEOFLA_1-1_SHP_LAMB93_FR-ED111/COMMUNES/", layer="COMMUNE")
dep <- readOGR(dsn = "/home/arthur/Téléchargements/projet-R-DSA-2016/GEOFLA_1-1_SHP_LAMB93_FR-ED111/DEPARTEMENTS/", layer="DEPARTEMENT", stringsAsFactors=FALSE,encoding = "UTF-8")
pop <- as.data.table(tapply(X=com$POPULATION, INDEX = com$CODE_DEPT, FUN = sum), keep.rownames=TRUE)
colnames(pop) <- c("CODE_DEPT","POPULATION")
superficie <- as.data.table(tapply(X=com$SUPERFICIE, INDEX = com$CODE_DEPT, FUN = sum), keep.rownames=TRUE)
colnames(superficie) <- c("CODE_DEPT","SUPERFICIE")
dep@data <- inner_join(dep@data, pop)
dep@data <- inner_join(dep@data, superficie)
dep@data$POPULATION <- dep@data$POPULATION * 1000

On va ensuite récupérer les données d’accidents de la route

debut_url_base_accident = "https://www.data.gouv.fr/s/resources/"
acc_caract <- read.csv(file = paste(debut_url_base_accident,"base-de-donnees-accidents-corporels-de-la-circulation-sur-6-annees/20150806-155035/caracteristiques_2010.csv",sep=''),colClasses=c("com"="character","dep"="character"))
acc_caract <- rbind(acc_caract,read.csv(file = paste(debut_url_base_accident,"base-de-donnees-accidents-corporels-de-la-circulation-sur-6-annees/20150806-154723/caracteristiques_2011.csv",sep=""), colClasses=c("com"="character","dep"="character")))
acc_caract <- rbind(acc_caract,read.csv(file = paste(debut_url_base_accident,"base-de-donnees-accidents-corporels-de-la-circulation-sur-6-annees/20150806-154431/caracteristiques_2012.csv",sep=""), colClasses=c("com"="character","dep"="character")))
acc_caract <- rbind(acc_caract,read.csv(file = paste(debut_url_base_accident,"base-de-donnees-accidents-corporels-de-la-circulation-sur-6-annees/20150806-154105/caracteristiques_2013.csv",sep=""), colClasses=c("com"="character","dep"="character")))
acc_caract <- rbind(acc_caract,read.csv(file = paste(debut_url_base_accident,"base-de-donnees-accidents-corporels-de-la-circulation-sur-6-annees/20150806-153701/caracteristiques_2014.csv",sep=""), colClasses=c("com"="character","dep"="character")))
acc_caract <- rbind(acc_caract,read.csv(file = paste(debut_url_base_accident,"base-de-donnees-accidents-corporels-de-la-circulation/20160909-181230/caracteristiques_2015.csv",sep=""), colClasses=c("com"="character","dep"="character")))
acc_caract$dep[which(acc_caract$dep %in% "201")] <- "2A0"
acc_caract$dep[which(acc_caract$dep %in% "202")] <- "2B0"
acc_caract$dep <- substr(acc_caract$dep, 1, 2)

Maintenant, on peut compter, par année, par département (ou en aggrégeant, temporellement)

dep_with_nb_acc <- function(acc_caract, dep,nb_an=1,normalize=FALSE){
acc_nb_par_dep <- count(acc_caract,"dep")
acc_nb_par_dep$dep <- substr(acc_nb_par_dep$dep, 1, 2)
names(acc_nb_par_dep)[names(acc_nb_par_dep) == "dep"] <- "CODE_DEPT"
dep@data <- inner_join(dep@data, acc_nb_par_dep)
dep@data$SUPERFICIE <- dep@data$SUPERFICIE/100 # en km^2
dep@data$freq_par_hab <- (dep@data$freq/nb_an)/dep@data$POPULATION
m <- sum(dep@data$freq/nb_an)/sum(dep@data$POPULATION)
if(normalize) dep@data$freq_par_hab <- log(dep@data$freq_par_hab/m)
return(dep)
}

On va alors constituer nos bases,

data_plot <-  c(
"2010_2015" = dep_with_nb_acc(acc_caract, dep,nb_an = 6),
"2010_2015_n" = dep_with_nb_acc(acc_caract, dep,nb_an = 6,normalize=TRUE))

La première correspond juste au nombre d’accident, sur 6 ans, normalisé par la population (ce qui pourraît être vu comme une fréquence d’accident corporel)

zmax = max(data_plot[[1]]@data$freq_par_hab)
spplot(obj = data_plot$'2010_2015',"freq_par_hab",at = seq(0, zmax, by = zmax/10),main = "")

mais on peut aussi normaliser par la fréquence nationale, afin de faire ressortir les départements les plus dangereux. On est également passé sur une échelle logarithmique,

zmin = min(data_plot[[8]]@data$freq_par_hab)
zmax = max(data_plot[[8]]@data$freq_par_hab)
spplot(obj = data_plot$'2010_2015_n',"freq_par_hab",at = seq(zmin, zmax, by = (zmax-zmin)/10),main = "")

Visualiser l’évolution de la taxe d’habitation

Dans le cadre du projet R de la formation en Data Science pour l’Actuariat, on continue à explorer les données sur data.gouv, par exemple https://data.gouv.fr/fr/datasets/impots-locaux/, sur les impôts locaux. C’est ce que Nicolas Jaudel avait proposé,

library(sp)
tax_hab=read.csv("http://freakonometrics.free.fr/F7815_taxe_hab.csv", col.names=c("CODDPT", "LIBDPT","E2001", "E2002","E2003","E2004","E2005", "E2006","E2007","E2008","E2009", "E2010"), header=FALSE, dec=",", sep=";", stringsAsFactors = FALSE)

On va extraire une information intéressante, par exemple la variation du taux moyen pour la taxe d’habitation, sur 10 ans,

tax_hab$taux_evo=tax_hab$E2010/tax_hab$E2001-1

Pour la carte des départements, on peut aller récupérer un fond de carte sur biogeo.ucdavis.edu.

download.file("http://freakonometrics.free.fr/FRA_adm2.rds","FRA_adm2.rds")
FR=readRDS("FRA_adm2.rds")
donnees_carte=data.frame(FR@data)
donnees_carte=cbind(donnees_carte,"taux_evo"=0)
for (i in 1:96){
donnees_carte[i,"taux_evo"]=tax_hab[tax_hab$CODDPT==donnees_carte$CCA_2[i],'taux_evo']*100
}

On a ainsi toutes les données, et on peut faire une carte

library(cartography)
cols <- rev(carto.pal(pal1 = "red.pal",n1 = 10,pal2="green.pal",n2=10))
plot(FR, col = "grey", border = "gray1", bg="#A6CAE0",xlim=c(-5.2,12))
choroLayer(spdf = FR,
df = donnees_carte,
var = "taux_evo",
breaks = seq(-30,70,5),
col = cols,
legend.pos = "topright",
legend.title.txt = "Evolution",
legend.values.rnd = 2,
add = TRUE)

Le sport, une activité de riches ?

Allez, un petit billet léger pour bien commencer les vacances de Noël. Toujours sur le projet R de la formation en Data Science pour l’Actuariat, Cyril Legrand a proposé de fusionner deux bases, une sur le revenu moyen par foyer (fiscal) par commune, et une sur le nom de licenciés dans les clubs sportifs. Il est faire un peu de retraitement des codes INSEE pour la jointure, car par exemple Marseille est codée sur les arrondissements dans une des bases, sur la ville dans l’autre.

Continue reading Le sport, une activité de riches ?

A quelle distance d’une banque habite-t-on ?

Dans le cadre du projet de R de la formation en Data Science pour l’Actuariat, je vais continuer à mettre en ligne des morceaux de codes qui peuvent être utiles, dans un contexte spatial. Le dernier billet, sur cartographier le vote pour le Brexit, avait été repris (et bien amélioré) sur le site des voisins de rgeomatic. Aujourd’hui, je vais m’inspirer du travail d’Etienne Flichy qui mixe répartition de la population sur le territoire, et localisation des agences bancaires.

On parle des banques ici, mais si on a une base avec les coiffeurs, les boulangeries, etc, on peut faire la même chose ! (autant dire qu’on va pouvoir s’amuser quand la base sirene sera rendue ouverte – dans les semaines à venir). On va supposer que l’on a une base avec toutes les banques géocodées. Bon, pour l’exercice, on va utiliser la localisation des agences bancaires, en utilisant les données de cbanque.com. C’est assez facile d’aller scraper le site, quand on regarde la façon dont sont faites les pages, e.g. http://cbanque.com/pratique/agences/credit-cooperatif/35/. Là on récupère les adresses (postales) et on peut utiliser https://adresse.data.gouv.fr/csv/ (ou différents outils) pour géolocaliser les adresses.

Continue reading A quelle distance d’une banque habite-t-on ?

Cartographier le vote pour le Brexit

Je suis, ces temps-ci plongé dans les projets de R que j’avais donnés pour la formation Data Science pour l’Actuariat. Comme j’ai eu plein de choses intéressantes, je me suis dit que je pourrais tenter de reprendre dans des billets sur le blog (en plus, ça me permet de vérifier au passage que les codes tournent).

Pour commencer la série (oui, il y en aura d’autres), Flavien Thery a proposé une cartographie du vote pour le Brexit. La première étape était de récupérer un fond de carte (on va utiliser la seconde zone administrative ici)

library(sp)
library(raster)
download.file("http://biogeo.ucdavis.edu/data/gadm2.8/rds/GBR_adm2.rds","GBR_adm2.rds")
UK=readRDS("GBR_adm2.rds")
UK@data[159,"HASC_2"]="GB.NR"
plot(UK, xlim = c(-4,-2), ylim = c(50, 59), main="UK areas")

Cette carte est un peu surprenante… parce qu’on est sur le Royaume uni, ce qui inclu l’Angleterre, le Pays de Galles, l’Ecosse et l’Irlande du Nord. Et l’Irlande du Sud n’est pas sur la carte… Autant la rajouter (ça sera plus joli si on met du bleu pour l’eau)

download.file("http://biogeo.ucdavis.edu/data/gadm2.8/rds/IRL_adm0.rds","IRL_adm0.rds")
IRL=readRDS("IRL_adm0.rds")
plot(IRL,add=TRUE)

Mais on commence à rajouter l’Irlande, pourquoi pas la France, qui est juste en bas à droite, et qu’on devrait voir un peu…

download.file("http://biogeo.ucdavis.edu/data/gadm2.8/rds/FRA_adm0.rds","FRA_adm0.rds")
FR=readRDS("FRA_adm0.rds")
plot(FR,add=TRUE)

Bon, arrêtons là….

On peut ensuite récupérer les données du référendum sur le Brexit (que je stocke sur le blog histoire de gagner un peu de temps)

loc="https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2016/12/EU-referendum-result-data.csv"
referendum=read.csv(loc,header=TRUE,dec=".",sep=",",stringsAsFactors = FALSE)
referendum=referendum[c(3,6,13,14)]
library(plyr)
referendum=ddply(referendum,.(Region,HASC_code),summarise,Remain=sum(Remain),Leave=sum(Leave))

On peut retrouver que le Brexit a gagné, avec 51,89% des votes exprimés (ce qui est conforme à ce que dit wikipedia)

> sum(referendum$Leave)/(sum(referendum$Leave)+sum(referendum$Remain))
[1] 0.5189184

On peut regarder, région par région, si c’est le leave ou le remain qui a gagné, en utilisant

referendum=referendum[c(referendum$Region!="Northern Ireland"),]
referendum=referendum[c(referendum$HASC_code!="Gibraltar"),]
row.names(referendum)=seq(1,nrow(referendum),1)
leave_or_remain=cbind(referendum,"Brexit?"=0)
leave_or_remain[,"Brexit?"]=ifelse(leave_or_remain$Remain<leave_or_remain$Leave,rgb(1,0,0,.7),rgb(0,0,1,.7))
map_data=data.frame(UK@data)
map_data=cbind(map_data,"Brexit"=0)
for (i in 1:nrow(map_data)){
if(map_data[i,"NAME_1"]=="Northern Ireland"){
map_data[i,"Brexit"]="blue"}else{
map_data[i,"Brexit"]=as.character(leave_or_remain[leave_or_remain$HASC_code==map_data$HASC_2[i],'Brexit?'])
}
}
plot(UK, col = map_data$Brexit, border = "gray1", xlim = c(-4,-2), ylim = c(50, 59), main="How the UK has voted?", bg="#A6CAE0")
plot(IRL, col = "lightgrey", border = "gray1",add=TRUE)
plot(FR, col = "lightgrey", border = "gray1",add=TRUE)

legend(-1,59,c("Leave","Remain"),fill=c("red","blue"),bty="n")

(on rajoute une petite légende pour être plus clair). Mais on peut aller plus loin en représentant le pourcentage obtenu, par région, par le camp du leave. Pour ça, on peut utiliser le package cartography

library(cartography)
cols <- carto.pal(pal1 = "red.pal",n1 = 5,pal2="green.pal",n2=5)
map_data=data.frame(UK@data)
map_data=cbind(map_data,"Percentage_Remain"=0)
for (i in 1:nrow(map_data)){
if(map_data[i,"NAME_1"]=="Northern Ireland"){
map_data[i,"Percentage_Remain"]=55.78}else{         map_data[i,"Percentage_Remain"]=100*leave_or_remain[leave_or_remain$HASC_code==map_data$HASC_2[i],"Remain"]/(leave_or_remain[leave_or_remain$HASC_code==map_data$HASC_2[i],"Remain"]+leave_or_remain[leave_or_remain$HASC_code==map_data$HASC_2[i],"Leave"])
}
}
plot(UK, col = "grey", border = "gray1", xlim = c(-4,-2), ylim = c(50, 59),bg="#A6CAE0")
plot(IRL, col = "lightgrey", border = "gray1",add=TRUE)
plot(FR, col = "lightgrey", border = "gray1",add=TRUE)
choroLayer(spdf = UK,
df = map_data,
var = "Percentage_Remain",
breaks = seq(0,100,10),
col = cols,
legend.pos = "topright",
legend.title.txt = "",
legend.values.rnd = 2,
add = TRUE)

(là encore, une légende permet d’aider à lire la carte). Amusant, n’est ce pas ?

R in Insurance, 2017

Following the successfull conferences in London (2013, 2014, 2016) and in Amsterdam (2015), the next edition will take place in Paris. The R in insurance 2017 will take place in ENSAE on June 8.

This one-day conference will focus again on applications in insurance and actuarial science that use R, the lingua franca for statistical computation. The intended audience of the conference includes both academics and practitioners who are active or interested in the applications of R in insurance. The two invited speakers are Katrien Antonio (KU Leuven) and Julie Seguela (Covea). It will be a nice event !

La centralisation française, et la répartition de la population sur le territoire

Suite à mon billet d’hier soir sur la distribution en France, j’ai eu plusieurs commentaires – sur Twitter – qui notaient que ce n’était pas surprenant que la France soit aussi concentrée, en terme de population, compte tenu de l’importance de la centralisation en France (par opposition à l’Allemagne, par exemple).  Et comme le note wikipedia, à la page sur le centralisme, “depuis la Révolution française (et même avant), la France est un État très centraliste“.

Il y a quelques années, Mattia Bunel (aka @mattiabunel) avait rédigé un joli mémoire sur l’influence des contraintes environnementales sur la distribution de la population française entre 1793 et 1999. Au delà du travail de rédaction, il y a surtout un gros travail sur les données. Mattia a remis en forme les données de cassini.ehess.fr. Et Mattia a eu la gentille de m’envoyer ses données, avec en ligne les villages en France, leur superficie, et leur population, à plusieurs dates, entre 1793 et 1999. Pour récupérer la superficie (Mattia a passer du temps à mettre à jour, pour les villages qui ont fusionné, en particulier) le code est tout simplement

> base=read.csv("/Cassini/export.csv")
> dim(base)
[1] 41409    40
> S=base$Superficie
> n=nchar(as.character(S))
> S=substr(S,1,n-3)
> Surface=as.numeric(gsub(" ", "", as.character(S), fixed = TRUE))
> idx=which(!is.na(Surface))
> B=base[idx,]
> dim(B)
[1] 36576    40

On a ainsi les 36,000 communes françaises. Ensuite, on peut extraire, par année, la population. Histoire d’avoir un calcul de l’indice de Gini qui soit cohérent avec ce que j’avais fait hier – dans Non-Uniform Population Density in some European Countries – l’idée est la suivante : supposons qu’il y ait 2 villages, un de superficie 4, et de population 4, et un autre de superficie 1, et de population 3. Dans le billet d’hier, je raisonnais pas unité spatiale (en l’occurrence un petit carré)

et on va faire pareil ici. Autrement dit, on va répartir la population uniformément au sein du village de taille 4, ce qui va faire 4 “villages unitaires” avec 1 habitant chacun,

On a ainsi 5 unités, 4 avec une population de 1, et 1 avec une population de 3. De cet échantillon – {1,1,1,1,3} – on peut calculer l’indice de Gini. La fonction pour extraire la population, et calculer l’indice de Gini est

> library(ineq)
> LC=function(annee){
+   nom=paste("X",annee,sep="")
+   P=B[,nom]
+   P=gsub(" ", "", as.character(P), fixed = TRUE)
+   Pop=as.numeric(substr(P,3,nchar(P)))
+   P=Pop[which(!is.na(Pop))]
+   S=Surface[which(!is.na(Pop))]
+   D=P/S
+   S1=round(S/20)
+   X=rep(NA,sum(S1))
+   s=0
+   for(i in 1:length(S1)){
+     X[s+1:S1[i]]=D[i]
+     s=s+S1[i]
+   }
+   Gini(X)
+ }

Si on calcule notre indice de Gini pour toutes les années, on obtient

> Y=names(base)[7:ncol(base)]
> Y=as.numeric(substr(Y,2,5))
> gini=Vectorize(LC)(Y)
> plot(Y,gini,type="b")

On retrouve un indice de Gini qui dépasse 0.7 aujourd’hui (que nous avions obtenu hier, avec une méthodologie assez différente), mais surtout, on voit que l’indice de Gini n’a cessé d’augmenter, depuis la révolution française… L’interprétation rapide serait que le centralisme français ne cesse pas d’augmenter, depuis 200 ans.