Réaliser un plan avec R et superposer des données

Utilisation de la librairie spatiale de R, sp.

La librairie sp permet de jongler avec des données spatiales au format rectangulaire (dataframe : juste des colonnes).

Elle présente des fonctions intéressantes que l'on peut relier à d'autres librairies (rgeos, cartography) pour réaliser tous les types de cartes.

1- Premier pas pour tracer une carte

1.1 - Charger des données et les apprêter

Charger des données qui caractérisent différents points : exemple de la Meuse

library("sp")

data("meuse.grid")

colnames(meuse.grid)

plot(meuse.grid$x,meuse.grid$y,pch=3)

Pour l'instant, on ne dispose que d'un ensemble de points qui couvre l'ensemble de la Meuse.

Il faudrait remplacer ces points par des pixels pour avoir un raster.

On notera aussi des coordonnées étranges qui justifient de changer de système géodésique si on désire des longitudes/latitudes classiques.

encadré - Le principe de la fonction gridded : remplacer les points par des pixels

# just 9 points on a grid:

x <- c(1,1,1,2,2,2,3,3,3)

y <- c(1,2,3,1,2,3,1,2,3)

xy <- cbind(x,y)

S <- SpatialPoints(xy)

class(S)

plot(S)

gridded(S) <- TRUE

gridded(S)

class(S)

summary(S)

plot(S)

# Conversion en objet spatial

coordinates(meuse.grid) <- ~x+y

is(meuse.grid) ; class(meuse.grid)

plot(meuse.grid)

# Création d'un pixel pour chaque point

# Faire un grilllage... gridded()

gridded(meuse.grid) <- TRUE

 plot(meuse.grid["dist"]) # Il faut indiquer une dimension qui servira à la coloration, ici le paramètre dist

1.2 - Réalisation des comptages/moyennes par secteur

Chargeons maintenant des données avec des points ponctuels : ici des analyses sur des prélèvement géolocalisés dans la Meuse :

data("meuse")

# meuse

# ce jeu de données donne des différents résultats d'analyse pour

# différents endroits de la Meuse

colnames(meuse)

meuse$cadmium

# La fonction coordinates va convertir la dataframe Meuse

# en "SpatialPointsDataFrame" avec regroupement des coordonnées x, y

# dans une colonne coordinates

coordinates(meuse) <- ~x+y

is(meuse)

Pour ces différents points, les présences de certains éléments (cadmium, zinc, etc) ont été analysées.

Segmenter une carte en régions en fonction de valeurs ou d'indice.

Ici, on pourrait découper la Meuse en 4 zones selon les valeurs du paramètres dist.

On va d'abord trier les valeurs de dist en 4 catégories.

# Utiliser la fonction cut pour classer les valeurs du paramètres dist dans les

# différents intervalles de 0 à 0,25 - 0,25 à 0,5 - 0,5 à 0,75... etc.

i = cut(meuse.grid$dist, c(0,.25,.5,.75,1), include.lowest = TRUE)

library("rgeos")

# Découper la Meuse en fonction des valeurs de i

ab = gUnaryUnion(as(meuse.grid, "SpatialPolygons"), i) # On peut essayer aussi meuse.grid$part.a à la place de i

plot(ab)

Combiner maintenant les résultats d'analyse des éléments avec la carte découpé que j'ai...

plot(ab)

plot(meuse,add=T,lwd=2,col="red")

Il serait bien d'avoir la composition moyenne de chacune de ces 4 régions. Problème : j'ai plusieurs mesure par aire. Il va falloir que je calcule donc automatiquement la moyenne pour chaque secteur. C'est le rôle de la fonction aggregate().

x = aggregate(meuse[c("zinc")], ab, mean)

x = aggregate(meuse[c("zinc", "copper","cadmium")], ab, mean)

spplot(x)

On indique à aggregate() 3 paramètres : la donnée à agréger, la segmentation correspondant et la méthode (ici, on fait la moyenne par secteur : mean).

On voit ici un petit problème... Automatiquement, la barre de coloration utilise la même échelle...

1.3 - Travailler la mise en forme

Un code qui récapitule ce qui a été fait ci-dessus en un simple copier-coller :

library("sp")

data("meuse.grid")

coordinates(meuse.grid) <- ~x+y

gridded(meuse.grid) <- TRUE

data("meuse")

coordinates(meuse) <- ~x+y

i = cut(meuse.grid$dist, c(0,.25,.5,.75,1), include.lowest = TRUE)

library("rgeos")

# Découper la Meuse en fonction des valeurs de i

ab = gUnaryUnion(as(meuse.grid, "SpatialPolygons"), i) # On peut essayer aussi meuse.grid$part.a à la place de i

x = aggregate(meuse[c("cadmium")], ab, mean)

spplot(x)

Définir des limites en x et en y

spplot(x, "cadmium",

  cuts=10,

  xlim=c(175000,185000),

  ylim=c(327000,335000))# définit le nombre d'intervalles de coloration

Changer la palette de couleurs et le découpage de la barre de légende

library(RColorBrewer) # pour brewer.pal( )

mypalette <- brewer.pal(8, 'Reds') # définit des gammes de couleurs

spplot(x, "cadmium",

  cuts=7,# définit le nombre d'intervalles de coloration

 col.regions = mypalette)

Compiler toutes les couches afin de faciliter l'affichage

scale = list("SpatialPolygonsRescale", layout.scale.bar(),

      offset = c(178600,332490), scale = 500, fill=c("transparent","black"))

text1 = list("sp.text", c(178600,332590), "0")

text2 = list("sp.text", c(179100,332590), "500 m")

arrow = list("SpatialPolygonsRescale", layout.north.arrow(),

     offset = c(178750,332000), scale = 400)

couches = list(scale,text1,text2,arrow)

mypalette <- brewer.pal(8, 'RdYlBu') # dé nit des gammes de couleurs

spplot(x, "cadmium", do.log=T,

 key.space=list(x=0.1,y=0.93,corner=c(0,1)),

 sp.layout=couches,

 main = "Zinc (top soil)",

 col.regions = mypalette)

spplot(x, "cadmium", #do.log=T,

 cut=7,

 sp.layout=couches,

 main = "Cadmium)",

 col.regions = mypalette)

2- Dessiner un plan manuellement avec sp

Il est possible de réaliser des plans spatiaux avec R project. Le pack "sp" (cf. librairie sp ou package sp) permet en effet d'ajouter des fonctions de gestion de l'espace à R.

1- Installer la librairie sp qui permet de gérer les données spatiales

install.packages("sp") ; library("sp")

2- Dessiner des polygones

x1 <- c(0, 72,72,0) ; y1 <- c(0,0,16.2,16.2) # Coordonnées x et y d'une salle carrée

c1 <- data.frame(x1, y1)

p1 <- Polygon(coords = c1, hole = F);

P1 <- Polygons(srl = list(p1), ID = "Aile de batiment") # Conversion en polygone

3- Combiner les polygones

Dans le cas où on aurait un deuxième polygone P2.

SP <- SpatialPolygons(Srl = list(P1,P2))

4. Tracer les polygones

plot(SP,col=temperatures,border = "white")

Pour aller plus loin avec sp.

Exemple de plan d'un bâtiment illustrant la température dans les différentes salles

Cet exemple permet de visualiser les températures moyennes relevées en hiver dans un bâtiment universitaire.

(Réalisé avec la participation de D. Paretour, L. Blois et R. Legleau)

Plan d'un bâtiment réalisé avec le logiciel R

install.packages("sp")

library("sp")

#############################################################################################

Compiler toutes les salles d'un bâtiment sous formes de polygones.

#############################################################################################

# Aile de Génie Biologique

x1 <- c(0, 72,72,0) ; y1 <- c(0,0,16.2,16.2) ; (c1 <- data.frame(x1, y1));p1 <- Polygon(coords = c1, hole = F);P1 <- Polygons(srl = list(p1), ID = "Aile de génie biologique")

# Couloir de génie biologique

x2 <- c(0, 72,72,0) ; y2 <- c(6.8,6.8,8.5,8.5) ; (c2 <- data.frame(x2, y2)) ;p2 <- Polygon(coords = c2, hole = F);P2 <- Polygons(srl = list(p2), ID = "Couloir de génie biologique")

# Toilettes 1

x3 <- c(0,1.5,1.5,0) ; y3 <- c(0,0,6.8,6.8) ; (c3 <- data.frame(x3, y3)) ;p3 <- Polygon(coords = c3, hole = F);P3 <- Polygons(srl = list(p3), ID = "Toilettes 1")

# cage escalier 1

x4 <- c(1.5,5,5,1.5) ; y4 <- c(0,0,6.8,6.8) ; (c4 <- data.frame(x4, y4));p4 <- Polygon(coords = c4, hole = F);P4 <- Polygons(srl = list(p4), ID = "cage escalier 1")

# prof 4

x5 <- c(5,8,8,5) ; y5 <- c(0,0,6.8,6.8) ; (c5 <- data.frame(x5, y5)) ;p5 <- Polygon(coords = c5, hole = F);P5 <- Polygons(srl = list(p5), ID = "prof 4")

# salle informatique 1

x6 <- c(8,11.3,11.3,8) ; y6 <- c(0,0,6.8,6.8) ; (c6 <- data.frame(x6, y6)) ;p6 <- Polygon(coords = c6, hole = F);P6 <- Polygons(srl = list(p6), ID = "salle informatique 1")

# salle informatique 2

x7 <- c(11.3,17.5,17.5,11.3) ; y7 <- c(0,0,6.8,6.8) ; (c7 <- data.frame(x7, y7));p7 <- Polygon(coords = c7, hole = F);P7 <- Polygons(srl = list(p7), ID = "salle informatique 2")

# couloir sortie 1

x8 <- c(17.5,20.7,20.7,17.5) ; y8 <- c(0,0,6.8,6.8) ; (c8 <- data.frame(x8, y8)) ;p8 <- Polygon(coords = c8, hole = F);P8 <- Polygons(srl = list(p8), ID = "couloir sortie 1")

# reunions professeurs

x9 <- c(20.7,24,24,20.7) ; y9 <- c(0,0,6.8,6.8) ; (c9 <- data.frame(x9, y9)) ;p9 <- Polygon(coords = c9, hole = F);P9 <- Polygons(srl = list(p9), ID = "reunions professeurs")

# biochimie

x10 <- c(24,33.5,33.5,24) ; y10 <- c(0,0,6.8,6.8) ; (c10 <- data.frame(x7, y7));p10 <- Polygon(coords = c10, hole = F);P10 <- Polygons(srl = list(p10), ID = "biochimie")

# reserves 1

x11 <- c(33.5,36.8,36.8,33.5) ; y11 <- c(0,0,6.8,6.8) ; (c11 <- data.frame(x11, y11)) ;p11 <- Polygon(coords = c11, hole = F);P11 <- Polygons(srl = list(p11), ID = "reserves 1")

# travaux 1

x12 <- c(36.8,43,43,36.8) ; y12 <- c(0,0,6.8,6.8) ; (c12 <- data.frame(x12, y12)) ;p12 <- Polygon(coords = c12, hole = F);P12 <- Polygons(srl = list(p12), ID = "travaux 1")

# microbiologie 1

x13 <- c(43,52.7,52.7,43) ; y13 <- c(0,0,6.8,6.8) ; (c13 <- data.frame(x13, y13));p13 <- Polygon(coords = c13, hole = F);P13 <- Polygons(srl = list(p13), ID = "microbiologie 1")

# laverie

x14 <- c(54.3,57.5,57.5,54.3) ; y14 <- c(0,0,2.3,2.3) ; (c14 <- data.frame(x14, y14)) ;p14 <- Polygon(coords = c14, hole = F);P14 <- Polygons(srl = list(p14), ID = "laverie")

# preparation microbio

x15 <- c(52.7,57.5,57.5,52.7) ; y15 <- c(0,0,4.5,4.5) ; (c15 <- data.frame(x15, y15)) ;p15 <- Polygon(coords = c15, hole = F);P15 <- Polygons(srl = list(p15), ID = "preparation microbio")

# vestiaire microbio

x16 <- c(52.7,57.5,57.5,52.7) ; y16 <- c(4.5,4.5,6.8,6.8) ; (c16 <- data.frame(x16, y16));p16 <- Polygon(coords = c16, hole = F);P16 <- Polygons(srl = list(p16), ID = "vestiaire microbio")

# microbiologie 2

x17 <- c(57.5,67,67,57.5) ; y17 <- c(0,0,6.8,6.8) ; (c17 <- data.frame(x17, y17)) ;p17 <- Polygon(coords = c17, hole = F);P17 <- Polygons(srl = list(p17), ID = "microbiologie 2")

# cage escalier 2

x18 <- c(67,70,70,67) ; y18 <- c(0,0,6.8,6.8) ; (c18 <- data.frame(x18, y18)) ;p18 <- Polygon(coords = c18, hole = F);P18 <- Polygons(srl = list(p18), ID = "cage escalier 2")

# toilettes 2

x19 <- c(70,72,72,70) ; y19 <- c(0,0,6.8,6.8) ; (c19 <- data.frame(x19, y19));p19 <- Polygon(coords = c19, hole = F);P19 <- Polygons(srl = list(p19), ID = "toilettes 2")

# prof 1

x20 <- c(0,3.2,3.2,0) ; y20 <- c(8.5,8.5,16.2,16.2) ; (c20 <- data.frame(x20, y20)) ;p20 <- Polygon(coords = c20, hole = F);P20 <- Polygons(srl = list(p20), ID = "prof 1")

# prof 2

x21 <- c(3.2,6.5,6.5,3.2) ; y21 <- c(8.5,8.5,16.2,16.2) ; (c21 <- data.frame(x21, y21)) ;p21 <- Polygon(coords = c21, hole = F);P21 <- Polygons(srl = list(p21), ID = "prof 2")

# prof 3

x22 <- c(6.5,9.8,9.8,6.5) ; y22 <- c(8.5,8.5,16.2,16.2) ; (c22 <- data.frame(x22, y22));p22 <- Polygon(coords = c22, hole = F);P22 <- Polygons(srl = list(p22), ID = "prof 3")

# secretariat

x23 <- c(9.8,14.5,14.5,9.8) ; y23 <- c(8.5,8.5,16.2,16.2) ; (c23 <- data.frame(x23, y23)) ;p23 <- Polygon(coords = c23, hole = F);P23 <- Polygons(srl = list(p23), ID = "secretariat")

# direction

x24 <- c(14.5,17.5,17.5,14.5) ; y24 <- c(8.5,8.5,16.2,16.2) ; (c24 <- data.frame(x24, y24)) ;p24 <- Polygon(coords = c24, hole = F);P24 <- Polygons(srl = list(p24), ID = "direction")

# couloir sortie 2

x25 <- c(17.5,22.5,22.5,17.5) ; y25 <- c(8.5,8.5,16.2,16.2) ; (c25 <- data.frame(x25, y25));p25 <- Polygon(coords = c25, hole = F);P25 <- Polygons(srl = list(p25), ID = "couloir sortie 2")

# chimie 1

x26 <- c(22.5,32,32,22.5) ; y26 <- c(8.5,8.5,16.2,16.2) ; (c26 <- data.frame(x26, y26)) ;p26 <- Polygon(coords = c26, hole = F);P26 <- Polygons(srl = list(p26), ID = "chimie 1")

# reserves 2

x27 <- c(32,35,35,32) ; y27 <- c(8.5,8.5,16.2,16.2) ; (c27 <- data.frame(x27, y27)) ;p27 <- Polygon(coords = c27, hole = F);P27 <- Polygons(srl = list(p27), ID = "reserves 2")

# chimie 2

x28 <- c(35,44.5,44.5,35) ; y28 <- c(8.5,8.5,16.2,16.2) ; (c28 <- data.frame(x28, y28));p28 <- Polygon(coords = c28, hole = F);P28 <- Polygons(srl = list(p28), ID = "chimie 2")

# travaux dirigés

x29 <- c(44.5,49.5,49.5,44.5) ; y29 <- c(8.5,8.5,16.2,16.2) ; (c29 <- data.frame(x29, y29)) ;p29 <- Polygon(coords = c29, hole = F);P29 <- Polygons(srl = list(p29), ID = "travaux dirigés")

# biologie animale

x30 <- c(49.5,59,59,49.5) ; y30 <- c(8.5,8.5,16.2,16.2) ; (c30 <- data.frame(x30, y30)) ;p30 <- Polygon(coords = c30, hole = F);P30 <- Polygons(srl = list(p30), ID = "biologie animale")

# reserves 3

x31 <- c(59,62,62,59) ; y31 <- c(8.5,8.5,16.2,16.2) ; (c31 <- data.frame(x31, y31));p31 <- Polygon(coords = c31, hole = F);P31 <- Polygons(srl = list(p31), ID = "reserves 3")

# biologie végétale

x32 <- c(62,72,72,62) ; y32 <- c(8.5,8.5,16.2,16.2) ; (c32 <- data.frame(x32, y32)) ;p32 <- Polygon(coords = c32, hole = F);P32 <- Polygons(srl = list(p32), ID = "biologie végétale")

# Conversion en polygone spatial

SP <- SpatialPolygons(Srl = list(P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13,P14,P15,P16,P17,P18,P19,P20,P21,P22,P23,P24,P25,P26,P27,P28,P29,P30,P31,P32))

Définir des températures pour les différentes salles P1 à P32.

temperatures = c(rep(20,20),21,26,25,24,rep(20,5),18,22,18)

temperatures[temperatures<19.5] <- "#0080FF"

temperatures[temperatures>=19.5&temperatures<=20.5] <- "#2EFE9A"

temperatures[temperatures>20.5&temperatures<=22] <- "#F5D0A9"

temperatures[temperatures>22] <- "#FE2E2E"

# A explorer rampe <- colorRampPalette(c("#0080FF", "#2EFE9A", "#F5D0A9","#FE2E2E"))

Tracer le plan

plot(SP,col=temperatures,border = "white")

Légender

text(mean(x1),mean(y1),"Couloir",adj=c(0.5,0.5),cex=0.5)

text(mean(x4),mean(y4),"Cage\nescalier",adj=c(0.5,0.5),cex=0.5)

text(mean(x5),mean(y5),"Bureau\nenseignant",adj=c(0.5,0.5),cex=0.5)

text(mean(x21),mean(y21),"Bureau\nenseignant",adj=c(0.5,0.5),cex=0.5)

text(mean(x22),mean(y22),"Bureau\nenseignant",adj=c(0.5,0.5),cex=0.5)

text(mean(x23),mean(y23),"Secrétariat",adj=c(0.5,0.5),cex=0.5)

text(mean(x6),mean(y6),"Bureau\nenseignant",adj=c(0.5,0.5),cex=0.5)

text(mean(x7),mean(y7),"Salle\ninformatique",adj=c(0.5,0.5),cex=0.5)

text(mean(x10),mean(y10),"Biochimie",adj=c(0.5,0.5),cex=0.5)

text(mean(x24),mean(y24),"Direction",adj=c(0.5,0.5),cex=0.5)

Mise à jour de la page septembre 2020