#************************************************************************************************ #### Script du code pour le mémoire ------------------------------------------------------------- #************************************************************************************************ #************************************************************************************************ #### Ouverture des données et des packages + mise au propre ------------------------------------- #************************************************************************************************ library(readxl) library(tidyverse) library(dplyr) library(ggplot2) library(httr) library(readr) library(jsonlite) library(maps) library(plotly) library(lubridate) library(tidyr) library(stringr) library(RColorBrewer) library(purrr) library(viridis) data <- read_excel(path="Documents/donnees/Agences_2023_2024 copie.xlsx", col_names=T) head(data) # renommer les noms des colonnes pour enlever les espaces et que les noms soient + pratiques data <- data %>% rename(code_BIOMAE = `Code BIOMAE`, date = `Date J0`, campagne = `Campagne`, nom_stat = `Libellé station`, code_stat = `Code station`, plomb = `Plomb (Pb)`, nickel = `Nickel (Ni)`, cadmium = `Cadmium (Cd)`) # mettre les métaux en numérique et remplacer les , par des . au cas où data <- data %>% mutate(plomb = as.numeric(gsub(",", ".", plomb)), na.rm=T) data <- data %>% mutate(nickel = as.numeric(gsub(",", ".", nickel)), na.rm=T) data <- data %>% mutate(cadmium = as.numeric(gsub(",", ".", cadmium)), na.rm=T) #************************************************************************************************ #### Chercher les données manquantes sur Naïades: Calcium, Sodium, Magnésium, Longitude, Latitude ---- #************************************************************************************************ station=unique(data$code_stat) # Codes SANDRE: Ca, Na, Mg codes_ions <- c(1374, 1375, 1372) # Récupération par groupes de 10 stations (sinon l'URL est trop grande) groupes <- split(station, ceiling(seq_along(station) / 10)) naiades_list <- list() # liste vide # Boucle for pour aller chercher par groupe de stations dans l'URL for (i in seq_along(groupes)) { groupe <- groupes[[i]] stations_liste <- paste(groupe, collapse = ",") url <- paste0( "https://hubeau.eaufrance.fr/api/v2/qualite_rivieres/analyse_pc?", "code_station=", stations_liste, "&code_parametre=1374,1375,1372,1302", "&date_debut_prelevement=2023-01-01", "&date_fin_prelevement=2024-12-31", "&size=20000", "&format=json", # ne fonctionne pas avec le format .csv "&latitude", "&longitude") reponse <- GET(url) if (status_code(reponse) %in% c(200, 206)) { # on veut pas un code 400, on veut 206: contenu partiel ou 200: succès contenu <- fromJSON(content(reponse, as = "text", encoding = "UTF-8")) # on lit le fichier JSON naiades_list[[i]] <- contenu$data } } # Combiner tous les groupes naiades_raw <- bind_rows(naiades_list) # Garder uniquement les colonnes utiles naiades <- naiades_raw %>% select(code_station, date_prelevement, libelle_parametre, resultat, latitude, longitude) %>% group_by(code_station, libelle_parametre, latitude, longitude) %>% summarise(resultat = median(resultat, na.rm=T)) %>% pivot_wider(names_from = libelle_parametre, values_from = resultat, values_fn = median) %>% rename( code_stat = code_station, calcium = Calcium, sodium = Sodium, magnesium = `Magnésium`, pH = `Potentiel en Hydrogène (pH)`) # Jointure sur Code Station data_complet <- data %>% left_join(naiades, by = "code_stat") View(data_complet) #************************************************************************************************ #### Stats descriptives et graphiques de visualisation des données ------------------------------ #************************************************************************************************ # faire mon summary pour avoir les stats descriptives par station et par métaux summary_metaux <- data %>% pivot_longer( cols = c(plomb, nickel, cadmium), names_to = "metal", values_to = "concentration") %>% group_by(metal) %>% summarise(n = sum(!is.na(concentration)), min = min(concentration, na.rm = TRUE), q025 = quantile(concentration, 0.025, na.rm = TRUE), q25 = quantile(concentration, 0.25, na.rm = TRUE), moyenne = mean(concentration, na.rm = TRUE), mediane = median(concentration, na.rm = TRUE), ecart_type = sd(concentration, na.rm = TRUE), q75 = quantile(concentration, 0.75, na.rm = TRUE), q975 = quantile(concentration, 0.975, na.rm = TRUE), max = max(concentration, na.rm = TRUE), IC95_bas = t.test(concentration[!is.na(concentration)], conf.level = 0.95)$conf.int[1], IC95_haut = t.test(concentration[!is.na(concentration)], conf.level = 0.95)$conf.int[2]) print(summary_metaux) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Distribution groupée 3 en 1 en boxplot data_long <- data_complet %>% pivot_longer( cols = c(plomb, nickel, cadmium), names_to = "metal", values_to = "concentration" ) %>% filter(!is.na(concentration)) calculer_labels <- function(y) { stats <- quantile(y, probs = 0.5, na.rm = TRUE) return(data.frame( y = stats, label = round(stats, 3))) } BBAC <- data.frame(metal = c("nickel", "plomb", "cadmium"), bbac = c(0.4, 0.31, 0.025)) # Graphique ggplot(data_long, aes(x = metal, y = concentration, fill = metal)) + geom_boxplot(alpha = 0.7, outlier.colour = "grey1") + stat_summary(fun.data = calculer_labels, geom = "text", aes(label = after_stat(label)), hjust = -3, size = 5, color = "darkred", fontface = "bold") + geom_point(data = BBAC, aes(x = metal, y = bbac), shape = 4, color = "red", size = 6, stroke = 2, inherit.aes = FALSE) + labs(x = "Metal", y = "Concentration (mg/kg PF)") + scale_y_log10() + ylim(c(0, 1.5)) + theme_minimal() + scale_fill_manual(values = c("plomb" = "#21908CFF", "nickel" = "#B63679FF", "cadmium" = "#F98C0AFF")) + theme(axis.title = element_text(size = 16), axis.text = element_text(size = 16), legend.title = element_text(size = 16, face = "bold"), legend.text = element_text(size = 16)) #************************************************************************************************ #### Graphiques de la France avec les concentrations -------------------------------------------- #************************************************************************************************ france <- map_data("france") # Groupes de métaux: 0 -- BBAC -- seuil 25% -- seuil 50% -- seuil 75% -- au-delà seuil_plomb <- c(0, 0.31, 0.36525, 0.436, 0.6625, Inf) seuil_nickel <- c(0, 0.4, 0.456, 0.526, 0.64, Inf) seuil_cadmium <- c(0, 0.025, 0.032, 0.044, 0.076, Inf) # pour faire les légendes: on arrondit les nombres pour alléger label_plomb <- c("< 0.31", "0.31 - 0.37", "0.37 - 0.44", "0.44 - 0.66", "> 0.66") label_nickel <- c("< 0.40", "0.40 - 0.46", "0.46 - 0.53", "0.53 - 0.64", "> 0.64") label_cadmium <- c("< 0.025", "0.025 - 0.032", "0.032 - 0.044", "0.044 - 0.076", "> 0.076") # Maps des ions # Calcium ggplot() + geom_polygon(data = france, aes(x = long, y = lat, group = group), fill = "grey93", color = "white") + geom_point(data = naiades %>% filter(!is.na(calcium)) %>% arrange(calcium), aes(x = longitude, y = latitude, color = calcium), size = 3) + scale_color_gradient(low = "moccasin", high = "lightsalmon4") + labs(color = "Calcium (mg/L)") + theme_void() + theme(legend.title = element_text(size = 18), legend.text = element_text(size = 16)) # Magnesium ggplot() + geom_polygon(data = france, aes(x = long, y = lat, group = group), fill = "grey93", color = "white") + geom_point(data = naiades %>% filter(!is.na(magnesium)) %>% arrange(magnesium), aes(x = longitude, y = latitude, color = magnesium), size = 3) + scale_color_gradient(low = "lightsteelblue", high = "#39558CFF") + labs(color = "Magnésium (mg/L)") + theme_void() + theme(legend.title = element_text(size = 18), legend.text = element_text(size = 16)) # Sodium ggplot() + geom_polygon(data = france, aes(x = long, y = lat, group = group), fill = "grey93", color = "white") + geom_point(data = naiades %>% filter(!is.na(sodium)) %>% arrange(sodium), aes(x = longitude, y = latitude, color = sodium), size = 3) + scale_color_gradient(low = "lightpink", high = "pink4") + labs(color = "Sodium (mg/L)") + theme_void() + theme(legend.title = element_text(size = 18), legend.text = element_text(size = 16)) # Maps des métaux # Plomb data_plomb <- data_complet %>% filter(!is.na(plomb)) %>% mutate(classe = cut(plomb, breaks = seuil_plomb, labels = label_plomb, include.lowest = TRUE, right = TRUE)) ggplot() + geom_polygon(data = france, aes(x = long, y = lat, group = group), fill = "grey91", color = "white") + geom_point(data = data_plomb %>% arrange(classe), aes(x = longitude, y = latitude, color = classe), size = 3) + scale_color_manual(values = c("< 0.31" = "#AADC32FF", "0.31 - 0.37" = "#5DC863FF", "0.37 - 0.44" = "#21908CFF", "0.44 - 0.66" = "#3B528BFF", "> 0.66" = "#440154FF"), name = "Plomb (mg/kg PF)") + theme_void() + theme(legend.position = "right", legend.title = element_text(size = 18), legend.text = element_text(size = 16)) # Nickel data_nickel <- data_complet %>% filter(!is.na(nickel)) %>% mutate(classe = cut(nickel, breaks = seuil_nickel, labels = label_nickel, include.lowest = TRUE, right = TRUE)) ggplot() + geom_polygon(data = france, aes(x = long, y = lat, group = group), fill = "grey91", color = "white") + geom_point(data = data_nickel %>% arrange(classe), aes(x = longitude, y = latitude, color = classe), size = 3) + scale_color_manual(values = c("< 0.40" = "#FEC287FF", "0.40 - 0.46" = "#FB8861FF", "0.46 - 0.53" = "#B63679FF", "0.53 - 0.64" = "#51127CFF", "> 0.64" = "#000004FF"), name = "Nickel (mg/kg PF)") + theme_void() + theme(legend.position = "right", legend.title = element_text(size = 18), legend.text = element_text(size = 16)) # Cadmium data_cadmium <- data_complet %>% filter(!is.na(cadmium)) %>% mutate(classe = cut(cadmium, breaks = seuil_cadmium, labels = label_cadmium, include.lowest = TRUE, right = TRUE)) ggplot() + geom_polygon(data = france, aes(x = long, y = lat, group = group), fill = "grey91", color = "white") + geom_point(data = data_cadmium %>% arrange(classe), aes(x = longitude, y = latitude, color = classe), size = 3) + scale_color_manual(values = c("< 0.025" = "#F9C932FF", "0.025 - 0.032" = "#F98C0AFF", "0.032 - 0.044" = "#BB3754FF", "0.044 - 0.076" = "#56106EFF", "> 0.076" = "#000004FF"), name = "Cadmium (mg/kg PF)") + theme_void() + theme(legend.position = "right", legend.title = element_text(size = 18), legend.text = element_text(size = 16)) #************************************************************************************************ #### Boxplot: Calculs et graphes de distribution en fonction de la concentration en ions ---- #************************************************************************************************ # Calcium data_complet <- data_complet %>% mutate(classe_calcium = cut(calcium, breaks = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 130, 160, 190, 220), labels = c("0-10", "10-20", "20-30", "30-40", "40-50", "50-60", "60-70", "70-80", "80-90", "90-100", "100-130", "130-160", "160-190", "190-220"), include.lowest = TRUE)) # Magnésium data_complet <- data_complet %>% mutate(classe_magnesium = cut(magnesium, breaks = c(0, 3, 6, 10, 15, 20, Inf), labels = c("0-3", "3-6", "6-10", "10-15", "15-20", "20+"), include.lowest = TRUE)) # Sodium data_complet <- data_complet %>% mutate(classe_sodium = cut(sodium, breaks = c(0, 5, 10, 15, 20, 25, 30, Inf), labels = c("0-5", "5-10", "10-15", "15-20", "20-25", "25-30", "30+"), include.lowest = TRUE)) # Calcul des effectifs effectif_cd1 <- data_complet %>% group_by(classe_calcium) %>% summarise(n = n(), y_pos = quantile(cadmium, 0, na.rm = TRUE), .groups = "drop") effectif_cd2 <- data_complet %>% group_by(classe_magnesium) %>% summarise(n = n(), y_pos = quantile(cadmium, 0, na.rm = TRUE), .groups = "drop") effectif_cd3 <- data_complet %>% group_by(classe_sodium) %>% summarise(n = n(), y_pos = quantile(cadmium, 0, na.rm = TRUE), .groups = "drop") effectif_ni1 <- data_complet %>% group_by(classe_calcium) %>% summarise(n = n(), y_pos = quantile(nickel, 0, na.rm = TRUE), .groups = "drop") effectif_ni2 <- data_complet %>% group_by(classe_magnesium) %>% summarise(n = n(), y_pos = quantile(nickel, 0, na.rm = TRUE), .groups = "drop") effectif_ni3 <- data_complet %>% group_by(classe_sodium) %>% summarise(n = n(), y_pos = quantile(nickel, 0, na.rm = TRUE), .groups = "drop") effectif_pb <- data_complet %>% group_by(classe_calcium) %>% summarise(n = n(), y_pos = quantile(plomb, 0, na.rm = TRUE), .groups = "drop") # Plomb degrade <- colorRampPalette(brewer.pal(9, "YlGnBu"))(15) ggplot(data_complet %>% filter(!is.na(classe_calcium), !is.na(plomb)), aes(x = classe_calcium, y = plomb, fill = classe_calcium)) + geom_boxplot(alpha = 0.7, outlier.colour = "grey1") + labs(x = "Classe de Calcium (mg/L)", y = "Plomb (mg/kg PF)", fill = "Classe Ca") + geom_text(data = effectif_pb, aes(x = classe_calcium, y = 0, label = paste0("n=", n)), vjust = -0.5, size = 4.5, fontface = "bold") + theme_minimal() + ylim(c(0,0.80)) + scale_fill_manual(values = degrade) + theme(axis.title = element_text(size = 14), axis.text = element_text(size = 12), legend.title = element_text(size = 14, face = "bold"), legend.text = element_text(size = 14)) # Nickel en fonction du Calcium degrade <- colorRampPalette(brewer.pal(9, "PuRd"))(15) ggplot(data_complet %>% filter(!is.na(classe_calcium), !is.na(nickel)), aes(x = classe_calcium, y = nickel, fill = classe_calcium)) + geom_boxplot(alpha = 0.7, outlier.colour = "grey1") + labs(x = "Classe de Calcium (mg/L)", y = "Nickel (mg/kg PF)", fill = "Classe Ca") + geom_text(data = effectif_ni1, aes(x = classe_calcium, y = 0, label = paste0("n=", n)), vjust = -0.5, size = 4.5, fontface = "bold") + theme_minimal() + ylim(c(0,0.75)) + scale_fill_manual(values = degrade) + theme(axis.title = element_text(size = 14), axis.text = element_text(size = 12), legend.title = element_text(size = 14, face = "bold"), legend.text = element_text(size = 14)) # Nickel en fonction du Magnésium degrade <- colorRampPalette(brewer.pal(9, "PuRd"))(6) ggplot(data_complet %>% filter(!is.na(classe_magnesium), !is.na(nickel)), aes(x = classe_magnesium, y = nickel, fill = classe_magnesium)) + geom_boxplot(alpha = 0.7, outlier.colour = "grey1") + labs(x = "Classe de Magnésium (mg/L)", y = "Nickel (mg/kg PF)", fill = "Classe Mg") + geom_text(data = effectif_ni2, aes(x = classe_magnesium, y = 0, label = paste0("n=", n)), vjust = -0.5, size = 4.5, fontface = "bold") + theme_minimal() + ylim(c(0,0.75)) + scale_fill_manual(values = degrade) + theme(axis.title = element_text(size = 14), axis.text = element_text(size = 12), legend.title = element_text(size = 14, face = "bold"), legend.text = element_text(size = 14)) # Nickel en fonction du Sodium degrade <- colorRampPalette(brewer.pal(9, "PuRd"))(7) ggplot(data_complet %>% filter(!is.na(classe_sodium), !is.na(nickel)), aes(x = classe_sodium, y = nickel, fill = classe_sodium)) + geom_boxplot(alpha = 0.7, outlier.colour = "grey1") + labs(x = "Classe de Sodium (mg/L)", y = "Nickel (mg/kg PF)", fill = "Classe Na") + geom_text(data = effectif_ni3, aes(x = classe_sodium, y = 0, label = paste0("n=", n)), vjust = -0.5, size = 4.5, fontface = "bold") + theme_minimal() + ylim(c(0,0.75)) + scale_fill_manual(values = degrade) + scale_fill_manual(values = degrade) + theme(axis.title = element_text(size = 14), axis.text = element_text(size = 12), legend.title = element_text(size = 14, face = "bold"), legend.text = element_text(size = 14)) # Cadmium en fonction du calcium degrade <- colorRampPalette(brewer.pal(9, "YlOrBr"))(15) ggplot(data_complet %>% filter(!is.na(classe_calcium), !is.na(cadmium)), aes(x = classe_calcium, y = cadmium, fill = classe_calcium)) + geom_boxplot(alpha = 0.7, outlier.colour = "grey1") + labs(x = "Classe de Calcium (mg/L)", y = "Cadmium (mg/kg PF)", fill = "Classe Ca") + geom_text(data = effectif_cd1, aes(x = classe_calcium, y = 0, label = paste0("n=", n)), vjust = -0.5, size = 4.5, fontface = "bold") + theme_minimal() + scale_fill_manual(values = degrade) + ylim(c(0, 0.15)) + scale_fill_manual(values = degrade) + theme(axis.title = element_text(size = 14), axis.text = element_text(size = 12), legend.title = element_text(size = 14, face = "bold"), legend.text = element_text(size = 14)) # Cadmium en fonction du magnésium degrade <- colorRampPalette(brewer.pal(9, "YlOrBr"))(6) ggplot(data_complet %>% filter(!is.na(classe_magnesium), !is.na(cadmium)), aes(x = classe_magnesium, y = cadmium, fill = classe_magnesium)) + geom_boxplot(alpha = 0.7, outlier.colour = "grey1") + labs(x = "Classe de Magnesium (mg/L)", y = "Cadmium (mg/kg PF)", fill = "Classe Mg") + geom_text(data = effectif_cd2, aes(x = classe_magnesium, y = 0, label = paste0("n=", n)), vjust = -0.5, size = 4.5, fontface = "bold") + theme_minimal() + scale_fill_manual(values = degrade) + ylim(c(0, 0.15)) + scale_fill_manual(values = degrade) + theme(axis.title = element_text(size = 14), axis.text = element_text(size = 12), legend.title = element_text(size = 14, face = "bold"), legend.text = element_text(size = 14)) # Cadmium en fonction du sodium degrade <- colorRampPalette(brewer.pal(9, "YlOrBr"))(7) ggplot(data_complet %>% filter(!is.na(classe_sodium), !is.na(cadmium)), aes(x = classe_sodium, y = cadmium, fill = classe_sodium)) + geom_boxplot(alpha = 0.7, outlier.colour = "grey1") + labs(x = "Classe de Sodium (mg/L)", y = "Cadmium (mg/kg PF)", fill = "Classe Na") + geom_text(data = effectif_cd3, aes(x = classe_sodium, y = 0, label = paste0("n=", n)), vjust = -0.5, size = 4.5, fontface = "bold") + theme_minimal() + scale_fill_manual(values = degrade) + ylim(c(0, 0.15)) + scale_fill_manual(values = degrade) + theme(axis.title = element_text(size = 14), axis.text = element_text(size = 12), legend.title = element_text(size = 14, face = "bold"), legend.text = element_text(size = 14)) #************************************************************************************************ #### Modèle BLM pour retrouver Ceau à partir des données terrain #************************************************************************************************ source("Documents/script_BLM_memoire.R") # comme ça je peux utiliser le modèle BLM sans avoir à le réécrire ici #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Paramètres t_expo <- 7 # durée d'exposition en jours #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Appliquer le modèle BLM inverse pour chaque concentration de Ca/Mg/Na # conversion en mol pour le modèle Ca_mol_stat = data_complet$calcium / 40.078 / 1000 Mg_mol_stat = data_complet$magnesium / 24.305 / 1000 Na_mol_stat = data_complet$sodium / 22.990 / 1000 H_mol_stat = 10^(-data_complet$pH) # Cadmium Ceau_Cd_mol = pmap_dbl( list(ci = data_complet$cadmium, ca = Ca_mol_stat, mg = Mg_mol_stat, na = Na_mol_stat, hh = H_mol_stat), function(ci, ca, mg, na, hh) { if (is.na(ci) | is.na(ca) | is.na(hh)) return(NA_real_) calculer_Ceau(ci, params_Cd, ca, mg, na, hh, t_expo) } ) # reconversion en µg Ceau_Cd_ugL = Ceau_Cd_mol * 112.41 * 1e6 # Nickel Ceau_Ni_mol = pmap_dbl( list(ci = data_complet$nickel, ca = Ca_mol_stat, mg = Mg_mol_stat, na = Na_mol_stat, hh = H_mol_stat), function(ci, ca, mg, na, hh) { if (is.na(ci) | is.na(ca) | is.na(hh)) return(NA_real_) calculer_Ceau(ci, params_Ni, ca, mg, na, hh, t_expo) } ) Ceau_Ni_ugL = Ceau_Ni_mol * 58.693 * 1e6 # Plomb Ceau_Pb_mol = pmap_dbl( list(ci = data_complet$plomb, ca = Ca_mol_stat, mg = Mg_mol_stat, na = Na_mol_stat, hh = H_mol_stat), function(ci, ca, mg, na, hh) { if (is.na(ci) | is.na(ca) | is.na(hh)) return(NA_real_) calculer_Ceau(ci, params_Pb, ca, mg, na, hh, t_expo) } ) Ceau_Pb_ugL = Ceau_Pb_mol * 207.2 * 1e6 # Rajouter les colonnes dans data_complet data_complet <- data_complet %>% mutate( Ceau_Cd_ugL = Ceau_Cd_mol * 112.41 * 1e6, Ceau_Ni_ugL = Ceau_Ni_mol * 58.693 * 1e6, Ceau_Pb_ugL = Ceau_Pb_mol * 207.2 * 1e6) #************************************************************************************************ #### Récupération sur Naïades des concentrations métalliques dans l'eau ------------------------ #************************************************************************************************ # Codes SANDRE des métaux dans l'eau codes_metaux <- c( 1388, # Cadmium 1386, # Nickel 1382) # Plomb # Récupération par groupes de 10 stations groupes <- split(station, ceiling(seq_along(station) / 10)) metaux_list <- list() for (i in seq_along(groupes)) { groupe <- groupes[[i]] stations_str <- paste(groupe, collapse = ",") url <- paste0( "https://hubeau.eaufrance.fr/api/v2/qualite_rivieres/analyse_pc?", "code_station=", stations_str, "&code_parametre=", paste(codes_metaux, collapse = ","), "&code_support=3", "&fraction=3", "&date_debut_prelevement=2023-01-01", "&date_fin_prelevement=2024-12-31", "&size=20000", "&format=json" ) reponse <- GET(url) if (status_code(reponse) %in% c(200, 206)) { contenu <- fromJSON(content(reponse, as = "text", encoding = "UTF-8")) if (length(contenu$data) > 0) { metaux_list[[i]] <- contenu$data } } } # Combiner metaux_eau_raw <- bind_rows(metaux_list) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Nettoyage de metaux_eau_raw metaux_eau <- metaux_eau_raw %>% select(code_stat = code_station, date = date_prelevement, metal = libelle_parametre, concentration = resultat, unite = symbole_unite, code_remarque) %>% mutate( date = as.Date(date), mois = lubridate::month(date, label = TRUE, abbr = TRUE), annee = lubridate::year(date), # code_remarque = 10 veut dire "< seuil de quantification" # on les garde mais on les signale sous_LQ = code_remarque == 10) %>% filter(metal %in% c("Cadmium", "Nickel", "Plomb")) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Statistiques descriptives stats_eau <- metaux_eau %>% group_by(metal, unite) %>% summarise(n = n(), n_sous_LQ = sum(sous_LQ, na.rm = TRUE), pct_sous_LQ = round(100 * n_sous_LQ / n, 1), min = min(concentration, na.rm = TRUE), q25 = quantile(concentration, 0.25, na.rm = TRUE), mediane = median(concentration, na.rm = TRUE), moyenne = mean(concentration, na.rm = TRUE), q75 = quantile(concentration, 0.75, na.rm = TRUE), max = max(concentration, na.rm = TRUE), .groups = "drop") print(stats_eau) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Graphique 1 Distribution des concentrations (boxplot) NQE <- data.frame(metal = c("Nickel", "Plomb", "Cadmium"), nqe = c(4.0, 1.2, 0.08)) ggplot(metaux_eau, aes(x = metal, y = concentration, fill = metal)) + geom_boxplot(alpha = 0.7, outlier.colour = "grey1") + stat_summary(fun = median, geom = "text", aes(label = sprintf("%.3f", after_stat(y))), hjust = -3, vjust = -0.5, size = 5, color = "darkred", fontface = "bold") + geom_point(data = NQE, aes(x = metal, y = nqe), shape = 4, color = "red", size = 6, stroke = 2, inherit.aes = FALSE) + scale_y_log10() + labs(x = "Metal", y = "Concentration (µg/L) (échelle log)") + theme_minimal() + scale_fill_manual(values = c("Plomb" = "#21908CFF", "Nickel" = "#B63679FF", "Cadmium" = "#F98C0AFF")) + theme(legend.position = "none", plot.title = element_text(face = "bold", hjust = 0.5), axis.title = element_text(size = 16), axis.text = element_text(size = 16)) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Graphique 2 Evolution temporelle avec NQE NQE <- data.frame(metal = c("Cadmium", "Nickel", "Plomb"), nqe = c(0.08, 4.0, 1.2), # µg/L NQE MA nqe_label = c("NQE Cd = 0.08", "NQE Ni = 4.0", "NQE Pb = 1.2")) ggplot(metaux_eau, aes(x = date, y = concentration, color = metal)) + geom_point(alpha = 0.4, size = 1.5) + geom_smooth(method = "loess", se = TRUE, linewidth = 1) + geom_hline(data = NQE, aes(yintercept = nqe), linetype = "dashed", color = "red", linewidth = 0.8) + geom_text(data = NQE, aes(x = as.Date("2023-02-01"), y = nqe, label = nqe_label), color = "red", vjust = -0.5, size = 3) + facet_wrap(~ metal, scales = "free_y") + scale_color_manual(values = c("Cadmium" = "#F98C0AFF", "Nickel" = "#B63679FF", "Plomb" = "#21908CFF")) + labs(x = "Date", y = "Concentration (µg/L)") + theme_minimal() + theme(legend.position = "none", plot.title = element_text(face = "bold", hjust = 0.5)) + theme(axis.title = element_text(size = 16), axis.text = element_text(size = 16), legend.title = element_text(size = 16, face = "bold"), legend.text = element_text(size = 16)) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Graphique 3 Variabilité mensuelle (saisonnalité) ggplot(metaux_eau, aes(x = mois, y = concentration, fill = metal)) + geom_boxplot(alpha = 0.7) + facet_wrap(~ metal, scales = "free_y") + scale_fill_manual(values = c("Cadmium" = "#F98C0AFF", "Nickel" = "#B63679FF", "Plomb" = "#21908CFF")) + labs(x = "Mois", y = "Concentration (µg/L) (échelle log)") + scale_y_log10() + theme_minimal() + theme(axis.title = element_text(size = 16), axis.text = element_text(size = 12), legend.title = element_text(size = 16, face = "bold"), legend.text = element_text(size = 14)) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Médiane par station et par métal (avec filtre LQ) # LQ : Cadmium = 0.01 µg/L, Nickel = 0.5 µg/L, Plomb = 0.1 µg/L # On remplace par NA les valeurs sous LQ SANS supprimer le reste de la ligne metaux_eau_med <- metaux_eau %>% mutate(concentration = case_when( metal == "Cadmium" & concentration <= 0.01 ~ NA_real_, metal == "Nickel" & concentration <= 0.5 ~ NA_real_, metal == "Plomb" & concentration <= 0.1 ~ NA_real_, TRUE ~ concentration)) %>% group_by(code_stat, metal) %>% summarise(Ceau_med_ugL = median(concentration, na.rm = TRUE), n = n(), .groups = "drop") %>% pivot_wider(names_from = metal, values_from = c(Ceau_med_ugL, n)) head(metaux_eau_med) # Jointure avec data_complet data_complet <- data_complet %>% left_join(metaux_eau_med %>% select(code_stat, Ceau_med_ugL_Cadmium, Ceau_med_ugL_Nickel,Ceau_med_ugL_Plomb), by = "code_stat") #************************************************************************************************ #### Comparaison Ceau prédite (BLM) vs Ceau mesurée (Naïades) ---------------------------------- #************************************************************************************************ # Graphique mesurée vs prédite: même logique que le QQ plot # Cadmium data_complet %>% filter(!is.na(Ceau_Cd_ugL), !is.na(Ceau_med_ugL_Cadmium)) %>% ggplot(aes(x = Ceau_med_ugL_Cadmium, y = Ceau_Cd_ugL)) + geom_point(color = "#F98C0AFF", size = 3, alpha = 0.7) + geom_abline(slope = 1, intercept = 0, # droite 1:1 linetype = "dashed", color = "grey40") + labs(x = "Ceau mesurée Naïades (µg/L)", y = "Ceau prédite BLM (µg/L)") + ylim(c(0,0.25)) + xlim(c(0,0.5)) + theme_minimal() + theme(plot.title = element_text(face = "bold", hjust = 0.5)) + theme(axis.title = element_text(size = 16), axis.text = element_text(size = 16), legend.title = element_text(size = 16, face = "bold"), legend.text = element_text(size = 16)) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Graphique densité Ceau prédite (BLM) vs réelle avec NQE nqe <- data.frame(metal = c("Cadmium", "Nickel", "Plomb"), NQE = c(0.08, 4.0, 1.2)) data_long_Ceau <- data_complet %>% filter(!is.na(Ceau_Cd_ugL), !is.na(Ceau_Ni_ugL), !is.na(Ceau_Pb_ugL)) %>% select(Ceau_Cd_ugL, Ceau_Ni_ugL, Ceau_Pb_ugL, Ceau_med_ugL_Cadmium, Ceau_med_ugL_Nickel, Ceau_med_ugL_Plomb) %>% # Remplacer par NA les valeurs sous seuil LQ sans supprimer toute la ligne mutate( Ceau_med_ugL_Cadmium = ifelse(Ceau_med_ugL_Cadmium <= 0.01, NA, Ceau_med_ugL_Cadmium), Ceau_med_ugL_Nickel = ifelse(Ceau_med_ugL_Nickel <= 0.5, NA, Ceau_med_ugL_Nickel), Ceau_med_ugL_Plomb = ifelse(Ceau_med_ugL_Plomb <= 0.1, NA, Ceau_med_ugL_Plomb)) %>% pivot_longer( cols = everything(), names_to = "variable", values_to = "valeur") %>% mutate(metal = case_when( grepl("Cd|Cadmium", variable) ~ "Cadmium", grepl("Ni|Nickel", variable) ~ "Nickel", grepl("Pb|Plomb", variable) ~ "Plomb"), source = case_when(grepl("med", variable) ~ "Réelle", TRUE ~ "Prédite (BLM)")) %>% filter(!is.na(valeur), !is.infinite(valeur), valeur > 0) ggplot(data_long_Ceau, aes(x = valeur, fill = source, color = source)) + geom_density(alpha = 0.4, linewidth = 0.8) + geom_vline(data = nqe, aes(xintercept = NQE), linetype = "dashed", color = "red", linewidth = 1, inherit.aes = FALSE) + geom_text(data = nqe, aes(x = NQE, y = Inf, label = paste0("NQE=", NQE, " µg/L")), color = "red", angle = 90, vjust = -0.3, hjust = 1.2, size = 5, inherit.aes = FALSE) + facet_wrap(~ metal, scales = "free") + scale_fill_manual(values = c("Réelle" = "grey50", "Prédite (BLM)" = "#440154FF")) + scale_color_manual(values = c("Réelle" = "grey50", "Prédite (BLM)" = "#440154FF")) + scale_x_log10(labels = scales::label_comma(accuracy = 0.01)) + labs(x = "Concentration dans l'eau (µg/L) (échelle log)", y = "Densité", fill = "Source", color = "Source") + theme_minimal() + theme(strip.text = element_text(size = 18, face = "bold"), axis.title = element_text(size = 16), axis.text = element_text(size = 16), legend.title = element_text(size = 16, face = "bold"), legend.text = element_text(size = 16)) #************************************************************************************************ #### Modèle TK (comparaison avec TK-BLM) ------------------------------------------------------- #************************************************************************************************ # Paramètres TK params_Ni_TK <- list( ku = 0.33, # L/g/j ke = 0.13, # j-1 Ceau = 3.78) # µg/L (0.00378 mg/L x 1000) params_Pb_TK <- list( ku = 1.01, # L/g/j ke = 0.33, # j-1 Ceau = 10.0) # µg/L (0.01 mg/L x 1000) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Résolution du modèle TK modele_TK <- function(t, state, params) { Ci <- state["Ci"] dCi <- params$ku * params$Ceau - params$ke * Ci return(list(dCi)) } sol_Ni_TK <- ode(y = c(Ci = 0.1892), times = temps, func = modele_TK, parms = params_Ni_TK) sol_Pb_TK <- ode(y = c(Ci = 0), times = temps, func = modele_TK, parms = params_Pb_TK) tab_Ni_TK <- as.data.frame(sol_Ni_TK) %>% mutate(metal = "Nickel (TK)", modele = "TK") tab_Pb_TK <- as.data.frame(sol_Pb_TK) %>% mutate(metal = "Plomb (TK)", modele = "TK") tab_Ni_BLM <- tab_Ni %>% mutate(modele = "TK-BLM") tab_Pb_BLM <- tab_Pb %>% mutate(modele = "TK-BLM") #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Inversion TK Ceau = f(Ci) calculer_Ceau_TK <- function(Ci_t, params, t, Ci0 = 0) { Ceau_ugL <- (Ci_t - Ci0 * exp(-params$ke * t)) * params$ke / (params$ku * (1 - exp(-params$ke * t))) return(Ceau_ugL) } Ci_gamme <- seq(0.01, 5, by = 0.01) # µg/g PS # Calcul Ceau pour chaque Ci Ceau_TK_Ni <- sapply(Ci_gamme, function(ci) calculer_Ceau_TK(ci, params_Ni_TK, t_expo, Ci0 = 0.1892)) Ceau_TK_Pb <- sapply(Ci_gamme, function(ci) calculer_Ceau_TK(ci, params_Pb_TK, t_expo, Ci0 = 0)) tab_inv_Ni_TK <- data.frame(Ci_ugkg = Ci_gamme, Ceau_ugL = Ceau_TK_Ni, metal = "Nickel (TK)") tab_inv_Pb_TK <- data.frame(Ci_ugkg = Ci_gamme, Ceau_ugL = Ceau_TK_Pb, metal = "Plomb (TK)") #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Appliquer le modèle TK inverse sur les données terrain data_complet <- data_complet %>% mutate( # Nickel TK Ceau_Ni_TK_ugL = (nickel - 0.1892 * exp(-params_Ni_TK$ke * t_expo)) * params_Ni_TK$ke / (params_Ni_TK$ku * (1 - exp(-params_Ni_TK$ke * t_expo))), # Plomb TK Ceau_Pb_TK_ugL = (plomb - 0 * exp(-params_Pb_TK$ke * t_expo)) * params_Pb_TK$ke / (params_Pb_TK$ku * (1 - exp(-params_Pb_TK$ke * t_expo)))) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Graphique densité TK vs réelle avec NQE data_long_TK <- data_complet %>% select(Ceau_Ni_TK_ugL, Ceau_Pb_TK_ugL, Ceau_med_ugL_Nickel, Ceau_med_ugL_Plomb) %>% mutate( Ceau_med_ugL_Nickel = ifelse(Ceau_med_ugL_Nickel < 0.5, NA, Ceau_med_ugL_Nickel), Ceau_med_ugL_Plomb = ifelse(Ceau_med_ugL_Plomb < 0.1, NA, Ceau_med_ugL_Plomb) ) %>% pivot_longer(cols = everything(), names_to = "variable", values_to = "valeur") %>% mutate(metal = case_when( grepl("Ni|Nickel", variable) ~ "Nickel", grepl("Pb|Plomb", variable) ~ "Plomb"), source = case_when( grepl("med", variable) ~ "Réelle", TRUE ~ "Prédite (TK)")) %>% filter(!is.na(valeur), !is.infinite(valeur), valeur > 0) ggplot(data_long_TK, aes(x = valeur, fill = source, color = source)) + geom_density(alpha = 0.4, linewidth = 0.8) + geom_vline(data = nqe %>% filter(metal != "Cadmium"), aes(xintercept = NQE), linetype = "dashed", color = "red", linewidth = 1, inherit.aes = FALSE) + geom_text(data = nqe %>% filter(metal != "Cadmium"), aes(x = NQE, y = Inf, label = paste0("NQE=", NQE, " µg/L")), color = "red", angle = 90, vjust = -0.3, hjust = 1.2, size = 5, inherit.aes = FALSE) + facet_wrap(~ metal, scales = "free") + scale_fill_manual(values = c("Réelle" = "grey60", "Prédite (TK)" = "#440154FF")) + scale_color_manual(values = c("Réelle" = "grey60", "Prédite (TK)" = "#440154FF")) + scale_x_log10() + labs(x = "Concentration dans l'eau (µg/L) (échelle log)", y = "Densité", fill = "Source", color = "Source") + theme_minimal() + theme(strip.text = element_text(size = 18, face = "bold"), axis.title = element_text(size = 16), axis.text = element_text(size = 16), legend.title = element_text(size = 16, face = "bold"), legend.text = element_text(size = 16))