diff --git a/R/calcMonthlyMean.R b/R/calcMonthlyMean.R index f43b1369a5add617698db7963ede557b0e8ff1a6..2410a2ed0f413aae4a23dc6a6c1df59f5cb92d5b 100644 --- a/R/calcMonthlyMean.R +++ b/R/calcMonthlyMean.R @@ -42,7 +42,8 @@ calcMonthlyMeanQuantile <- function(x, prob, threshold = 0.8) { dfMY[, 1] <- lubridate::month(dfMY[, 1]) l <- lapply(seq.int(12), function(mo) { - apply(dfMY[dfMY[, 1] == mo, -1], 2, quantile, probs = prob, type = 8) + apply(dfMY[dfMY[, 1] == mo, -1], 2, + quantile, probs = prob, type = 8, na.rm = TRUE) }) dfMMQ <- do.call(rbind, l) rownames(dfMMQ) <- paste0("M", seq.int(12)) diff --git a/R/calcQMNA.R b/R/calcQMNA.R index bde27d547f1bfdb830078f335e81d4552cba0238..60be68764ef1afbb347f254064e4d835226019d9 100644 --- a/R/calcQMNA.R +++ b/R/calcQMNA.R @@ -199,6 +199,7 @@ calcQLogNn <- function(x, return_period) { x <- x[!is.na(x)] nbY <- length(x) + if (nbY == 0) return(NA) Freq <- 1 / return_period nbxnul <- length(x[x == 0]) diff --git a/R/convertUnitQ.R b/R/convertUnitQ.R index 709e84b51c5e53a7f1cb5f6a2a7af0a4833820cf..381dd99d2298d9fcf6ba8d43018bbe84a4e154a2 100644 --- a/R/convertUnitQ.R +++ b/R/convertUnitQ.R @@ -56,5 +56,5 @@ convertUnitQ.data.frame <- function(x, areas, from, to, ...) { out <- matrix(out, ncol = 1) colnames(out) <- ids } - cbind(x[, 1], as.data.frame(out)) + cbind(DatesR = x[, 1], as.data.frame(out)) } diff --git a/R/saveFlowDB.R b/R/saveFlowDB.R index a0766fa10dbce6dcfe4309f7c984fa93239e3751..398ad2e00b68d2724f6f9c0b0316981fae47b25d 100644 --- a/R/saveFlowDB.R +++ b/R/saveFlowDB.R @@ -20,34 +20,40 @@ saveFlowDB <- function(path, InputsModel, OutputsModel, indicators_periods = NUL names(BS_reseau) <- c("CdSiteHydro", "lambert2.x", "lambert2.y", "area", "nom", "id_aval", "distance_aval", "model") readr::write_tsv(BS_reseau, file.path(path, "stations.tsv")) - # Simulated flows - Qsim <- attr(OutputsModel, "Qm3s") - Qsim$DatesR <- format(Qsim$DatesR, format = "%Y-%m-%d") - readr::write_tsv(Qsim, file.path(path, "ts_Qsim.tsv") - ) + if (inherits(OutputsModel, "GRiwrmOutputsModel")) { + # Simulated flows + Qsim <- attr(OutputsModel, "Qm3s") + Qsim$DatesR <- format(Qsim$DatesR, format = "%Y-%m-%d") + readr::write_tsv(Qsim, file.path(path, "ts_Qsim.tsv")) - # Simulated contributing flows (GR outputs) - listQGR <- lapply(names(OutputsModel), function(id) { - if (is.null(OutputsModel[[id]]$QsimDown)) { - return(OutputsModel[[id]]$Qsim_m3) - } else { - areas <- InputsModel[[id]]$BasinAreas - return(OutputsModel[[id]]$QsimDown * areas[length(areas)] * 1e3) - } - }) # m3/time-step - Qcontrib <- do.call(cbind, listQGR) / attr(InputsModel, "TimeStep") # m3/s - colnames(Qcontrib) <- names(OutputsModel) - Qcontrib <- cbind(DatesR = Qsim$DatesR, as.data.frame(Qcontrib)) + # Simulated contributing flows (GR outputs) + listQGR <- lapply(names(OutputsModel), function(id) { + if (is.null(OutputsModel[[id]]$QsimDown)) { + return(OutputsModel[[id]]$Qsim_m3) + } else { + areas <- InputsModel[[id]]$BasinAreas + return(OutputsModel[[id]]$QsimDown * areas[length(areas)] * 1e3) + } + }) # m3/time-step + Qcontrib <- do.call(cbind, listQGR) / attr(InputsModel, "TimeStep") # m3/s + colnames(Qcontrib) <- names(OutputsModel) + Qcontrib <- cbind(DatesR = Qsim$DatesR, as.data.frame(Qcontrib)) - readr::write_tsv(Qcontrib, file.path(path, "ts_Qcontrib.tsv")) + readr::write_tsv(Qcontrib, file.path(path, "ts_Qcontrib.tsv")) + Q <- attr(OutputsModel, "Qm3s") + } else { + Qobs <- OutputsModel + Qobs$DatesR <- format(Qobs$DatesR, format = "%Y-%m-%d") + readr::write_tsv(Qobs, file.path(path, "ts_Qobs.tsv")) + Q <- OutputsModel + } # Global hydrological indicators if (!is.null(indicators_periods)) { lapply(indicators_periods, function(period) { # global indicators - Qsim <- attr(OutputsModel, "Qm3s") - Qsim <- Qsim[Qsim$DatesR >= period[1] & Qsim$DatesR <= period[2], ] - dfInd <- calcHydrologicalIndicators(Qsim) + Q <- Q[Q$DatesR >= period[1] & Q$DatesR <= period[2], ] + dfInd <- calcHydrologicalIndicators(Q) y <- format(period, "%Y") readr::write_tsv(dfInd, file.path(path, @@ -55,7 +61,7 @@ saveFlowDB <- function(path, InputsModel, OutputsModel, indicators_periods = NUL paste(y, collapse = "-"), ".tsv"))) # Monthly regime - QMIA <- calcMonthlyInterannualMean(Qsim) + QMIA <- calcMonthlyInterannualMean(Q) QMIA <- cbind(Id = colnames(QMIA), t(QMIA)) readr::write_tsv(as.data.frame(QMIA), file.path(path, @@ -63,7 +69,7 @@ saveFlowDB <- function(path, InputsModel, OutputsModel, indicators_periods = NUL paste(y, collapse = "-"), ".tsv"))) # Monthly 5 years dry regime - QM5<- calcMonthlyMeanQuantile(Qsim, 0.2) + QM5<- calcMonthlyMeanQuantile(Q, 0.2) QM5 <- cbind(Id = colnames(QM5), t(QM5)) readr::write_tsv(as.data.frame(QM5), file.path(path, diff --git a/bookdown/02-calage.Rmd b/bookdown/02-calage.Rmd index b35cf47fb155156d42e1ba6fd0b8ea77694a1ad5..ad2ad362fc1fca788f2d9fda94e9c52499591764 100644 --- a/bookdown/02-calage.Rmd +++ b/bookdown/02-calage.Rmd @@ -154,7 +154,7 @@ plotRatio(Qsim, Qobs, calcQJXA10, "R-QJXA10") ## Débits caractéristiques -Nous calculons ici l'ensemble des indicateurs hydrologiques utilisés dans les fiches Explore 2070: +Nous calculons ici l'ensemble des indicateurs hydrologiques utilisés dans les fiches Explore 2070 et le rapport entre les valeurs observées et calculées : ```{r, eval=cfg$data$write_results} indicators <- c("QA", @@ -175,7 +175,7 @@ calcQcaract <- function(indicator, Qobs, Qsim) { areas <- griwrm$area names(areas) <- griwrm$id Qobs_m3s <- convertUnitQ(Qobs, areas, from = "mm", to = "m3/s") -Qsim_m3s <- convertUnitQ(Qsim, areas, from = "mm", to = "m3/s") +Qsim_m3s <- attr(OM_with, "Qm3s") l <- lapply(indicators, calcQcaract, Qobs = Qobs_m3s, Qsim = Qsim_m3s) mQc <- do.call(cbind, l) dfQc <- cbind(data.frame(Id = rownames(mQc)), mQc) @@ -186,5 +186,27 @@ write.table(dfQc, row.names = FALSE) ``` +On enregistre ici les indicateurs hydrologiques pour les débits observés et simulés au même format que celui qui sera utilisé pour analyser les débits naturels historiques et projetés. + +```{r} +indicators_periods <- list( + calage = as.POSIXct(c(cfg$calibration$date$start, cfg$calibration$date$end), tz = "UTC"), + reference = as.POSIXct(cfg$hydroclim$drias$periods$ref, tz = "UTC") +) + +InputsModel <- CreateInputsModel(BasinsObs_with) +saveFlowDB(path = getDataPath(cfg$calibration$path, "Qobs", cfg = cfg), + InputsModel = InputsModel, + OutputsModel = Qobs_m3s, + indicators_periods = indicators_periods, + cfg = cfg) +saveFlowDB(path = getDataPath(cfg$calibration$path, "Qsim", cfg = cfg), + InputsModel = InputsModel, + OutputsModel = OM_with, + indicators_periods = indicators_periods, + cfg = cfg) +``` + + Ces indicateurs sont téléchargeables à l'adresse : https://nextcloud.inrae.fr/s/adinzGa3AmLEXnZ/download?path=%2F02-calibration&files=Q_indicators.tsv