Skip to content

GitLab

Commit 2a9a711b authored by Dorchies David's avatar Dorchies David
Browse files

feat(delta): calc all indicators on drias2020 

Refs #12
parent 5cd2a6ad
Pipeline #32039 failed with stages
in 31 minutes and 33 seconds
Hide whitespace changes
Inline Side-by-side
NAMESPACE
View file @ 2a9a711b
......@@ -23,6 +23,7 @@ export(addReservoirsGRiwrm)
export(addReservoirsQ)
export(areaPoly)
export(calcMonthlyInterannualMean)
export(calcMonthlyInterannualSum)
export(calcMonthlyMean)
export(calcMonthlyMeanQuantile)
export(calcQA)
......
R/calcMonthlyMean.R
View file @ 2a9a711b
......@@ -3,6 +3,8 @@
#' @description
#' Compute montly means, interannual monthly means and quantile of interannual monthly means.
#'
#' Compute also the mean of monthly cumulated data for precipitations and ETP
#'
#' @param x a [data.frame] with a first column containing dates in [POSIXt] format and one [numeric] column per gauging station
#' @inheritParams calcQMNA
#' @param prob [numeric], non-exceedance probability
......@@ -16,8 +18,8 @@ calcMonthlyInterannualMean <- function(x, threshold = 0.8) {
dfMY <- calcMonthlyMean(x, threshold)
dfM <- SeriesAggreg(dfMY, Format = "%m", rep("mean_na_rm", ncol(x)-1))
dfM2 <- as.matrix(dfM[, -1, drop = FALSE])
rownames(dfM2) <- lubridate::month(dfM[, 1])
dfM2[as.character(seq.int(12)), ]
rownames(dfM2) <- paste0("M", lubridate::month(dfM[, 1]))
dfM2[paste0("M", seq.int(12)), ]
}
#' @rdname calcMonthlyMean
......@@ -43,6 +45,16 @@ calcMonthlyMeanQuantile <- function(x, prob, threshold = 0.8) {
apply(dfMY[dfMY[, 1] == mo, -1], 2, quantile, probs = prob, type = 8)
})
dfMMQ <- do.call(rbind, l)
rownames(dfMMQ) <- seq.int(12)
rownames(dfMMQ) <- paste0("M", seq.int(12))
dfMMQ
}
#' @rdname calcMonthlyMean
#' @export
calcMonthlyInterannualSum <- function(x) {
PM <- SeriesAggreg(x, Format = "%Y%m", rep("sum", ncol(x) - 1))
dfM <- SeriesAggreg(PM, Format = "%m", rep("mean", ncol(x) - 1))
dfM2 <- as.matrix(dfM[, -1, drop = FALSE])
rownames(dfM2) <- paste0("M", lubridate::month(dfM[, 1]))
dfM2[paste0("M", seq.int(12)), ]
}
R/plot_seine_map.R
View file @ 2a9a711b
......@@ -9,7 +9,6 @@
#' @return
#' @export
#'
#' @examples
plot_seine_map <- function(r, breaks, title, credits = "IGN, Inrae", ...) {
library(mapsf)
......
R/pmedian.R
View file @ 2a9a711b
......@@ -16,11 +16,8 @@ pmedian <- function(...) {
cells <- sapply(l, function(x) x[i])
median(cells)
})
if (is.matrix(l[[1]])) {
m <- matrix(m, ncol = ncol(l[[1]]))
colnames(m) <- colnames(l[[1]])
} else {
names(m) <- names(l[[1]])
}
m <- matrix(m, ncol = ncol(l[[1]]))
colnames(m) <- colnames(l[[1]])
rownames(m) <- rownames(l[[1]])
return(m)
}
bookdown/06-evolution_indicateurs_hydrologiques.Rmd
View file @ 2a9a711b
......@@ -9,9 +9,30 @@ L'évolution est représentée sous la forme du rapport entre l'indicateur calcu
La période de référence correspond à la période 1976-2005 et la période future choisie ici correspond à la période "fin de siècle" 2071-2100.
## Lecture des séries temporelles
```{r}
selectPeriod <- function(df, period) {
periods <- lapply(cfg$hydroclim$drias$periods, as.POSIXct, tz = "UTC")
if (!period %in% names(periods)) {
stop("`period`should be in ",
paste(sprintf("\"%s\"", names(periods)), collapse = ", "))
}
selectDates <- df$DatesR >= periods[[period]][1] &
df$DatesR <= periods[[period]][2]
return(df[selectDates, ])
}
```
### Données climatiques
Les indicateurs sont calculés séparément pour les deux scénarios d'émission RCP 4.5 et RCP 8.5 et une synthèse des différents couples GCM/RCM est réalisée.
```{r}
rcps <- cfg$hydroclim$drias$rcp[-1]
scenarios <- gsub("/", "_", cfg$hydroclim$drias$scenarios)
scenariosX <- rep(scenarios, length(rcps))
rcpsX <- rep(rcps, each = length(scenarios))
```
## Lecture des données climatiques
Les données climatiques sont stockées dans les objets `BasinsObs` (Voir `?loadBasinsObs` pour le chargement). Le nom des fichiers contenant ces données est formé comme suit&nbsp;: `BasinObs_[RCP]_[GCM_RCM].RDS`. Ce qui se traduit en R par&nbsp:
......@@ -23,18 +44,16 @@ BasinsObs <- loadBasinsObs(file, cfg = cfg)
On peut charger toutes les données climatiques de tous les scénarios GCM/RCM dans une liste avec le script suivant&nbsp;:
```{r}
rcps <- "rcp4.5"
scenarios <- gsub("/", "_", cfg$hydroclim$drias$scenarios)
scenarios <- rep(scenarios, length(rcps))
loadAllBasinsObs <- function(rcp, scenario) {
message("Processing ", rcp, " scenario ", scenario, "...")
file <- paste0(paste("BasinObs", rcp, scenario, sep = "_"), ".RDS")
loadBasinsObs(file, cfg = cfg)
}
AllBasinsObs <- mapply(loadAllBasinsObs, rcp = rcps, scenario = scenarios, SIMPLIFY = FALSE)
names(AllBasinsObs) <- paste(rcps, scenarios, sep = " - ")
AllBasinsObs <- mapply(loadAllBasinsObs,
rcp = rcpsX,
scenario = scenariosX,
SIMPLIFY = FALSE)
names(AllBasinsObs) <- paste(rcpsX, scenariosX, sep = " - ")
```
On réarrange les données pour avoir une liste de data.frame pour chaque variable étudiée:
......@@ -47,68 +66,93 @@ formatObs <- function(BasinObs, item) {
P_drias <- lapply(AllBasinsObs, formatObs, item = "P")
E_drias <- lapply(AllBasinsObs, formatObs, item = "E")
T_drias <- lapply(AllBasinsObs, formatObs, item = "Temp")
rm(AllBasinsObs)
```
Chaque élément de la liste `AllBasinsObs` est une liste structurée de classe *BasinsObs* (Voir `?createBasinsObs`).
### Débits simulés
Les données sont stockées au format TSV (i.e. texte avec séparateur tabulation) dans le dossier&nbsp;: `cfg$Qnat$path/Drias2020/Qnat-v1/[RCP]/[GCM_RCM]/`
On le récupère avec la commande suivante&nbsp;:
```r
path <- getDataPath(cfg$Qnat$path, "Drias2020/Qnat-v1", rcp, scenario, cfg = cfg)
```
Les débits simulés aux différents noeuds du modèle semi-distribué peuvent être chargés à l'aide de la fonction suivante&nbsp;:
```{r}
loadDriasQnat <- function(rcp_scenario) {
message("Read Qsim for ", rcp_scenario)
readQsim(cfg$Qnat$path, "Drias2020/Qnat-v1", rcp_scenario, cfg = cfg)
}
QsimDrias <- lapply(paste(rcps, scenarios, sep="/"), loadDriasQnat)
names(QsimDrias) <- paste(rcps, scenarios, sep = " - ")
```
## Listes des indicateurs
### Données moyennes mensuelles et annuelles
### Indicateurs mensuels
Pour le climat&nbsp;:
Pour les données climatiques&nbsp;:
- Cumul pluviométrie (mm)
- Moyenne températures (°C)
- Cumul ETP (mm)
```{r}
# Précipitation moyenne mensuelle pour la période de référence
calcMonthlyInterannualSum(selectPeriod(P_drias[[1]][, 1:3], "ref"))
# Température moyenne mensuelle pour l'horizon fin de siècle
calcMonthlyInterannualMean(selectPeriod(T_drias[[1]][, 1:3], "end"))
# Calcul pour tous les scénarios et périodes
calcPeriod <- function(period, df, calcFUN) {
ind <- calcFUN(selectPeriod(df, period))
t(ind)
}
calcAll <- function(data, calcFUN) {
lapply(data, function(df) {
periods <- names(cfg$hydroclim$drias$periods)
names(periods) <- periods
lapply(as.list(periods), calcPeriod, df = df, calcFUN = calcFUN)
})
}
P_month <- calcAll(P_drias, calcMonthlyInterannualSum)
T_month <- calcAll(T_drias, calcMonthlyInterannualMean)
E_month <- calcAll(E_drias, calcMonthlyInterannualSum)
```
Pour les débits&nbsp;:
- Débit moyen
- Débit moyen sec de période de retour 5 ans
Ces données ont été calculées et enregistrées lors de la simulation des débits.
```{r}
monthlyMean <- function(df) {
dfOut <- SeriesAggreg(df, "%m", rep("mean", ncol(df) - 1))
cbind(month = lubridate::month(dfOut$DatesR), as.matrix(dfOut[, -1]))
loadIndicators <- function(rcp, scenario, indicator) {
periods <- names(cfg$hydroclim$drias$periods)
names(periods) <- periods
lapply(periods, function(period) {
file <- paste0(indicator, "_",
substr(cfg$hydroclim$drias$periods[[period]][1], 1, 4),
"-",
substr(cfg$hydroclim$drias$periods[[period]][2], 1, 4),
".tsv")
file <- getDataPath(
cfg$Qnat$path,
"Drias2020/Qnat-v1",
rcp,
scenario,
file,
cfg = cfg
)
df <- read.csv(file, sep = "\t")
m <- as.matrix(df[, -1])
rownames(m) <- df[, 1]
m
})
}
selectPeriod <- function(df, period) {
mPeriods <- data.frame(period = c("ref", "proche", "milieu", "fin"),
start = as.POSIXct(c("1976-01-01", "2021-01-01", "2041-01-01", "2071-01-01"), tz = "UTC"),
end = as.POSIXct(c("2005-12-31", "2050-12-31", "2070-12-31", "2100-12-31"), tz = "UTC"))
period <- period
if (!period %in% mPeriods$period) {
stop("`period`should be equal to ", paste(sprintf("\"%s\"", mPeriods$period), collapse = ", "))
}
selectDates <- df$DatesR >= mPeriods$start[mPeriods$period == period] &
df$DatesR <= mPeriods$end[mPeriods$period == period]
return(df[selectDates, ])
loadAllIndicators <- function(indicator) {
ind <- mapply(rcp = rcpsX,
scenario = scenariosX,
indicator = indicator,
loadIndicators,
SIMPLIFY = FALSE)
names(ind) <- paste(rcpsX, scenariosX, sep = " - ")
return(ind)
}
monthlyMean(selectPeriod(P_drias[[1]][, 1:3], "ref"))
Q_month <- loadAllIndicators("Q_monthly")
Q_month5 <- loadAllIndicators("Q_monthly_5years")
```
### Indicateurs synthétiques
Pour les étiages&nbsp;: VCN10, VCN30 et QMNA pour les périodes de retour 2 ans, 5 ans et 10 ans.
......@@ -120,20 +164,7 @@ Pour les faibles/forts débits&nbsp;: les quantiles 95% et 10% de débits journa
Exemples de calculs pour les débits simulés pour le scénario `r names(QsimDrias)[1]` des stations La Seine à Vernon et la Seine à Paris Austerlitz sur la période de référence 1976-2005&nbsp;:
```{r}
Qsim <- QsimDrias[[1]] # Sélection du scénario
# Sélection de la plage de date et des colonnes date et 1ère station
selectDates <- Qsim$DatesR >= as.POSIXct("1976-01-01") &
Qsim$DatesR <= as.POSIXct("2005-12-31")
QsimRef <- Qsim[selectDates, names(Qsim) %in% c("DatesR", "H7900010", "H5920010")]
# QMNA5
calcQMNAn(QsimRef, 5)
# QJXA 10
calcQJXAn(QsimRef, 10)
# VCN10-2
calcVCNn(QsimRef, 10, 2)
Q_indicators <- loadAllIndicators("Q_indicators")
```
......@@ -148,22 +179,27 @@ Les données utilisées en entrées sont&nbsp;:
Les tableaux de synthèse fournissent la valeur minimale, médiane et maximale des évolutions parmi les couples scénario/modèle climatiques.
```{r}
calcDelta <- function(sim, FUN, period, delta = "*", notFirstCol = FALSE, ...) {
calcDelta <- function(ind, rcp, period, delta = "*") {
if (!delta %in% c("*", "+")) stop("`delta`should be equal to \"*\" or \"+\"")
simRef <- lapply(sim, selectPeriod, period = "ref")
vRef <- lapply(simRef, FUN, ...)
simFut <- lapply(sim, selectPeriod, period = period)
vFut <- lapply(simFut, FUN, ...)
deltas <- lapply(names(sim), function(scenario) {
if (delta == "*") {
m <- (vFut[[scenario]] - vRef[[scenario]]) / vRef[[scenario]] * 100
scenarios <- names(ind)
names(scenarios) <- scenarios
ind <- lapply(scenarios, function(x) {
if (grepl(rcp, x, fixed = TRUE)) {
ind[[x]]
} else {
m <- vFut[[scenario]] - vRef[[scenario]]
NULL
}
if (notFirstCol) {
m[, 1] <- vRef[[scenario]][, 1]
})
ind <- ind[!sapply(ind,is.null)]
deltas <- lapply(names(ind), function(scenario) {
dfPer <- ind[[scenario]][[period]]
dfRef <- ind[[scenario]]$ref
if (delta == "*") {
m <- (dfPer - dfRef) / dfRef * 100
} else {
m <- dfPer - dfRef
}
return(m)
})
list(
min = do.call(pmin, deltas),
......@@ -171,36 +207,42 @@ calcDelta <- function(sim, FUN, period, delta = "*", notFirstCol = FALSE, ...) {
max = do.call(pmax, deltas)
)
}
tableDeltaStation <- function(station, delta) {
l <- lapply(delta, function(m) {
m[station, ]
})
do.call(cbind, l)
}
```
Exemple pour les pluies moyennes mensuelles:
Exemple pour les pluies moyennes mensuelles pour les 6 premières stations:
```{r}
deltaP <- calcDelta(P_drias, monthlyMean, period = "fin", notFirstCol = TRUE)
lapply(deltaP, function(x) x[, 1:3])
deltaPM <- calcDelta(P_month, rcp = "rcp4.5", period = "end")
lapply(deltaPM, function(x) t(head(x)))
```
Exemple pour la température moyenne mensuelle:
Exemple pour la température moyenne mensuelle à Paris:
```{r}
deltaT <- calcDelta(T_drias, monthlyMean, period = "fin", delta = "+", notFirstCol = TRUE)
lapply(deltaT, function(x) x[, 1:3])
deltaTM <- calcDelta(T_month, rcp = "rcp4.5", period = "end", delta = "+")
knitr::kable(tableDeltaStation("H8012010", deltaTM), digits = 1)
```
**N.B.:** les données climatiques correspondent ici au données moyennes du sous-bassin versant et pas le bassin versant entier. Il faudrait agréger les données des bassins amont pour avoir une moyenne du bassin pour chaque station.
Exemple sur le QMNA5:
Exemple sur le débit mensuel moyen à Paris
```{r}
deltaQMNA5 <- calcDelta(QsimDrias, calcQMNAn, period = "fin", return_period = 5)
lapply(deltaQMNA5, function(x) x[1:2])
deltaQM <- calcDelta(Q_month, rcp = "rcp4.5", period = "end")
knitr::kable(tableDeltaStation("H8012010", deltaQM), digits = 1)
```
Exemple sur le VCN10-2:
Exemple sur tous les indicateurs hydrologiques à Paris:
```{r}
deltaVCN10_2 <- calcDelta(QsimDrias, calcVCNn, period = "fin", k = 10, return_period = 2)
lapply(deltaVCN10_2, function(x) x[1:2])
deltaIndicators <- calcDelta(Q_indicators, rcp = "rcp4.5", period = "end")
knitr::kable(tableDeltaStation("H8012010", deltaIndicators), digits = 1)
```
Exemple sur le QJXA10:
......@@ -210,4 +252,3 @@ deltaQJXA10 <- calcDelta(QsimDrias, calcQJXAn, period = "fin", return_period = 2
lapply(deltaQJXA10, function(x) x[1:2])
```
man/calcMonthlyMean.Rd
View file @ 2a9a711b
......@@ -4,6 +4,7 @@
\alias{calcMonthlyInterannualMean}
\alias{calcMonthlyMean}
\alias{calcMonthlyMeanQuantile}
\alias{calcMonthlyInterannualSum}
\title{Compute monthly mean}
\usage{
calcMonthlyInterannualMean(x, threshold = 0.8)
......@@ -11,6 +12,8 @@ calcMonthlyInterannualMean(x, threshold = 0.8)
calcMonthlyMean(x, threshold = 0.8)
calcMonthlyMeanQuantile(x, prob, threshold = 0.8)
calcMonthlyInterannualSum(x)
}
\arguments{
\item{x}{a \link{data.frame} with a first column containing dates in \link{POSIXt} format and one \link{numeric} column per gauging station}
......@@ -24,4 +27,6 @@ calcMonthlyMeanQuantile(x, prob, threshold = 0.8)
}
\description{
Compute montly means, interannual monthly means and quantile of interannual monthly means.
Compute also the mean of monthly cumulated data for precipitations and ETP
}
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