Commit 386521f6 authored by Dorchies David's avatar Dorchies David
Browse files

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

parent d1c688b1
Pipeline #33634 passed with stages
in 47 minutes and 14 seconds
......@@ -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,
......
......@@ -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,85 @@ 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)
```
## Ajustements lois de Gumbel période 1976-2005
```{r, results = 'asis'}
invisible(lapply(unique(QJXA_ref$scenario), bilan_scenario, QJXA = QJXA_ref, alphas = alphas))
```
## Ajustements lois de Gumbel période 2021-2050
```{r, results = 'asis'}
invisible(lapply(unique(QJXA_near$scenario), bilan_scenario, QJXA = QJXA_near, alphas = alphas))
```
## Ajustements lois de Gumbel période 2041-2070
```{r, results = 'asis'}
invisible(lapply(unique(QJXA_middle$scenario), bilan_scenario, QJXA = QJXA_middle, alphas = alphas))
```
## Ajustements lois de Gumbel période 2071-2100
```{r, results = 'asis'}
invisible(lapply(unique(QJXA_end$scenario), bilan_scenario, QJXA = QJXA_end, alphas = alphas))
```
## Ajustements lois de Gumbel période 1976-2100
```{r, results = 'asis'}
invisible(lapply(unique(QJXA_all$scenario), bilan_scenario, QJXA = QJXA_all, t_max = 150, alphas = alphas))
```
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment