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.

Ensuite, il nous faut savoir où habitent les gens. On peut utiliser les communes, mais ce n’est pas assez précis. Plus intéressant, on peut utiliser les zones iris. On récupère les informations relatives à la population en utilisant

mondir="/home/arthur/"
for (x in c("a","b","c","d","e")){
urlf<-paste("http://telechargement.insee.fr/fichiersdetail/RP2012/txt/RP2012_LOGEMTz",x,"_txt.zip",sep="")
zipf<-paste(mondir,"/",x,".zip",sep="")
download.file(urlf,zipf)
unzip(zipf,files=unzip(zipf,list=TRUE)$Name[1],overwrite = TRUE)
file.remove(zipf)
}

et on récupère les contours pour les fonds de carte sur

url_map<-"http://www.data.gouv.fr/s/resources/contour-des-iris-insee-tout-en-un/20150428-161348/iris-2013-01-01.zip"
mapf<-paste(getwd(),"/map.zip",sep="")
download.file(url_map,mapf)
unzip(mapf,overwrite = TRUE)

On remet un peu en forme, pour avoir des belles données

vectvar<-c("COMMUNE","IRIS","IPONDL","INPER")
newdata<-data.frame(row.names = vectvar)
for (x in c("a","b","c","d","e")){
data<-read.csv(paste("FD_LOGEMTZ",toupper(x),"_2012.txt",sep=""),header=TRUE,sep=";",dec=".")
newdata<-rbind(newdata,data[,vectvar])
print(x) # permet de visualiser où on en est du réassemblage
}
write.table(newdata,file="tablepop.txt",sep=";",dec=".",col.names=TRUE,row.names = FALSE)

Il faut corriger certaines dénominations, dans la base

newdata$DCOMIRIS<-as.character(newdata$IRIS)
newdata[newdata$IRIS=="ZZZZZZZZZ","DCOMIRIS"]<-paste(
as.character(newdata[newdata$IRIS=="ZZZZZZZZZ","COMMUNE"]),"0000",sep="")
newdata$NPER<-as.integer(as.character(newdata$INPER))

On a maintenant des NA dans la population, et on met des 0 à la place

newdata$NPER[is.na(newdata$NPER)]<-0
pop<-aggregate(
newdata[newdata$INPER!="Y","IPONDL"]*newdata[newdata$INPER!="Y","NPER"],
by=list(DCOMIRIS=newdata[newdata$INPER!="Y","DCOMIRIS"]), FUN=sum)
names(pop)[2]<-"NH"

Cette fois, on doit être bon. On vérifie la population, dans notre base de données

> sum(pop$NH)
[1] 63759971

Ce qui correspond à la population totale. On n’a perdu personne (pour l’instant…). On peut visualiser la densité de population, par zone iris, et positionner les centroïdes des zones iris (pour les calculs des distances, on placera toute la population au centre)

library(rgdal)
France<-readOGR(dsn="./.",layer="iris-2013-01-01")
France<-subset(France,!duplicated(France@data$DCOMIRIS)) # dédoublonne la carte
levels(France@data$NOM_COM)<-c(levels(France@data$NOM_COM),"inconnu")
France@data$NOM_COM[is.na(France@data$NOM_COM)]<-"inconnu"

On indexe ensuite toutes les zones iris,

library(rgeos)
France@data$DEP<-substr(as.character(France@data$DEPCOM),1,2)
France@data$NUM<-row.names(France@data)
cent<-gCentroid(France,byid=TRUE) ## récupère les centroides
irisdata<-cbind(France@data,cent@coords)
names(irisdata)[10:11]<-c("x_centroid","y_centroid")
irisdata$SURF<-as.vector(sapply(slot(France,"polygons"),slot,"area"))
irisdata$DCOMIRIS<-as.character(irisdata$DCOMIRIS)

et on rajoute la population

irisdata<-merge(irisdata,pop,by=c("DCOMIRIS"),all.x=TRUE)
irisdata[,10:13][is.na(irisdata[,10:13])]<-0

On peut vérifier la population totale qu’il nous reste,

> sum(irisdata$NH)
[1] 62238867

on a perdu un bon million de personnes… mais on va faire avec (pour l’instant). Pour visualiser notre population, et nos banques, on va se mettre sur une région, par exemple Bourgogne-Franche Comté (2,7 millions d’habitants), sinon c’est vraiment trop gros. En plus, la base de localisation des agences bancaires se limite à cette région (soit un peu plus de 3000 agences… j’ai mis une extraction dans la base en ligne, pour ceux qui veulent jouer avec le code)

France<-merge(France,irisdata[,9:13],by=c("NUM"),all.x=TRUE)
library(classInt)
library(RColorBrewer)
France@data$DENSITE<-France@data$NH/France@data$SURF
seuilcoul<-classIntervals(France@data$DENSITE,5,style="quantile",dataPrecision = 1)
France@data$COULEUR<-findColours(seuilcoul,brewer.pal(5,"Blues"))
zone=c("21","25","39","58","70","71","89","90")
library(sp)
agences<-read.csv("http://freakonometrics.hypotheses.org/files/2016/12/geocodage_agences.csv",header=TRUE,sep=";",dec=".")
agc<-SpatialPoints(agences[,11:12])
agc@proj4string@projargs<-'+proj=longlat +datum=WGS84'
agc<-spTransform(agc,France@proj4string)

Pour avoir quelque chose de joli, on peut faire du leaflet (au début du billet, mais également en ligne, en grand format)

m<-leaflet() %>%
addTiles() %>%
addPolygons(data=France[France@data$DEP%in%zone,],weight=1,
fillColor=France[France@data$DEP%in%zone,]@data$COULEUR,
opacity=1,fillOpacity=0.6,color=grey(0.5)) %>%
addCircleMarkers(lng=agc$LON,lat=agc$LAT,radius=1,color="red")
m

ou si on zoom un peu

Ensuite, pour chaque centroïde de zone iris (où on concentre la population), on peut calculer la distance à la banque la plus proche.

library(spatstat)
irisbfc=irisdata[irisdata$DEP %in% zone,]
agences$RESEAU=as.character(agences$RESEAU)
centroid=SpatialPoints(irisbfc[,10:11])
centroid@proj4string@projargs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
matdist=spDists(centroid,agc,longlat=TRUE) # calcule toutes les distances iris / agences
dst=apply(matdist,1,min)
pop=irisbfc$NH
prop=function(d) sum(pop[dst<=d])/sum(pop)
u=seq(0,10,by=.1)
v=Vectorize(prop)(u)
plot(u,v*100)


Autrement dit, ici, 80% des gens habitent à moins de 3km (à vol d’oiseau) d’une agence bancaire. Et seulement 1% sont à plus de 6,5km, en Bourgogne-France Comté en tous cas…. Bien sûr, on pourrait faire mieux, en utilisant les distances réelles (mais on en parlera bientôt…)

---

---