Commit a3f3a1ab authored by Dorchies David's avatar Dorchies David
Browse files

feat(evolution): add delta calculation

- add pmedian (from pmin and pmax)

Refs #12
parent a33aede5
Pipeline #31923 passed with stages
in 54 minutes and 9 seconds
......@@ -44,6 +44,7 @@ export(max_na_rm)
export(mean_na_rm)
export(min_na_rm)
export(not_na_count)
export(pmedian)
export(readDrias2020)
export(readQsim)
export(saveFlowDB)
......
......@@ -15,7 +15,6 @@
#' @source Chazot, Sébastien, Perrin, Charles, Jean-Philippe Vidal, Eric Sauquet, Mathilde Chauveau, et Nathalie Rouchy. 2012. « Explore 2070 - Lot Hydrologie de surface - A2 - Résultats : Fiches, Cartes et Graphes ». Ministère de l’écologie, du développement durable, des transports et du logement.
#'
#' @examples
calcQA <- function(x, threshold = 0.8, ...) {
UseMethod("calcQA", x)
}
......
#' Median value
#'
#' @description
#' Returns the (regular or parallel) median of the input values.
#'
#' `pmedian` takes one or more vectors as arguments which should have same dimensions.
#'
#' @param ... Several [matrix] of same dimensions
#'
#' @return a [matrix] of same dimension as the ones in `...` with median value calculated for each cell.
#' @export
#'
pmedian <- function(...) {
l <- list(...)
m <- sapply(seq.int(length(l[[1]])), function(i) {
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]])
}
return(m)
}
......@@ -23,7 +23,7 @@ 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 <- cfg$hydroclim$drias$rcp[-1]
rcps <- "rcp4.5"
scenarios <- gsub("/", "_", cfg$hydroclim$drias$scenarios)
scenarios <- rep(scenarios, length(rcps))
......@@ -34,6 +34,19 @@ loadAllBasinsObs <- function(rcp, scenario) {
}
AllBasinsObs <- mapply(loadAllBasinsObs, rcp = rcps, scenario = scenarios, SIMPLIFY = FALSE)
names(AllBasinsObs) <- paste(rcps, scenarios, sep = " - ")
```
On réarrange les données pour avoir une liste de data.frame pour chaque variable étudiée:
```{r}
formatObs <- function(BasinObs, item) {
cbind(DatesR = BasinObs$DatesR, as.data.frame(BasinObs[[item]]))
}
P_drias <- lapply(AllBasinsObs, formatObs, item = "P")
E_drias <- lapply(AllBasinsObs, formatObs, item = "E")
T_drias <- lapply(AllBasinsObs, formatObs, item = "Temp")
```
Chaque élément de la liste `AllBasinsObs` est une liste structurée de classe *BasinsObs* (Voir `?createBasinsObs`).
......@@ -74,7 +87,25 @@ Pour les débits&nbsp;:
- Débit moyen de période de retour 5 ans
```{r}
monthlyMean <- function(df) {
dfOut <- SeriesAggreg(df, "%m", rep("mean", ncol(df) - 1))
cbind(month = lubridate::month(dfOut$DatesR), as.matrix(dfOut[, -1]))
}
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, ])
}
monthlyMean(selectPeriod(P_drias[[1]][, 1:3], "ref"))
```
......@@ -117,11 +148,66 @@ 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}
calcComparison <- function(obs, sim, scenarios, FUN, ...) {
calcDelta <- function(sim, FUN, period, delta = "*", notFirstCol = FALSE, ...) {
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
} else {
m <- vFut[[scenario]] - vRef[[scenario]]
}
if (notFirstCol) {
m[, 1] <- vRef[[scenario]][, 1]
}
return(m)
})
list(
min = do.call(pmin, deltas),
med = do.call(pmedian, deltas),
max = do.call(pmax, deltas)
)
}
```
Exemple pour les pluies moyennes mensuelles:
```{r}
deltaP <- calcDelta(P_drias, monthlyMean, period = "fin", notFirstCol = TRUE)
lapply(deltaP, function(x) x[, 1:3])
```
Exemple pour la température moyenne mensuelle:
```{r}
deltaT <- calcDelta(T_drias, monthlyMean, period = "fin", delta = "+", notFirstCol = TRUE)
lapply(deltaT, function(x) x[, 1:3])
```
**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:
```{r}
deltaQMNA5 <- calcDelta(QsimDrias, calcQMNAn, period = "fin", return_period = 5)
lapply(deltaQMNA5, function(x) x[1:2])
```
Exemple sur le VCN10-2:
```{r}
deltaVCN10_2 <- calcDelta(QsimDrias, calcVCNn, period = "fin", k = 10, return_period = 2)
lapply(deltaVCN10_2, function(x) x[1:2])
```
Exemple sur le QJXA10:
```{r}
deltaQJXA10 <- calcDelta(QsimDrias, calcQJXAn, period = "fin", return_period = 2)
lapply(deltaQJXA10, function(x) x[1:2])
```
## Lecture des données
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/pmedian.R
\name{pmedian}
\alias{pmedian}
\title{Median value}
\usage{
pmedian(...)
}
\arguments{
\item{...}{Several \link{matrix} of same dimensions}
}
\value{
a \link{matrix} of same dimension as the ones in \code{...} with median value calculated for each cell.
}
\description{
Returns the (regular or parallel) median of the input values.
\code{pmedian} takes one or more vectors as arguments which should have same dimensions.
}
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