feat(calage): compute all hydrologicial indicators

- fix errors due to data gaps

Refs #12

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))

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])

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))
 }

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,

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&nbsp;:
 https://nextcloud.inrae.fr/s/adinzGa3AmLEXnZ/download?path=%2F02-calibration&files=Q_indicators.tsv