fix(calage): hydrologic indicators in m3/s instead of mm/time step

Refs #5

NAMESPACE

@@ -18,6 +18,7 @@ S3method(calcQMNAn,data.frame)
 S3method(calcVCN,matrix)
 S3method(calcVCN,numeric)
 S3method(calcVCNn,data.frame)
+S3method(convertUnitQ,data.frame)
 export(addReservoirsGRiwrm)
 export(addReservoirsQ)
 export(areaPoly)
@@ -41,6 +42,7 @@ export(calcVCN30_10)
 export(calcVCN30_2)
 export(calcVCN30_5)
 export(calcVCNn)
+export(convertUnitQ)
 export(createBasinsObs)
 export(getAprioriIds)
 export(getDataPath)

R/convertUnitQ.R

+#' Convert flows units
+#'
+#' @param x a [data.frame] with a first column with the time stamp in [POSIXct] format and folling columns with flows
+#' @param ...
+#'
+#' @return [data.frame]  with a first column with the time stamp in [POSIXct] format and folling columns with converted flows
+#'
+#' @rdname ConvertUnitQ
+#' @export
+#'
+#' @examples
+#' data("L0123001", package = "airGR")
+#' dfQmm <- data.frame(DatesR = BasinObs$DatesR, L0123001 = BasinObs$Qmm)
+#' dfQls <- convertUnitQ.data.frame(dfQmm, c(L0123001 = BasinInfo$BasinArea), from = "mm", to = "l/s")
+#' head(dfQls)
+#' head(BasinObs$Qls)
+convertUnitQ <- function(x, ...) {
+    UseMethod("convertUnitQ", x)
+}
+
+
+#' @rdname ConvertUnitQ
+#' @param areas named [vector], basin areas in km², `names` must correspond to station ids in columns names of `x`
+#' @param from flow unit of `x`, see details
+#' @param to flow unit of returned values
+#'
+#' @details
+#' Available units are:
+#'
+#' - "mm" for mm per time step
+#' - "m3" for m3 per time step
+#' - "m3/s" for m3 per second
+#' - "l/s" for liter per second
+#'
+#' The time step is determined with the time difference between rows in `x`.
+#'
+#' @export
+convertUnitQ.data.frame <- function(x, areas, from, to, ...) {
+    ids <- names(x)[-1]
+    areas <- areas[ids]
+    time_step <- as.numeric(x[2,1] - x[1,1], units = "secs")
+    # Remove "/" from unit
+    from <- gsub("/", "", from)
+    to <- gsub("/", "", to)
+    # Factor to convert from m3/time step to...
+    factors <- list(mm = 1E-3 / areas,
+                    m3 = 1,
+                    m3s = 1 / time_step,
+                    ls = 1000 / time_step)
+    out <- apply(x[, -1, drop = FALSE], 1, function(row) {
+        row  * factors[[to]] / factors[[from]]
+    })
+    if (is.matrix(out)) {
+        out <- t(out)
+    } else {
+        out <- matrix(out, ncol = 1)
+        colnames(out) <- ids
+    }
+    cbind(x[, 1], as.data.frame(out))
+}

bookdown/02-calage.Rmd

@@ -4,7 +4,8 @@
 library(seinebasin2)
 cfg <- loadConfig()
 knitr::opts_chunk$set(
-  fig.width = 8
+  fig.width = 8, 
+  fig.asp = 1
 )
 ```
 
@@ -28,11 +29,13 @@ Le calage est alors effectué en prenant en compte les options par défaut qui s
 - utiliser la régularisation des paramètres en utilisant le bassin amont le plus contributif comme donneur de paramètres a priori.
 
 ```{r, eval=cfg$calibration$eval, warning=FALSE, message=FALSE}
-OC <-
-  Calibration(
-    BasinsObs,
-    paramFile = "calibration_without_reservoirs.txt",
-  )
+if (cfg$calibration$eval) {
+    OC <-
+      Calibration(
+        BasinsObs,
+        paramFile = "calibration_without_reservoirs.txt",
+      )
+}
 ```
 
 ## Calage des stations du bassin avec prise en compte des réservoirs {#calage-avec-reservoirs}
@@ -50,18 +53,19 @@ BasinsObs_with$Q <- addReservoirsQ(BasinsObs_with)
 L'opération de calage peut ensuite avoir lieu avec prise en compte de ces influences.
 
 ```{r, eval=cfg$calibration$eval, warning=FALSE, message=FALSE}
-OC <-
-  Calibration(
-    BasinsObs_with,
-    paramFile = "calibration_with_reservoirs.txt",
-  )
+if (cfg$calibration$eval) {
+    OC <-
+      Calibration(
+        BasinsObs_with,
+        paramFile = "calibration_with_reservoirs.txt",
+      )
+}
 ```
 
 ## Comparaison des débits simulés avec et sans prise en compte des réservoirs
 
 Il faut lancer la simulation des deux modèles pour les sortir cote à cote. La fonction `seinebasin2::RunModel` ne nécessite qu'on lui fournisse les données d'entrée et le paramétrage du modèle pour fournir les chroniques de débits simulés.
 
-
 ```{r, message=FALSE}
 OM_without <- RunModel(BasinsObs, "calibration_without_reservoirs.txt")
 OM_with <- RunModel(BasinsObs_with, "calibration_with_reservoirs.txt")
@@ -119,7 +123,6 @@ plotRatio <- function(Qsim, Qobs, FUN, title) {
     breaks <- c(decQuant, incQuant)
     
     plot_seine_map(r, breaks, title)
-
 }
 ```
 
@@ -153,7 +156,7 @@ plotRatio(Qsim, Qobs, calcQJXA10, "R-QJXA10")
 
 Nous calculons ici l'ensemble des indicateurs hydrologiques utilisés dans les fiches Explore 2070:
 
-```{r}
+```{r, eval=cfg$data$write_results}
 indicators <- c("QA", 
                 "VCN10_2", "VCN30_2", "QMNA2",
                 "VCN10_5", "VCN30_5", "QMNA5",
@@ -169,7 +172,11 @@ calcQcaract <- function(indicator, Qobs, Qsim) {
                      paste0("R-", indicator))
     return(m)
 }
-l <- lapply(indicators, calcQcaract, Qobs = Qobs, Qsim = 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")
+l <- lapply(indicators, calcQcaract, Qobs = Qobs_m3s, Qsim = Qsim_m3s)
 mQc <- do.call(cbind, l)
 dfQc <- cbind(data.frame(Id = rownames(mQc)), mQc)
 write.table(dfQc, 
@@ -179,5 +186,5 @@ write.table(dfQc,
           row.names = FALSE)
 ```
 
-Ces indicateurs sont téléchargeable à l'adresse&nbsp;:
+Ces indicateurs sont téléchargeables à l'adresse&nbsp;:
 https://owncloud.dorch.fr/index.php/s/pCPZvY4lk6xGC8m/download?path=%2F02-calibration&files=Q_indicators.tsv

man/ConvertUnitQ.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/convertUnitQ.R
+\name{convertUnitQ}
+\alias{convertUnitQ}
+\alias{convertUnitQ.data.frame}
+\title{Convert flows units}
+\usage{
+convertUnitQ(x, ...)
+
+\method{convertUnitQ}{data.frame}(x, areas, from, to, ...)
+}
+\arguments{
+\item{x}{a \link{data.frame} with a first column with the time stamp in \link{POSIXct} format and folling columns with flows}
+
+\item{...}{}
+
+\item{areas}{named \link{vector}, basin areas in km², \code{names} must correspond to station ids in columns names of \code{x}}
+
+\item{from}{flow unit of \code{x}, see details}
+
+\item{to}{flow unit of returned values}
+}
+\value{
+\link{data.frame}  with a first column with the time stamp in \link{POSIXct} format and folling columns with converted flows
+}
+\description{
+Convert flows units
+}
+\details{
+Available units are:
+\itemize{
+\item "mm" for mm per time step
+\item "m3" for m3 per time step
+\item "m3/s" for m3 per second
+\item "l/s" for liter per second
+}
+
+The time step is determined with the time difference between rows in \code{x}.
+}
+\examples{
+data("L0123001", package = "airGR")
+dfQmm <- data.frame(DatesR = BasinObs$DatesR, L0123001 = BasinObs$Qmm)
+dfQls <- convertUnitQ.data.frame(dfQmm, c(L0123001 = BasinInfo$BasinArea), from = "mm", to = "l/s")
+head(dfQls)
+head(BasinObs$Qls)
+}