Cartographie des écoles de Bruxelles avec R

Bien que R puisse même faire le café, nous nous contenterons aujourd'hui de produire une carte de Bruxelles et d'y placer des écoles.

Préparation de l'espace de travail

# Réinitialiser Espace de travail
rm(list=ls())

# Définir le dossier de travail
root.dir="C:/Users/.../CartoR"
setwd(root.dir)

# Activer les libraries
library(maptools)
library(foreign)
library(lattice)
library(latticeExtra)

Pensez à définir le chemin menant à votre propre espace de travail et à télécharger les packages nécessaires pour la production d'une carte.

Importation des données écoles

ecoles=read.csv(paste(root.dir,"Ecoles.txt", sep="/"),
 header=T, sep=",", dec = ".")

Le fichier Ecole.txt est disponible ici. Téléchargez le fichier et décompressez-le dans votre espace de travail. Il s'agit d'un échantillon de 100 écoles bruxelloises. Si les coordonnées sont réelles, les données ont été simulées.

Une carte n'est rien de plus qu'un repére cartésien

xyplot(Y~X, ecoles,
 aspect="iso",
 panel = function(x, y, ...) {
  panel.xyplot(x, y,
   pch=21,
   ...)
})

Nous générons un graphique où chaque point est placé selon ses coordonnées XY. Nous spécifions l'option "iso" pour l'argument aspect afin d'obtenir un repére orthonormé. Notons, à ce stade, que la fonction xyplot appelle en fait la fonction panel.xyplot.

carte 1

Suppression du cadre

theme=trellis.par.get()
theme$layout.widths[c("axis.key.padding",
 "right.padding", "left.padding", "ylab.right", "axis.right"
 )]=0
theme$layout.widths["right.padding"]=.5
theme$layout.widths["key.right"]=1.5
theme$layout.heights[c("key.axis.padding", "main.key.padding",  "top.padding", "bottom.padding","ylab.top", "axis.top"
 )]=0
theme$axis.line$alpha=0
xyplot(Y~X, ecoles,
 aspect="iso", par.settings=theme,
 panel = function(x, y, ...) {
  panel.xyplot(x, y,
   pch=21,
  ...)
 },
 scales=list(x=list(at=NULL), y=list(at=NULL)),
 xlab="", ylab=""
)

Il est possible d'appliquer des thémes aux graphiques lattice (donc également la fonction xyplot). Trellis.par.get() permet de récupérer le théme par défaut et de l'attribuer à l'objet theme. Ceci fait, on peut le modifier et l'appliquer à notre graphique grace à l'option par.settings. Ici, nous nous contentons de retirer le cadre. Nous finissons cette action avec les arguments scales, xlab, ylab.

carte 2

Grosseur des points proportionnelle à la taille de l'école

Ajust=3.5
ecoles[,"Point.cex"]=ecoles[,"Taille"]/(Ajust*min(ecoles[,"Taille"]))
xyplot(Y~X, ecoles,
 aspect="iso", par.settings=theme,
 panel = function(x, y, ...) {
  panel.xyplot(x, y,
   pch=21,
   cex=ecoles[,"Point.cex"],
  ...)
 },
 scales=list(x=list(at=NULL), y=list(at=NULL)),
 xlab="", ylab=""
)

Nous ajoutons une colonne Point.cex à la data frame ecoles. Celle-ci nous permettra de définir la taille de chaque point en précisant l'argument cex de la fonction panel.xyplot. La taille des points est calculée automatiquement sur base de la taille minimum afin qu'ils soient tous visibles. Si les points ne sont pas assez gros, nous pouvons augmenter la valeur attribuée à l'objet Ajust.

carte 3

Ajouter une légende

xyplot(Y~X, ecoles,
 aspect="iso", par.settings=theme,
 panel = function(x, y, ...) {
  panel.xyplot(x, y,
   pch=21,
   cex=ecoles[,"Point.cex"],
  ...)
 },
 scales=list(x=list(at=NULL), y=list(at=NULL)),
 xlab="", ylab="",
 legend=list(
  inside=list(fun=draw.key(key=list(title="Nombre\nd'éléves",
   cex=.6,
   points=list(pch=21,
    cex=c(200, 400, 600, 800)/(Ajust*min(ecoles[,"Taille"])),
    col=rgb(51,51,51, maxColorValue=255),
    fill=rgb(255,255,255, maxColorValue=255)),
   text=list(c("200", "400", "600", "800"),
    padding.text = 10,
    height = 10,
    cex=.8),
   padding.text=4
  )),
  x = .01, y = .3)
 )
)

carte 4

Couleur des points relative à l'ISE

nbClass=6
class.pal=c("#005A32", "#41AB5D", "#A1D99B", "#9ECAE1", "#4292C6", "#084594")
class.pal=data.frame(Class=1:nbClass,
 Point.col=class.pal,
 Point.fill=paste(class.pal, "B3", sep=""))

Nous spécifions le nombre de couleurs différentes que nous voulons. Nous choisissons 6 couleurs (palette) pour la lisibilité et les précisons au format hexadécimal. Pour chaque couleur Point.col, nous définissons une version transparente Point.fill.

class.break=floor(quantile(1:sum(ecoles[,"Taille"]), seq(0,1, by=1/nbClass), names=F))
tmp.a=cumsum(ecoles[order(ecoles[,"ISE"]),"Taille"])
tmp.b=c(sapply(class.break[1:nbClass], function (x) which(tmp.a==min(tmp.a[(tmp.a>x)]))), length(tmp.a))
ecoles[order(ecoles[,"ISE"]),"Class"]=cut(tmp.a, tmp.a[tmp.b], labels=F, include.lowest=T)

Nous voulons à présent diviser nos écoles selon leur ISE (indice socio-économique) afin qu'elles appartiennent à l'une des 6 catégories de couleur. Différentes possibilités existent pour catégoriser nos écoles. Ici, nous faisons en sorte que chaque classe contienne un nombre d'étudiants équivalent à la maniére des classes ISE de l'encadrement différencié.

ecoles=merge(ecoles, class.pal, by="Class")
class.lab=do.call("rbind",by(ecoles, ecoles[,"Class"], simplify = F, function(x) {
 paste(round(min(x["ISE"]), 2), round(max(x["ISE"]), 2), sep=";")
}))

Ceci fait, nous attribuons à chaque école sa couleur et créons les labels pour la légende.

xyplot(Y~X, ecoles,
 aspect="iso", par.settings=theme,
 panel = function(x, y, ...) {
  panel.xyplot(x, y,
   pch=21,
   cex=ecoles[,"Point.cex"],
   col=as.character(ecoles[,"Point.col"]),
   fill=as.character(ecoles[,"Point.fill"]),
  ...)
 },
 scales=list(x=list(at=NULL), y=list(at=NULL)),
 xlab="", ylab="",
 legend=list(
  inside=list(fun=draw.key(key=list(title="Nombre\nd'éléves",
   cex=.6,
   points=list(pch=21,
    cex=c(200, 400, 600, 800)/(Ajust*min(ecoles[,"Taille"])),
    col=rgb(51,51,51, maxColorValue=255),
    fill=rgb(255,255,255, maxColorValue=255)),
   text=list(c("200", "400", "600", "800"),
    padding.text = 10,
    height = 10,
    cex=.8),
   padding.text=4
   )),
   x = .01, y = .3),
  inside=list(fun=draw.key(key=list(title="ISE\n(Ecoles)",
   cex=.6,
   points=list(pch=21,
    cex=2,
    col=as.character(class.pal[,"Point.col"]),
    fill=as.character(class.pal[,"Point.fill"])),
   text=list(class.lab, cex=.8),
   padding.text=3
   )),
   x = .01, y = .99)
))

carte 5

Ajouter fond de carte

bruxelles = readShapeSpatial("UrbAdm_Sd.shp")

La fonction readShapeSpatial permet de lire les fichiers Shape. Nous avons préalablement téléchargé les données UrbIS sur le site du CIRB.

bruxelles@data[,"NbLigne"]=as.numeric(row.names(bruxelles@data))

A l'intérieur de ce fichier shape se cache une data frame reprenant les attributs des différents secteurs (provenant du fichier "UrbAdm_Sd.dbf"). Comme il est important de conserver l'ordre initial de cette data frame, nous récupérons les numéros de ligne dans une nouvelle colonne.

xyplot(Y~X, ecoles,
 aspect="iso", par.settings=theme,
 panel = function(x, y, ...) {
  panel.polygonsplot(x=coordinates(bruxelles)[,1],
   y=coordinates(bruxelles)[,2],
   z=bruxelles@data[,"NbLigne"],
   subscripts=1:nrow(bruxelles@data),
   grid.polygons=bruxelles,
   sp.layout=NULL,
   col="#888888",
   col.regions="#FFFFFF",
   ...)
  panel.xyplot(x, y,
   pch=21,
   cex=ecoles[,"Point.cex"],
   col=as.character(ecoles[,"Point.col"]),
   fill=as.character(ecoles[,"Point.fill"]),
  ...)
 },
 scales=list(x=list(at=NULL), y=list(at=NULL)),
 xlab="", ylab="",
 xlim=c(141100, 158000), ylim=c(161400, 178300),
 legend=list(
  inside=list(fun=draw.key(key=list(title="Nombre\nd'éléves",
   cex=.6,
   points=list(pch=21,
    cex=c(200, 400, 600, 800)/(Ajust*min(ecoles[,"Taille"])),
    col=rgb(51,51,51, maxColorValue=255),
    fill=rgb(255,255,255, maxColorValue=255)),
   text=list(c("200", "400", "600", "800"),
    padding.text = 10,
    height = 10,
    cex=.8),
   padding.text=4
   )),
   x = .01, y = .3),
  inside=list(fun=draw.key(key=list(title="ISE\n(Ecoles)",
   cex=.6,
   points=list(pch=21,
    cex=2,
    col=as.character(class.pal[,"Point.col"]),
    fill=as.character(class.pal[,"Point.fill"])),
   text=list(class.lab, cex=.8),
   padding.text=3
   )),
   x = .01, y = .99)
))

C'est ici que lattice montre toute sa flexibilité. Il est possible de récupérer les fonctions des autres types de graphiques pour les utiliser dans notre xyplot. La fonction appelée par spplot (méthode pour représenter les données spatiales) est panel.polygonsplot. Nous définissons également les limites inférieures et supérieures des axes pour que l'entiéreté de la carte apparaisse.

carte 6

Utiliser les données secteur pour colorier le fond de carte

secteurs=read.csv(paste(root.dir,"Secteurs.txt", sep="/"),  header=T, sep=",", dec = ".")

Nous importons les données qui nous serviront à colorier le fond de carte. Il s'agit de la densité de population des 5 à 19 ans par secteurs.

bruxelles.com = read.dbf("UrbAdm_Mu.dbf", as.is = FALSE)
for (i in bruxelles.com[,"ID"]) {
 bruxelles@data[!is.na(bruxelles@data[,"MU_ID"])&(bruxelles@data[,"MU_ID"]==i), "MU_ID"] =
  as.character(bruxelles.com[(bruxelles.com[,"ID"]==i),"MUNC"])
}
bruxelles[is.na(bruxelles[["MU_ID"]]),][["NAME_FRE"]]
bruxelles@data[(bruxelles@data[["NAME_FRE"]]=="REGIES"),"MU_ID"]="21013"
bruxelles@data[(bruxelles@data[["NAME_FRE"]]=="GIBET"),"MU_ID"]="21006"
bruxelles@data[(bruxelles@data[["NAME_FRE"]]=="LUXEMBOURG (PLACE DE)"),"MU_ID"]="21009"
bruxelles@data[(bruxelles@data[["NAME_FRE"]]=="CITE-JARDIN"),"MU_ID"]="21010"
bruxelles@data[,"SDDC"]=paste(bruxelles@data[,"MU_ID"], bruxelles@data[,"SDDC"], sep="")
bruxelles@data=merge(bruxelles@data,
 secteurs,
 by.x="SDDC", by.y="IdSec", all.x=T, sort=F)

Malheureusement, l'identifiant du secteur dans la data frame secteurs se présente sous le format "Code INS"+"Code secteur" alors que bruxelles@data contient l'information dans deux colonnes: "MU_ID" et "SDDC". Bien que la colonne "SDDC" contienne bien le code secteur, la colonne "MU_ID" est un code à 4 chiffres. Comme la vie est bien faite, il nous suffit de charger le fichier "UrbAdm_Mu.dbf" afin de récupérer le code INS. Avant cela, nous corrigeons quelques erreurs dans la data frame et fusionnons ensuite nos données. En d'autres termes, tout èa pour ajouter une colonne "densité" à bruxelles@data.

nArea=20
bruxelles@data[,"Area"]=cut(bruxelles@data[,"Dens"],
 quantile(bruxelles@data[!is.na(bruxelles@data[,"Dens"]),"Dens"], seq(0,1,1/nArea)),
 labels=F, include.lowest=T)
bruxelles@data[is.na(bruxelles@data[,"Area"]),"Area"]=nArea+1
area.pal=c(rev(as.character(heat.colors(nArea, alpha = .8))),
 "#B1B1B1CC")

A nouveau, nous définissons manuellement les catégories et leurs couleurs. Cette fois-ci, nous optons pour 20 catégories de densité contenant le même nombre de secteurs. La 21e catégorie reprend les secteurs pour lesquels nous n'avons pas de donnée. Nous utilisons heat.colors pour générer une palette de 20 couleurs (que nous inversons) et ajoutons une 21e pour le secteur sans donnée.

bruxelles@data=bruxelles@data[order(bruxelles@data[,"NbLigne"]),]
row.names(bruxelles@data)=bruxelles@data[,"NbLigne"]

Comme il semble qu'il soit important de maintenir l'ordre de bruxelles@data, nous nous en assurons.

xyplot(Y~X, ecoles,
 aspect="iso", par.settings=theme,
 panel = function(x, y, ...) {
  panel.polygonsplot(x=coordinates(bruxelles)[,1],
   y=coordinates(bruxelles)[,2],
   z=bruxelles@data[,"Area"],
   at=seq(min(bruxelles@data[,"Area"])-.5,
    max(bruxelles@data[,"Area"])+.5,
    by=1),
   subscripts=1:nrow(bruxelles@data),
   grid.polygons=bruxelles,
   sp.layout=NULL,
   col="#888888",
   region=TRUE,
   col.regions=area.pal,
   ...)
  panel.xyplot(x, y,
   pch=21,
   cex=ecoles[,"Point.cex"],
   col=as.character(ecoles[,"Point.col"]),
   fill=as.character(ecoles[,"Point.fill"]),
  ...)
 },
 scales=list(x=list(at=NULL), y=list(at=NULL)),
 xlab="", ylab="",
 xlim=c(141100, 158000), ylim=c(161400, 178300),
 legend=list(
  inside=list(fun=draw.key(key=list(title="Nombre\nd'éléves",
   cex=.6,
   points=list(pch=21,
    cex=c(200, 400, 600, 800)/(Ajust*min(ecoles[,"Taille"])),
    col=rgb(51,51,51, maxColorValue=255),
    fill=rgb(255,255,255, maxColorValue=255)),
   text=list(c("200", "400", "600", "800"),
    padding.text = 10,
    height = 10,
    cex=.8),
   padding.text=4
   )),
   x = .01, y = .3),
  inside=list(fun=draw.key(key=list(title="ISE\n(Ecoles)",
   cex=.6,
   points=list(pch=21,
    cex=2,
    col=as.character(class.pal[,"Point.col"]),
    fill=as.character(class.pal[,"Point.fill"])),
   text=list(class.lab, cex=.8),
   padding.text=3
   )),
   x = .01, y = .99),
  right=list(fun=draw.colorkey(key=list(col=area.pal[1:nArea],
   labels=as.character(c(0,
    do.call("c",by(bruxelles@data, bruxelles@data[,"Area"], simplify = F,
     function(x) round(max(x["Dens"]), 0)))[-(nArea+1)])),
   at=1:(nArea+1),
   width=1.5,
   tck=0,
   axis.line=list(lwd=1)
  )))
))

Pour mettre à jour le graphique, nous devons préciser trois choses: le vecteur reprenant les catégories ("z"), les valeurs qui indiquent un changement de classe ("at") et la palette à utiliser ("col.regions"). Nous ajoutons également une légende.

carte 7

Finitions et impressions

graph=xyplot(Y~X, ecoles,
 aspect="iso", par.settings=theme,
 panel = function(x, y, ...) {
  panel.polygonsplot(x=coordinates(bruxelles)[,1],
   y=coordinates(bruxelles)[,2],
   z=bruxelles@data[,"Area"],
   at=seq(min(bruxelles@data[,"Area"])-.5,
    max(bruxelles@data[,"Area"])+.5,
    by=1),
   subscripts=1:nrow(bruxelles@data),
   grid.polygons=bruxelles,
   sp.layout=NULL,
   col="#888888",
   region=TRUE,
   col.regions=area.pal,
   ...)
  panel.xyplot(x, y,
   pch=21,
   cex=ecoles[,"Point.cex"],
   col=as.character(ecoles[,"Point.col"]),
   fill=as.character(ecoles[,"Point.fill"]),
  ...)
  panel.segments(x0=156500, x1=157500, y0=161500, y1=161500,
   col=rgb(51,51,51, maxColorValue=255),
   lwd=2)
  panel.text(x=157000, y=161750, labels="1km", cex=.65)
  panel.text(x=157800, y=176500, labels="Densité (Secteurs)", cex=.9, srt=270)
  panel.text(x=141500, y=161800, cex=.6, pos = 4,
   labels="Fond de carte : Brussels UrbIS© - Distribution & Copyright CIRB")
  panel.text(x=141500, y=161600, cex=.6, pos = 4,
   labels="Données enseignement (2011-12): ETNIC - Communauté franèaise")
 },
 scales=list(x=list(at=NULL), y=list(at=NULL)),
 xlab="", ylab="",
 xlim=c(141100, 158000), ylim=c(161400, 178300),
 legend=list(
  inside=list(fun=draw.key(key=list(title="Nombre\nd'éléves",
   cex=.6,
   points=list(pch=21,
    cex=c(200, 400, 600, 800)/(Ajust*min(ecoles[,"Taille"])),
    col=rgb(51,51,51, maxColorValue=255),
    fill=rgb(255,255,255, maxColorValue=255)),
   text=list(c("200", "400", "600", "800"),
    padding.text = 10,
    height = 10,
    cex=.8),
   padding.text=4
   )),
   x = .01, y = .3),
  inside=list(fun=draw.key(key=list(title="ISE\n(Ecoles)",
   cex=.6,
   points=list(pch=21,
    cex=2,
    col=as.character(class.pal[,"Point.col"]),
    fill=as.character(class.pal[,"Point.fill"])),
   text=list(class.lab, cex=.8),
   padding.text=3
   )),
   x = .01, y = .99),
  right=list(fun=draw.colorkey(key=list(col=area.pal[1:nArea],
   labels=as.character(c(0,
    do.call("c",by(bruxelles@data, bruxelles@data[,"Area"], simplify = F,
     function(x) round(max(x["Dens"]), 0)))[-(nArea+1)])),
   at=1:(nArea+1),
   width=1.5,
   tck=0,
   axis.line=list(lwd=1)
  )))
))

Nous ajoutons quelques éléments pour finaliser notre graphique: un titre, des références et une échelle. Nous enregistrons le graphique dans l'objet graph.

svg(file.path(root.dir,"map.svg"), width = 8, height = 8)
print(graph)
dev.off()
pdf(file.path(root.dir,"map.pdf"), width = 8, height = 8, onefile=FALSE)
print(graph)
dev.off()
png(file.path(root.dir,"map.png"),
 width = 6000, height = 6000, res = 720,
 bg = "transparent")
print(graph)
dev.off()

Il est à présent possible d'imprimer ce graphique sous différents formats (PNG, PDF ou SVG).

carte 8