feat: calcul des périodes de retour de crue à Paris

DESCRIPTION

@@ -16,7 +16,6 @@ Encoding: UTF-8
 LazyData: true
 Suggests:
     bookdown,
-    ggplot2,
     knitr,
     lubridate,
     prettymapr,
@@ -28,19 +27,21 @@ Suggests:
     testthat (>= 3.0.0),
     tidyr,
     mapsf,
-    tidyquant
+    tidyquant,
+    extRemes
 VignetteBuilder: knitr
 Depends:
     R (>= 2.10),
     airGRiwrm (>= 0.5.0.9002)
 Imports:
+    ggplot2,
     config,
     glue,
     httr,
     magrittr,
     ncdf4,
     purrr,
-    readr, 
+    readr,
     sf,
     stats,
     utils,

bookdown/07-caracterisation_crues_paris.Rmd

@@ -4,7 +4,7 @@
 library(seinebasin2)
 library(tidyr)
 library(ggplot2)
-library(InraeThemes)
+library(extRemes)
 cfg <- loadConfig()
 knitr::opts_chunk$set(
   fig.width = 8,
@@ -86,7 +86,90 @@ plotReturnPeriods(QJXA_end, "Période 2071-2100")
 
 ```{r}
 period <- c(cfg$hydroclim$drias$periods$ref[1], cfg$hydroclim$drias$periods$end[2])
-plotReturnPeriods(statQJXA(QsimDrias, period), 
+QJXA_all <- statQJXA(QsimDrias, period)
+plotReturnPeriods(QJXA_all, 
                   "Période 1976-2100")
 ```
 
+## Ajustement des lois de probabilité
+
+```{r}
+return.level.ci <- function(fit, t_max, alpha = 0.05) {
+  return.period <- seq(2, t_max)
+  estimate <- return.level(fit, return.period = return.period)
+  estimate[estimate < 0] <- 0
+  
+  l <- lapply(alpha, function(a) {
+    rl <- return.level(fit, return.period = return.period, do.ci = TRUE, alpha = a)
+    class(rl) <- "matrix"
+    rl <- rl[,-2]
+    rl[rl<0] <- 0
+    colnames(rl) <- c(sprintf("p%02.0f", a/2 * 100), sprintf("p%02.0f", (1-a/2) * 100))
+    rl
+  })
+  rl <- do.call(cbind, l)
+  rl <- cbind(estimate, rl)
+  rl <- rl[,sort(colnames(rl))]
+  rl <- cbind(return.period = return.period, rl)
+  as.data.frame(rl)
+}
+
+bilan_scenario <- function(scenario, QJXA, t_max = 100, alphas) {
+  cat("\n\n### Scénario ", scenario, "\n")
+  QJXA <- QJXA[QJXA$scenario == scenario, ]
+  fit <- fevd(x = Q, data = QJXA, type = "Gumbel", time.units = "years")
+  dfRl <- return.level.ci(fit, t_max, alphas)
+  p <- ggplot(dfRl, 
+         aes(x = return.period)) + 
+    geom_ribbon(aes(ymin = p02, ymax = p98, fill = "Incertitude 95%"), fill = "#66c1bf", alpha = 0.3, show.legend = T) +
+    geom_ribbon(aes(ymin = p15, ymax = p85), fill = "#008c8e", alpha = 0.3) +
+    geom_line(aes(y = estimate, color = "Estimation"), color = "#423089") + 
+    geom_point(data = QJXA, 
+               aes(x = t, y = Q, color = "Observations"), color = "#ed6e6c") +
+    labs(x = "Période de retour (années)",
+         y = "Débit (m3/s)",
+         title = paste("Scénario", scenario))
+  print(p)
+}
+```
+
+Il faut aussi choisir des niveaux d'incertitudes, nous prendrons les valeurs classiques de 70% et 95% d'incertitude.
+
+```{r}
+alphas <- c(0.05, 0.3)
+```
+
+
+```{r}
+bilan_scenario(unique(QJXA_ref$scenario)[1], QJXA_ref, alphas = alphas)
+```
+
+## Ajustements lois de Gumbel période 1976-2005
+
+```{r, results = 'asis'}
+invisible(lapply(unique(QJXA_ref$scenario), bilan_scenario, QJXA = QJXA_ref, return.period = return.period, alphas = alphas))
+```
+## Ajustements lois de Gumbel période 2021-2050
+
+```{r, results = 'asis'}
+invisible(lapply(unique(QJXA_near$scenario), bilan_scenario, QJXA = QJXA_near, return.period = return.period, alphas = alphas))
+```
+
+## Ajustements lois de Gumbel période 2041-2070
+
+```{r, results = 'asis'}
+invisible(lapply(unique(QJXA_middle$scenario), bilan_scenario, QJXA = QJXA_middle, return.period = return.period, alphas = alphas))
+```
+
+## Ajustements lois de Gumbel période 2071-2100
+
+```{r, results = 'asis'}
+invisible(lapply(unique(QJXA_end$scenario), bilan_scenario, QJXA = QJXA_end, return.period = return.period, alphas = alphas))
+```
+
+## Ajustements lois de Gumbel période 1976-2100
+
+```{r, results = 'asis'}
+invisible(lapply(unique(QJXA_all$scenario), bilan_scenario, QJXA = QJXA_all, return.period = c(return.period, 150), alphas = alphas))
+```
+