Extraction and analysis of data station by station
Objectives
In the tibble created in the previous example, a lot of stations are not suitable for computing basic statistical parameters (annual or interannual mean, median, standard deviation, etc.) because they provide only few results each year and are not representative of the seasonal cycle of nitrates. This means that the tibble resulting from the query may include a lot of unexploitable data and could be much smaller.
The data can be computed directly from the database station by station by including the query in a loop using a known list of stations. For each station, the number of results (available nitrates data) is calculated each year and compared to a threshold to determine if a statistical/graphical analysis can be computed or not.
Selection of stations available for analysis
In the following example, a list of the stations of the Côte d’Or department is created.
Total number of stations is:
We use the function get_qualite_rivieres_analyse() to
retrieve nitrate concentration values
(code_parametre = "1340") from 2000 to 2022.
nitrates_21 <- get_qualite_rivieres_analyse(
code_departement = "21",
date_debut_prelevement = "2000-01-01",
date_fin_prelevement = "2022-12-31",
code_parametre = "1340",
fields = c(
"code_station",
"libelle_station",
"libelle_fraction",
"date_prelevement",
"resultat",
"symbole_unite"
)
)We compute some annual statistics for each station:
nitrates_21 <- nitrates_21 %>%
mutate(date_prelevement = as.POSIXct(date_prelevement),
year = year(date_prelevement))
station_stats <- nitrates_21 %>%
group_by(code_station, libelle_station, year) %>%
summarise(nb_analyses = n(),
nitrate_mean = mean(resultat),
nitrate_p90 = quantile(resultat, probs = 0.9),
.groups = 'drop')
station_stats
#> # A tibble: 1,830 × 6
#> code_station libelle_station year nb_analyses nitrate_mean nitrate_p90
#> <chr> <chr> <dbl> <int> <dbl> <dbl>
#> 1 03000029 LE RUISSEAU DU FEU A… 2017 6 23.4 37.8
#> 2 03000029 LE RUISSEAU DU FEU A… 2021 6 28.3 36.5
#> 3 03000029 LE RUISSEAU DU FEU A… 2022 12 24.5 36.6
#> 4 03000033 La Seine à Courceau 2016 6 21.3 29.2
#> 5 03000033 La Seine à Courceau 2017 6 18.2 30.2
#> 6 03000033 La Seine à Courceau 2021 6 25.8 31.5
#> 7 03000033 La Seine à Courceau 2022 18 20 25
#> 8 03000070 Ruisseau de Jugny 2016 6 15.7 20.4
#> 9 03000070 Ruisseau de Jugny 2017 5 20.9 35.8
#> 10 03000070 Ruisseau de Jugny 2021 6 17.8 25.5
#> # ℹ 1,820 more rowsThe number of lines of the tibble station_stats
corresponds to the number of stations with at least one sample.
Stations with less than 10 values per year on average and less than 10 years of data are excluded of the statistical analysis.
valid_stations <- station_stats %>%
group_by(code_station, libelle_station) %>%
summarise(analyses_per_year = mean(nb_analyses), nb_years = n()) %>%
filter(analyses_per_year >= 10, nb_years >= 10)
valid_stations
#> # A tibble: 13 × 4
#> # Groups: code_station [13]
#> code_station libelle_station analyses_per_year nb_years
#> <chr> <chr> <dbl> <int>
#> 1 03001000 "LA SEINE A NOD-SUR-SEINE 1" 15.4 22
#> 2 03002100 "LA SEINE A CHARREY-SUR-SEINE 1" 11.4 22
#> 3 03006268 "L'OURCE A MENESBLE 1" 15.8 16
#> 4 03006271 "LA GROÈME A TERREFONDREE 1" 12.8 18
#> 5 03006590 "L'OURCE A AUTRICOURT 1" 13.9 22
#> 6 03034720 "LE SEREIN A VIEUX-CHATEAU 2" 14.9 18
#> 7 03036510 "L'ARMANÇON A MONTIGNY-SUR-ARMANCON … 15.9 16
#> 8 03040250 "L'OZE A DARCEY 1" 11.2 17
#> 9 06011000 "SAONE A AUXONNE 1" 13.2 22
#> 10 06014250 "CENT FONDS A SAULON-LA-RUE 1" 10.9 10
#> 11 06016500 "OUCHE A ECHENON " 11.2 22
#> 12 06017050 "SAONE A CHARREY-SUR-SAONE " 12.9 22
#> 13 06017070 "SAONE A SEURRE " 12.8 16The number of rows of the tibble valid_stations
corresponds to the number of stations with at least 10 samples per year
on average and 10 years of data.
Statistical analysis of samples
Then, we plot the annual distribution of nitrate levels with violin plots.
plot_nitrates <- function(code) {
station_details <- station_21[station_21$code_station == code, , drop = FALSE]
mean_samples <- valid_stations$analyses_per_year[valid_stations$code_station == code]
nitrates_station <- nitrates_21 %>% filter(code_station == code)
station_yearly_stats <- station_stats %>% filter(code_station == code)
p <- ggplot(nitrates_station, aes(x = as.factor(year), y = resultat)) +
labs(
x = "year",
y = "nitrates (mg/l)",
title = paste(
"station:",
station_details$code_station,
station_details$libelle_station
),
subtitle = paste(round(mean_samples, 1), "samples per year on average"),
caption = "mean and sd in blue, median + 1rst and 3rd quartile represented by dashed lines, percentile 90 in red"
) +
scale_x_discrete(labels = paste0(station_yearly_stats$year, "\nn=", station_yearly_stats$nb_analyses))
p <-
p +
geom_violin(
trim = TRUE,
scale = "width",
adjust = 0.5,
draw_quantiles = 0.9,
color = "red",
fill = "lightblue1"
) + # draw the violin and adds an horizontal red line corresponding to the quantile 90
geom_violin(
trim = TRUE,
scale = "width",
adjust = 0.5,
color = "black",
fill = "transparent"
) + # draw the same violin but with black lines and no fill
geom_violin(
trim = TRUE,
scale = "width",
adjust = 0.5,
draw_quantiles = c(0.25, 0.5, 0.75),
linetype = "dashed",
fill = "transparent"
) + # adds the median, the 1st and 3rd quartiles in dashed line
stat_summary(
fun.data = mean_sdl,
fun.args = list(mult = 1),
geom = "pointrange",
color = "blue4",
fill = "transparent"
) # adds the mean and the standard deviation in blue
}
