ASPE - Traitements sur les IPR
Pascal Irz
Le 20 mai 2022
Objectif
Ce document donne un exemple de chaîne de traitement des données Aspe pour présenter des IPR sur une région et effectuer des vérifications au moyen de l’ancienne macro Excel de calcul de cet indice.
Il suppose que les opérations d’importation depuis le dump ont été effectuées et les dataframe sauvegardés au format .RData, comme détaillé dans le tuto dédié.
Les fonctions du package aspe dédiées à ces traitements sont identifiées par le préfixe ipr_.
Chargement des packages et des données
library(aspe)
library(tidyverse)
load(file = "raw_data/tables_sauf_mei_2021_10_21_11_44_01.RData")La dernière date de saisie ou de modification d’une opération dans la base utilisée dans le présent tutoriel est le 07 mars 2022.
Grandes lignes de l’analyse
Structure de la base
La base Aspe comprend plus d’une centaine de tables dont la plupart contiennent des référentiels associant des codes à des modalités. Par exemple la table ref_protocole associe le code pro_id à la modalité “Pêche partielle par points (grand milieu)”. Ces tables sont dites “périphériques”.
Chaque table comprend une clé primaire (suffixe "_id"), identifiant unique de chacune de ses lignes.
Les tables qui constituent la colonne vertébrale de la base sont au nombre de six :
- station
- point_prelevement
- operation
- prelevement_elementaire
- lot_poissons
- mesure_individuelle
Ces tables sont liées de manières hiérarchique. Ainsi, chaque mesure individuelle se rapporte à un lot, qui se rapporte à un prélèvement élémentaire, qui se rapporte à une opération … jusqu’à la station.
Pour savoir sur quelle station a été capturé ce goujon de 87mm, il faut donc remonter toute la chaîne.
Principe
Pour simplifier les différents traitements et éviter d’avoir à reproduire à chaque requête de sélection toute la chaîne de liaison des tables de la colonne vertébrale, on peut construire un tableau de correspondance des clés primaires. On n’incluera toutefois pas la table mesure_individuelle car elle comprend des millions de lignes, ce qui alourdirait considérablement le tableau. Le traitement des mesures individuelles sera abordé ultérieurement. Pour la suite, ce tableau de correspondance sera nommé “passerelle” en ce qu’il fait la jonction entre les différentes composantes de la base.
passerelle <- mef_creer_passerelle()
names(passerelle)## [1] "sta_id" "pop_id" "ope_id" "pre_id" "lop_id"Comme cette table est assez importante (5402743 lignes) il est utile de la restreindre au périmètre qui nous intéresse. Ici, nous allons effectuer une sélection sur le département (codé par les deux premiers chiffres du code Insee de la commune) au moyen de la fonction mef_ajouter_dept().
passerelle <- passerelle %>%
mef_ajouter_dept() %>%
filter(dept %in% c(22, 29, 35, 56))NB Le package {aspe} permet aussi effectuer une sélection géographique, par exemple sur la base des contours d’un bassin. Dans ce cas, on utilise un polygone pour faire la sélection (ne sont conservées que les observations à l’intérieur).
On peut aussi sélectionner uniquement les opérations réalisées dans le cadre des réseaux, ici RCS, RHP et RRP, au moyen de la fonction mef_ajouter_objectif() et des intitulés de la table ref_objectif.
| obj_id | obj_libelle |
|---|---|
| 1 | Sauvetage - Transfert |
| 2 | DCE – Référence |
| 3 | RCS – Réseau de Contrôle de Surveillance |
| 4 | RRP – Réseau de Référence Pérenne |
| 5 | RCO – Réseau Contrôle opérationnel |
| 6 | RHP – Réseau Hydrobiologique Piscicole |
| 7 | RNB – Réseau National de Bassin |
| 8 | RNSORMCE – Réseau National de Suivi des Opérations de Restauration hydroMorphologiques des Cours d’Eau |
| 9 | Suivi des populations d’anguilles |
| 10 | Suivi des populations de saumons |
| 11 | Suivi des populations de truites |
| 12 | Étude |
| 13 | Suivi des cours d’eau intermittents |
| 14 | Suivi de restauration |
| 15 | RCA - Réseau de contrôle additionnel |
| 16 | RCRLB - Réseau de Contrôle des ME en Respect des objectifs environnementaux du bassin LB |
passerelle <- passerelle %>%
mef_ajouter_objectif() %>%
filter(obj_libelle %in% c("RCS – Réseau de Contrôle de Surveillance",
"RRP – Réseau de Référence Pérenne",
"RHP – Réseau Hydrobiologique Piscicole"))Bilan des IPR par station
On va utiliser la passerelle pour assurer le lien entre des informations contenues dans différentes tables.
ipr <- passerelle %>%
mef_ajouter_ipr() %>%
mef_ajouter_ope_date() %>%
filter(ope_date > lubridate::dmy("01/01/2005")) %>%
mef_ajouter_libelle() %>%
droplevels()Les premières lignes de ce tableau sont les suivantes :
ipr %>% head() %>% DT::datatable()Si l’on veut les mêmes données mais avec une colonne par année :
ipr_1c_par_an <- ipr_pivoter_1colonne_par_an(ipr_df = ipr)Les noms des colonnes sont alors :
names(ipr_1c_par_an)## [1] "pop_id" "2005" "2006" "2007" "2008" "2009" "2010" "2011"
## [9] "2012" "2013" "2014" "2015" "2016" "2017" "2018" "2019"
## [17] "2020" "2021"Pour exporter le tableau en format utilisable simplement avec Excel :
write.csv2(ipr_1c_par_an, file = "processed_data/ipr_bzh_pdl_large.csv",
row.names = FALSE, na = "")Mise en forme pour la macro Excel de calcul de l’IPR
Les notes IPR présentes dans la table operation_ipr de la base ASPE sont calculées par le SEEE. Auparavant, ce calcul était effectué par une macro Excel qui prenait en entrée les caractéristiques de la station, la surface pêchée et les captures par espèce. En sortie elle fournissait les 7 métriques et l’indice agrégé.
NB Les noms des colonnes sont importés depuis le fichier contenant la macro Excel, nommé dans cet exemple
"MacroIPR_Sortie.xlsx"et contenu dans le sous-répertoireraw_data.
Première étape : structuration du tableau.
data <- passerelle %>%
ipr_formater_pour_macro(date_debut = '01/01/2020')
names(data)## [1] "ope_id" "coursdo" "sta_libelle_sandre"
## [4] "ope_date" "opi_param_surf" "opi_param_sbv"
## [7] "opi_param_ds" "opi_param_lar" "opi_param_pent"
## [10] "opi_param_prof" "opi_param_alt" "opi_param_tjuillet"
## [13] "opi_param_tjanvier" "opi_param_bassin" "ANG"
## [16] "CHA" "CHE" "EPT"
## [19] "GOU" "LOF" "LPP"
## [22] "ROT" "VAI" "BRO"
## [25] "OCL" "PER" "PES"
## [28] "TAN" "ABL" "BBG"
## [31] "BRB" "BRE" "CCO"
## [34] "GAR" "GRE" "PCC"
## [37] "PSR" "SAN" "VAR"
## [40] "TRF" "EPI" "SAT"
## [43] "LPM" "TRM" "HBG"
## [46] "PCH" "ABH" "VAN"
## [49] "SPI" "FLE" "BOU"
## [52] "CAG" "BRX" "MUP"
## [55] "SIL" "IDE" "CAX"
## [58] "CMI"Seconde étape : Sélection et renommage des variables.
data <- data %>%
ipr_renommer_pour_macro(fichier_macro = "../raw_data/MacroIPR_Sortie.xlsx")
names(data)## [1] "N° de code ou de référence"
## [2] "Nom du cours d'eau"
## [3] "Nom de la station"
## [4] "Date de l'opération"
## [5] "v1"
## [6] "Surface échantillonnée (SURF)"
## [7] "Surface du bassin versant drainé (SBV)"
## [8] "Distance à la source (DS)"
## [9] "Largeur moyenne en eau (LAR)"
## [10] "Pente du cours d'eau (PEN)"
## [11] "Profondeur moyenne (PROF)"
## [12] "Altitude (ALT)"
## [13] "Température moyenne de juillet (TJUILLET)"
## [14] "Température moyenne de janvier (TJANVIER)"
## [15] "Unité hydrologique (HU)"
## [16] "v2"
## [17] "ABLab"
## [18] "ANGab"
## [19] "BAFab"
## [20] "BAMab"
## [21] "BLNab"
## [22] "BOUab"
## [23] "BBBab"
## [24] "BROab"
## [25] "CASab"
## [26] "CCOab"
## [27] "CHAab"
## [28] "CHEab"
## [29] "EPIab"
## [30] "EPTab"
## [31] "GARab"
## [32] "GOUab"
## [33] "GREab"
## [34] "HOTab"
## [35] "LOFab"
## [36] "LOTab"
## [37] "LPPab"
## [38] "OBRab"
## [39] "PCHab"
## [40] "PERab"
## [41] "PESab"
## [42] "ROTab"
## [43] "SANab"
## [44] "SATab"
## [45] "SPIab"
## [46] "TANab"
## [47] "TOXab"
## [48] "TRFab"
## [49] "VAIab"
## [50] "VANab"Exportation du tableau mis en forme.
write.csv2(data, "processed_data/aspe_format_macro.csv",
row.names = FALSE,
na = "")Quelques mises en forme
Evolution par station dans un département
Définition du département.
mon_dept <- '22'Création du sous-jeu de données filtré pour le graphique.
ipr_dept <- ipr %>%
filter(dept == mon_dept) %>% # filtrage sur le département
mef_filtrer_nb_mini_annees(nb_mini_annees = 7, # suppression des points avec moins de 7 années
var_id = "pop_id") %>%
select(pop_libelle,
pop_id,
dept,
annee,
ope_id,
ipr,
cli_libelle) %>%
distinct()
ipr_dept %>%
head() %>%
DT::datatable()Pour les graphiques IPR il est utile d’affecter à chaque classe de qualité un code couleur. On peut simplement ajouter une colonne au référentiel des classes avec la fonction ipr_completer_classes_couleur().
classe_ipr <- classe_ipr %>%
ip_completer_classes_couleur()Graphique pour une station :
gg_temp_ipr(df_ipr = ipr_dept,
var_id_sta = pop_libelle,
station_sel = c("TRIEUX à PLESIDY"),
var_ipr = ipr,
sup_500m = FALSE)
Par défaut l’axe des ordonnées est orienté vers le bas (il est inversé pr rapport à un graphique “classique”) donc plus les points sont en haut du graphique plus la qualité du milieu est élevée. On peut choisir de changer l’orientation de l’axe des ordonnées (IPR) avec l’argument inv_y = FALSE. Si le document de sortie est au format HTML, on peut aussi rendre le graphique interactif au moyen de l’argument interactif = TRUE.
gg_temp_ipr(df_ipr = ipr_dept,
var_id_sta = pop_libelle,
station_sel = c("TRIEUX à PLESIDY"),
var_ipr = ipr,
sup_500m = FALSE,
inv_y = FALSE,
interactif = TRUE)Si l’argument station_sel n’est pas renseigné, toutes les stations du dataframe sont toutes conservées.
mon_graphique <-
gg_temp_ipr(
df_ipr = ipr_dept,
var_id_sta = pop_libelle,
var_ipr = ipr,
sup_500m = FALSE
)
mon_graphique
Si l’on veut sauvegarder le graphique (ici dans le sous-répertoire "processed_data"), il faut utiliser la fonction ggsave() Le choix de l’extension choisie (“png”, “jpg”, “.bmp”) détermine le format de sortie.
Le graphique de sortie étant un un objet de classe ggplot, il peut être modifié par certaines des fonctions de ce package. Par exemple, les différences mineures entre l’affichage à l’écran et la sauvegarde font que la taille de police n’est pas toujours bien adaptée au format d’export. Dans ce cas, on peut modifier la taille de certains éléments texte comme ici le nom de la station :
mon_graphique <- mon_graphique +
theme(strip.text.x = element_text(size = 7))NB la synthaxe du package
ggplot2utilise pour enchaîner les manipulations l’opérateur+au lieu du%>%ailleurs …
Chemin de sauvegarde et nommage du fichier (concaténation avec la fonction paste0()) :
mon_chemin <- paste0("processed_data/ipr_", mon_dept, ".png")
mon_chemin## [1] "processed_data/ipr_22.png"Sauvegarde :
ggsave(filename = mon_chemin,
plot = mon_graphique,
width = 15, # largeur
height = 25, # hauteur
units = 'cm') # unité.Pour automatiser la production des graphiques sur un ensemble de départements, on peut utiliser les fonctions de la famille map du package purrr afin d’appliquer les fonctions non plus à chaque fois sur un objet (ex : le nom du département ou le graphique) mais sur une liste contenant plusieurs objets (ex : les noms des départements ou les graphiques).
mes_depts <- c('35', '22', '56', '29')
depts_libelles <- departement %>%
filter(dep_code_insee %in% mes_depts) %>%
pull(dep_libelle) %>%
as.character()
depts_libelles## [1] "Côtes-d'Armor" "Finistère" "Ille-et-Vilaine" "Morbihan"Filtrage du jeu de données aux départements concernés, en excluant les stations avec moins de 7 années de données.
data <- ipr %>%
filter(dept %in% mes_depts) %>%
mef_filtrer_nb_mini_annees(nb_mini_annees = 7,
var_id = "pop_id")On scinde le tableau pour en avoir un par département. L’objet data_par_dept est une liste contenant un dataframe par département.
data_par_dept <- split(x = data,
f = data$dept)Production des graphiques. La fonction map() prend en arguments la fonction à appliquer ainsi que deux listes contenant les arguments de la fonction :
- Fonction à appliquer :
.f = gg_temp_ipr - La liste des
dataframe:.x = data_par_dept
mes_graphiques <- map(.x = data_par_dept,
.f = gg_temp_ipr,
var_id_sta = pop_libelle,
var_ipr = ipr,
sup_500m = FALSE)L’objet mes_graphiques est une liste contenant autant d’éléments qu’il y a de graphiques, soit un par département. Pour accéder au xième de ces graphiques on le sélectionne dans la liste au moyen de [[x]]
mes_graphiques[2:4]## $`29`
##
## $`35`
##
## $`56`
Pour sauvegarder ces graphiques, on poursuit avec les fonctions map() quand on a une seule liste d’arguments à appliquer à la fonction et map2() quand il y en a deux.
# mise des polices à l'échelle
mes_graphiques <- map(.x = mes_graphiques,
.f = function(x) {x + theme(strip.text.x = element_text(size = 7))})
# construction des chemins et noms de fichiers
mes_chemins <- paste0("processed_data/ipr_", mes_depts, ".png")
mes_chemins
# sauvegarde
map2(.x = mes_chemins,
.y = mes_graphiques,
.f = ggsave,
width = 16, # largeur
height = 25, # hauteur
units = 'cm')Comme tout se déroule à merveille, vous pouvez vérifier la présence des fichiers image sauvegardés dans le sous-répertoire de sortie.
Pourcentage de stations en bon état
ipr_grapher_pc_bon(ipr_dept, titre = "22") +
labs(x = "")
Analyse par métrique
Principe
Il faut dans un premier temps collecter les données non seulement de l’indice agrégé mais de chacune de ses sept métriques. On peut s’intéresser aux valeurs des métriques elles-mêmes, mais aussi aux valeurs prédites par les modèles sous-jacents.
Ci-dessus, un dataframe passerelle a déjà été constitué pour effectuer les liens entre tables. Il est nécessaire de le filtrer pour ne conserver que les stations qui nous intéressent.
Assemblage du dataframe
Création d’un vecteur contenant les identifiants des opérations.
mes_operations <- ipr_dept %>%
pull(ope_id)Préparation du dataframe en repartant de la passerelle.
ipr_metriques <- passerelle %>%
select(-pre_id, -lop_id) %>% # variables inutiles
distinct() %>% # suppression des lignes dupliquées
mef_ajouter_ope_date() %>% # ajout date et année
filter(ope_id %in% mes_operations, annee > 2009) %>% # filtrage sur les identifiants d'opérations ope_id et année
mef_ajouter_metriques() # ajout des métriquesLes noms des colonnes du tableau sont les suivants :
names(ipr_metriques)## [1] "sta_id" "pop_id" "ope_id" "dept"
## [5] "obj_id" "obj_libelle" "ope_date" "annee"
## [9] "ner_theorique" "nel_theorique" "nte_theorique" "dit_theorique"
## [13] "dio_theorique" "dii_theorique" "dti_theorique" "ner_observe"
## [17] "nel_observe" "nte_observe" "dit_observe" "dio_observe"
## [21] "dii_observe" "dti_observe" "ner" "nel"
## [25] "nte" "dit" "dio" "dii"
## [29] "dti". . . .
Représentations graphiques sur un point
Préparation des éléments nécessaires
Mise en forme du tableau
La représentation graphique en radar va faire appel au package {fmsb} qui nécessite une mise en forme assez spécifique des données.
library(fmsb)Choix du point.
mon_pop <- 41964pop_data <- ipr_metriques %>%
mef_ajouter_libelle() %>%
select(sta_id,
pop_id,
pop_libelle,
annee,
ner:dti) %>%
distinct() %>% # pour éviter les doublons sur l'année qui sert de nom de lignes donc ne pê dupliquée
mef_ipr_radar(pop_id = mon_pop)Le tableau mis en forme ne contient comme colonnes que les métriques. Les lignes ont été nommées d’après les années. Deux lignes ont été insérées en haut du tableau, qui servent à préciser les bornes (mini et maxi) pour chacune des métriques. Avec la fonction ipr_mef_radar() le mini est fixé à zéro et le maxi est la valeur la plus élevée de l’ensemble des métriques pour l’ensemble des années.
knitr::kable(pop_data, row.names = TRUE)| ner | nel | nte | dit | dio | dii | dti | |
|---|---|---|---|---|---|---|---|
| 1 | 4.0000000 | 4.0000000 | 4.0000000 | 4.000000 | 4.000000 | 4.0000000 | 4.000000 |
| 2 | 0.0000000 | 0.0000000 | 0.0000000 | 0.000000 | 0.000000 | 0.0000000 | 0.000000 |
| 2010 | 1.4397068 | 0.4236727 | 0.6061896 | 3.498915 | 1.358605 | 0.4414578 | 1.741727 |
| 2012 | 1.4529206 | 0.4089054 | 0.3019119 | 3.864439 | 0.552212 | 0.5762340 | 1.277662 |
| 2016 | 0.2101564 | 0.4146136 | 0.5882401 | 2.831486 | 2.405983 | 0.2850565 | 1.881109 |
| 2014 | 1.4307231 | 0.4287589 | 0.6254237 | 2.562405 | 2.092888 | 0.1510891 | 2.636702 |
| 2018 | 1.4410185 | 0.4227217 | 1.8198722 | 3.124998 | 3.343828 | 0.4595565 | 2.723360 |
| 2020 | 1.4504344 | 0.4130288 | 1.7845554 | 2.978081 | 3.489961 | 0.6711200 | 1.052445 |
Titre du graphique
mon_titre <- ipr_metriques %>%
mef_ajouter_libelle() %>%
filter(pop_id == mon_pop) %>%
slice(1) %>%
pull(pop_libelle)Palette de couleur
Le package {RColorBrewer} propose des palettes de couleurs “standard” contrastées ou dégradées. Chaque couleur est identifiée par un code hexadécimal. Par exemple en utilisant un générateur en ligne on obtient la correspondance entre la vouleur et son code.
Ici nous utilisons la palette dégradée allant du rouge au bleu via le jaune "RdYlBu". Il faut une couleur par contour (= par année représentée) donc autant de couleurs que d’années.
nb_annees <- nrow(pop_data)
# Vecteurs des couleurs de ligne
couleurs_traits <- RColorBrewer::brewer.pal(n = nb_annees, name = "RdYlBu")
couleurs_traits## [1] "#D73027" "#F46D43" "#FDAE61" "#FEE090" "#E0F3F8" "#ABD9E9" "#74ADD1"
## [8] "#4575B4"On peut ensuite définir utiliser le package {GISTools} pour définir la palette pour les remplissages qui est simplement la même que celle des traits, mais avec de la transparence.
couleurs_remplissage <- couleurs_traits %>%
GISTools::add.alpha(alpha = 0.1)Graphique radar
Le graphique est produit par la fonction radarchart du package fmsb qui admet de très nombreux arguments permettant de personnaliser l’aspect des axes et de leurs graduations, des étiquettes, etc. Pour avoir accéder à l’aide de la fonction :
?radarchartMalheureusement cette aide est très sommaire et manque d’exemples.
radarchart(pop_data,
axistype = 2,
#custom polygon
pcol = couleurs_traits, pfcol = couleurs_remplissage, plwd = 4, plty = 1,
#custom the grid
cglcol = "grey", cglty = 1, axislabcol = "grey20",
#custom labels
vlcex = 0.8,
title = mon_titre)
legend(x = 0.7, y = 1.55,
legend = rownames(pop_data[-c(1,2),]),
bty = "n", pch = 20 ,
col = couleurs_traits,
text.col = "grey80", cex = 1.2, pt.cex = 3)
Diagramme en bâtons
pop_data <- pop_data %>%
slice(-(1:2)) %>%
rownames_to_column(var = "annee") %>%
mutate(annee = as.integer(annee)) %>%
pivot_longer(cols = -annee,
names_to = "metrique",
values_to = "valeur")
ggplot(data = pop_data,
aes(x = annee,
y = valeur)) +
geom_bar(stat = "identity",
fill = couleurs_traits[1]) +
labs(x = "",
y = "Valeur de la métrique",
title = mon_titre) +
facet_wrap(~metrique)