Séance 3 : Sous-ensembles et visualisation avec tidyverse
1 - Introduction
Au cours de cette séance, nous allons continuer à nous familiariser avec les packages du tidyverse et ggplot2 pour explorer des données spatiales à l’échelle d’un site. Nous utiliserons le même jeu de données jdd que dans les séances 1 et 2.
Lors de la dernière séance, nous avons utilisé la fonction plot() native (ou plutôt générique) de R.
Aujourd’hui nous allons utiliser la fonction ggplot(), issue de la bibliothèque du même nom
dédiée à la visualisation dans R mais aussi dans d’autres languages de programmation comme Python.
La fonction ggplot() est relativement intuitive, car elle fonctionne comme un système
de calques dans illustrator. Par exemple:
# utilisons le jeu de données d'exemple "iris"
head(iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosap <- ggplot(data=iris, # d'abord spécifier le jeu de données
aes(x=Petal.Width, # puis les variables x
y=Petal.Length, # et y
color=Species)) + # on peut faire varier les couleurs et
# les formes des points selon des variables
geom_point() + # geom_point pour un nuage de point,
# geom_bar pour des graph en barres,
# geom_boxplot pour boîtes à moustaches, etc.
theme_linedraw() # on peut aussi préciser un thème,
# et toutes sortes d'autres choses...
p # on a fait un objet graphique, on peut le projetter en tapant son nom
Note: en utilisant d’autres bibliothèques on peut ajouter des outils. Par
exemple la bibliothèque ggalt a une fonction geom_encircle() que l’on
peut appliquer à un sous-ensemble du jeu de données:
setosa <- iris %>% dplyr::filter(Species == "setosa")
virginica <- iris %>% dplyr::filter(Species == "virginica")
p +
geom_encircle(data = setosa, linetype = 1) +
geom_encircle(data = virginica, linetype = 1)
2 - Créer un plan et un carroyage
Avant toute chose, revenons à notre jeu de données de la dernière fois:
jdd <- as_tibble(read.csv("data/Data_exemple.csv", header = TRUE, sep = ","))
head(jdd)## # A tibble: 6 × 8
## numero type x_cm y_cm longueur_cm largeur_cm silex raccord
## <int> <chr> <int> <int> <dbl> <dbl> <chr> <int>
## 1 1 gravette 10 60 5 1.5 type A NA
## 2 2 gravette 30 50 5.9 1.8 type C NA
## 3 3 picardie 50 80 3.2 1.2 type B 1
## 4 4 percoir 260 100 4.7 2.1 type B NA
## 5 5 burin 240 120 6 2.4 type B 2
## 6 6 burin diedre 280 150 8 3.5 type A NATout a l’air en ordre ! Nous pouvons donc commencer :)
Dans un premier temps, nous allons essayer de reproduire un plan de répartition spatiale X, Y, déjà réalisé avec QGIS et Powerpoint. Vous pouvez consulter ce plan dans le dossier figures du projet, sous le nom de plan_exemple.pdf. L’objectif de cet exercice est de créer un code que vous pourrez ensuite adapter à d’autres données. Cela vous permettra:
1) de gagner un temps précieux, le code ne prenant que quelques secondes pour générer une figure ;
2) d’avoir une unité graphique dans votre document en utilisant toujours le même layout ;
3) d’organiser les différents éléments de vos figures toujours de la même façon.
Nous allons commencer par créer le carroyage du plan. D’après plan_exemple.pdf, la zone d’étude s’étend sur 16 carrés de 1 m de côté, organisés en 4 colonnes et 4 rangées. Vérifions tout d’abord les unités pour les coordonnées des points.
summary(jdd)## numero type x_cm y_cm longueur_cm
## Min. : 1 Length:33 Min. : 10.0 Min. : 50.0 Min. : 1.4
## 1st Qu.: 9 Class :character 1st Qu.: 50.0 1st Qu.: 90.0 1st Qu.: 4.5
## Median :17 Mode :character Median :210.0 Median :280.0 Median : 5.9
## Mean :17 Mean :174.8 Mean :217.6 Mean : 5.8
## 3rd Qu.:25 3rd Qu.:260.0 3rd Qu.:340.0 3rd Qu.: 7.4
## Max. :33 Max. :290.0 Max. :390.0 Max. :10.2
##
## largeur_cm silex raccord
## Min. :0.500 Length:33 Min. :1.0
## 1st Qu.:1.900 Class :character 1st Qu.:1.0
## Median :2.800 Mode :character Median :1.5
## Mean :2.976 Mean :1.5
## 3rd Qu.:3.700 3rd Qu.:2.0
## Max. :7.100 Max. :2.0
## NA's :29Il semble que les coordonnées soient en cm. En x, elles se situent entre 10 et 290 cm, et en y, entre 50 et 390 cm. Avec un carroyage basé sur des carrés de 1 m de côtés, les noeuds auront pour coordonnées (0;0), (0;100), (0;200), etc. Allons-y avec ggplot2 !
Tout d’abord, définissons les labels des axes :
carres.x <- c("A", "B", "C", "D")
carres.y <- c("1", "2", "3", "4")Créons à présent un objet ggplot avec des lignes tous les 100 cm en x et en y.
plan <- ggplot2::ggplot() +
# ces deux premières lignes définissent le thème graphique de la figure
theme_linedraw() +
theme(axis.ticks = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
panel.grid.major = element_line(),
panel.grid.minor = element_blank()) +
# ici, on définit les dimensions des axes : de 0 à 400 cm, avec la même échelle en x et en y (ratio = 1)
coord_fixed(ratio = 1, xlim = c(0, 400), ylim = c(0, 400)) +
# les deux lignes suivantes permettent de subdiviser les axes tous les 1 m
#scale_x_continuous(breaks = seq(0, 400, by = 100)) +
#scale_y_continuous(breaks = seq(0, 400, by = 100)) +
# enfin, ici on définit les labels des axes
ylab(element_blank()) +
xlab(element_blank()) +
annotate(geom = "text", x = c(seq(from = 50, to = 400, by = 100)), y = -10, label = carres.x) +
annotate(geom = "text", x = -10, y = c(seq(from = 50, to = 400, by = 100)), label = carres.y)
plan 
Nous pouvons à présent rajouter les points correspondant aux objets ainsi qu’un titre.
plan +
geom_point(data = jdd, mapping = aes(x = x_cm, y = y_cm)) +
labs(title = "Répartition spatiale des artéfacts - site de la Forêt")
La figure s’affiche dans l’onglet plot, nous pouvons donc suivre son évolution au fur et à mesure des modifications que l’on souhaite y apporter par la suite !
Ce plan reste toutefois peu informatif sur l’organisation des objets dans le site. Il est possible d’améliorer la figure en faisant varier la forme et la couleur des points en fonction des données qualitatives associées à chaque pièce.
plan +
geom_point(data = jdd, mapping = aes(x = x_cm, y = y_cm,
color = silex, shape = type,
size = 1)) + # on augmente un peu la taille des points pour plus de lisibilité !
labs(title = "Répartition spatiale des artéfacts - site de la Forêt")## Warning: The shape palette can deal with a maximum of 6 discrete values because
## more than 6 becomes difficult to discriminate; you have 10. Consider
## specifying shapes manually if you must have them.## Warning: Removed 12 rows containing missing values (`geom_point()`).
Malheureusement… la palette automatique ne contient pas assez de formes ! Nous allons donc devoir les définir manuellement dans les esthétiques de la figure. Pour cela nous allons attribuer à chaque valeur de type une forme.
Affichons pour commencer les levels de type (après l’avoir converti en factor).
jdd$type_fac <- factor(jdd$type,
# Généralement, les levels sont dans l'ordre alphabétique ; il est possible de changer leur ordre si vous le voulez !
levels = c("nucleus a eclats laminaires", "burin du Raysse",
"eclat", "lame", "chute de burin", "burin",
"burin diedre", "percoir", "picardie", "gravette"))
levels(jdd$type_fac)## [1] "nucleus a eclats laminaires" "burin du Raysse"
## [3] "eclat" "lame"
## [5] "chute de burin" "burin"
## [7] "burin diedre" "percoir"
## [9] "picardie" "gravette"On créé ensuite manuellement la palette de formes avec scale_shape_manual
plan +
geom_point(data = jdd, mapping = aes(x = x_cm, y = y_cm, color = silex,
fill = silex, shape = type_fac, size = 1)) +
scale_shape_manual(values = c("nucleus a eclats laminaires" = 7,
"burin du Raysse" = 12, "eclat" = 22, "lame" = 23,
"chute de burin" = 21, "burin" = 8,
"burin diedre" = 24, "percoir" = 10,
"picardie" = 14, "gravette" = 3)) +
labs(title = "Répartition spatiale des artéfacts - site de la Forêt")## Warning: Removed 1 rows containing missing values (`geom_point()`).
On peut aussi faire la même chose pour la palette de couleurs. Ici je vais utiliser les couleurs de la palette viridis.
# Ne pas oublier de convertir en factor !
jdd$silex_fac <- as.factor(jdd$silex)
levels(jdd$silex_fac)## [1] "type A" "type B" "type C"plan +
geom_point(data = jdd, mapping = aes(x = x_cm, y = y_cm, color = silex,
fill = silex, shape = type_fac, size = 1)) +
scale_shape_manual(values = c("nucleus a eclats laminaires" = 7,
"burin du Raysse" = 12, "eclat" = 22, "lame" = 23,
"chute de burin" = 21, "burin" = 8,
"burin diedre" = 24, "percoir" = 10,
"picardie" = 14, "gravette" = 3)) +
scale_color_manual(values = c("type A" = "#7AD151FF", "type B" = "#2A788EFF",
"type C" = "#440154FF")) +
scale_fill_manual(values = c("type A" = "#7AD151FF", "type B" = "#2A788EFF",
"type C" = "#440154FF")) +
labs(title = "Répartition spatiale des artéfacts - site de la Forêt")## Warning: Removed 1 rows containing missing values (`geom_point()`).
La figure est un peu plus facile à comprendre, mais le code pour la produire n’est pas très efficace… remarquez qu’il y a des répétitions (entre les l.178 et 180) et le fait de définir les palettes dans les fonctions complique la lecture.
Pour clarifier tout cela, nous pouvons créer des objets qui contiennent nos palettes. Nous n’aurons ensuite qu’à les appeler dans la fonction.
shape.pal <- c("nucleus a eclats laminaires" = 7, "burin du Raysse" = 12,
"eclat" = 22, "lame" = 23, "chute de burin" = 21, "burin" = 8,
"burin diedre" = 24, "percoir" = 10, "picardie" = 14,
"gravette" = 3)
color.pal <- c("type A" = "#7AD151FF", "type B" = "#2A788EFF", "type C" = "#440154FF")Pour alléger encore, on place les palettes et le titre dans l’objet plan de base.
plan.modif <- plan +
scale_shape_manual(values = shape.pal) +
scale_color_manual(values = color.pal) +
scale_fill_manual(values = color.pal) +
labs(title = "Répartition spatiale des artéfacts - site de la Forêt")
plan.modif +
geom_point(data = jdd, mapping = aes(x = x_cm, y = y_cm, color = silex,
fill = silex, shape = type_fac, size = 1)) # Même résultat !## Warning: Removed 1 rows containing missing values (`geom_point()`).
Modifions à présent la légende pour la rendre plus lisible
plan.modif +
geom_point(data = jdd, mapping = aes(x = x_cm, y = y_cm, color = silex,
fill = silex, shape = type_fac, size = 1)) +
# ici on va changer le titre des légendes
labs(color = "Matière première", shape = "Catégorie d'objet") +
# et là, nous pouvons supprimer les légendes que nous ne souhaitons pas voir apparaître dans la figure finale.
guides(fill = "none", size = "none")## Warning: Removed 1 rows containing missing values (`geom_point()`).
3 - Sélectionner des données à afficher
Dans cette nouvelle partie, nous allons apprendre comment sélectionner une partie des données, dans l’idée d’un requêtage permettant ensuite d’afficher une portion seulement des données initiales. Cette opération est très pratique par exemple pour réaliser des projections localisées sans avoir à modifier le tableau excel de base, ce qui est plus transparent et permet de limiter les erreurs possibles.
Découvrons comment mettre en oeuvre cette opération. Dans les quelques lignes suivantes, nous allons sélectionner les pièces qui ont fait l’objet d’un raccord ou d’un remontage afin de représenter la relation entre les pièces sur le plan.
Pour cela, nous allons créer un nouveau dataframe dans lequel nous allons “ranger” uniquement les lignes qui contiennent des données pour la variable raccord
# faisable avec la fonction is.na()
rac <- jdd[rowSums(is.na(jdd)) == 0,]
# aussi faisable avec la fonction complete.cases()
rac <- jdd[complete.cases(jdd$raccord), ]Les crochets permettent de dire à R que l’on va sélectionner des lignes (premier argument) ou des colonnes (deuxième argument) dans le data frame jdd. La fonction complete.cases() permet ensuite de repérer les cases qui ne sont pas vides. Cette ligne de code nous permet donc de sélectionner dans jdd les lignes pour lesquelles la valeur de raccord est renseignée.
Observons ce que cela donne :
head(jdd)## # A tibble: 6 × 10
## numero type x_cm y_cm longue…¹ large…² silex raccord type_…³ silex…⁴
## <int> <chr> <int> <int> <dbl> <dbl> <chr> <int> <fct> <fct>
## 1 1 gravette 10 60 5 1.5 type… NA gravet… type A
## 2 2 gravette 30 50 5.9 1.8 type… NA gravet… type C
## 3 3 picardie 50 80 3.2 1.2 type… 1 picard… type B
## 4 4 percoir 260 100 4.7 2.1 type… NA percoir type B
## 5 5 burin 240 120 6 2.4 type… 2 burin type B
## 6 6 burin diedre 280 150 8 3.5 type… NA burin … type A
## # … with abbreviated variable names ¹longueur_cm, ²largeur_cm, ³type_fac,
## # ⁴silex_fachead(rac)## # A tibble: 4 × 10
## numero type x_cm y_cm longu…¹ large…² silex raccord type_…³ silex…⁴
## <int> <chr> <int> <int> <dbl> <dbl> <chr> <int> <fct> <fct>
## 1 3 picardie 50 80 3.2 1.2 type… 1 picard… type B
## 2 5 burin 240 120 6 2.4 type… 2 burin type B
## 3 8 burin du Ray… 260 300 5.3 3.8 type… 1 burin … type B
## 4 9 chute de bur… 230 280 2 0.5 type… 2 chute … type A
## # … with abbreviated variable names ¹longueur_cm, ²largeur_cm, ³type_fac,
## # ⁴silex_facDans la data frame jdd, notez bien la présence de cases vides (NA) dans la colonne raccord. En revanche dans rac, seules les lignes pour lesquelles la valeur de raccord est renseignée (1 ou 2 en l’occurrence) sont présentes. Il n’y a donc que 4 pièces concernées.
Nous pouvons maintenant rajouter un lien entre les pièces en utilisant cette nouvelle data frame.
plan.modif +
geom_line(data = rac, mapping = aes(x = x_cm, y = y_cm, group = raccord),
linewidth = 0.5) + # l'argument 'group' permet de relier les pièces
# par groupe de mêmes valeurs
geom_point(data = jdd, mapping = aes(x = x_cm, y = y_cm, color = silex,
fill = silex, shape = type_fac, size = 1)) +
# légende
labs(color = "Matière première", shape = "Catégorie d'objet") +
guides(fill = "none", size = "none")## Warning: Removed 1 rows containing missing values (`geom_point()`).
Une autre façon de faire des sélections au sein d’une data frame est d’utiliser la fonction subset(). Cette fonction permet d’introduire des conditions plus complexes pour la sélection de lignes ou de colonnes, en utilisant des connecteurs logiques.
Par exemple, essayons de sélectionner les pièces situées entre 200 et 400 cm en x, et entre 300 et 500 cm en y.
jdd.zoom <- subset(jdd, x_cm >= 200 & x_cm <= 400
& y_cm >= 300 & y_cm <= 500)Pour créer jdd zoom, nous demandons R de renvoyer les lignes du data frame dont le x est supérieur ou égal (>=) à 200 ET (&) inférieur ou égal (<=) à 400 ET (&) dont le xy est supérieur ou égal (>=) à 300 ET (&) inférieur ou égal (<=) à 500.
Voyons le résultat sur le plan
plan.modif +
geom_point(data = jdd.zoom, mapping = aes(x = x_cm, y = y_cm, color = silex,
fill = silex, shape = type_fac, size = 1)) +
# titre et légende
labs(title = "Répartition spatiale des pièces de la zone nord",
color = "Matière première", shape = "Catégorie d'objet") +
guides(fill = "none", size = "none")
FELICITATIONS ! Vous avez fait vos premiers plots avec ggplot !!