From 799e719f84e43aa767493211ef60e3b6f2e698a6 Mon Sep 17 00:00:00 2001
From: Dorchies David <david.dorchies@inrae.fr>
Date: Wed, 3 Mar 2021 15:35:12 +0100
Subject: [PATCH] feat: add attribute Qm3s to GRiwrmOutputsModel

Closes #30
---
 R/CreateInputsModel.GRiwrm.R   | 26 +++++++++++++++++
 R/RunModel.GRiwrmInputsModel.R | 20 ++++++-------
 R/RunModel.Supervisor.R        |  1 +
 R/utils.R                      | 52 ++++++++++++++++++++++++++++++++++
 4 files changed, 89 insertions(+), 10 deletions(-)

diff --git a/R/CreateInputsModel.GRiwrm.R b/R/CreateInputsModel.GRiwrm.R
index 9f28b80..78287a8 100644
--- a/R/CreateInputsModel.GRiwrm.R
+++ b/R/CreateInputsModel.GRiwrm.R
@@ -20,6 +20,7 @@ CreateInputsModel.GRiwrm <- function(x, DatesR, Precip, PotEvap, Qobs, ...) {
       id, x, DatesR,Precip[,id], PotEvap[,id], Qobs, ...
     )
   }
+  attr(InputsModel, "TimeStep") <- getModelTimeStep(InputsModel)
   return(InputsModel)
 }
 
@@ -92,3 +93,28 @@ CreateOneGRiwrmInputsModel <- function(id, griwrm, DatesR, Precip, PotEvap, Qobs
 
   return(InputsModel)
 }
+
+
+#' Check time steps of the model of all the nodes and return the time step in seconds
+#'
+#' This function is called inside [CreateInputsModel.GRiwrm] for defining the time step of the big model.
+#'
+#' @param InputsModel a `GRiwrmInputsModel`
+#'
+#' @return A [numeric] representing the time step in seconds
+#'
+getModelTimeStep <- function(InputsModel) {
+  TS <- sapply(InputsModel, function(x) {
+    if (inherits(x, "hourly")) {
+      TimeStep <- 60 * 60
+    } else if (inherits(x, "daily")) {
+      TimeStep <- 60 * 60 * 24
+    } else {
+      stop("All models should be at hourly or daily time step")
+    }
+  })
+  if(length(unique(TS)) != 1) {
+    stop("Time steps of the model of all nodes should be identical")
+  }
+  return(unique(TS))
+}
diff --git a/R/RunModel.GRiwrmInputsModel.R b/R/RunModel.GRiwrmInputsModel.R
index ff6047f..c401406 100644
--- a/R/RunModel.GRiwrmInputsModel.R
+++ b/R/RunModel.GRiwrmInputsModel.R
@@ -14,21 +14,21 @@ RunModel.GRiwrmInputsModel <- function(x, RunOptions, Param, ...) {
   OutputsModel <- list()
   class(OutputsModel) <- c("GRiwrmOutputsModel", class(OutputsModel))
 
-  for(IM in x) {
-    message("RunModel.GRiwrmInputsModel: Treating sub-basin ", IM$id, "...")
+  for(id in names(x)) {
+    message("RunModel.GRiwrmInputsModel: Treating sub-basin ", x[[id]]$id, "...")
 
-    # Update x$Qupstream with simulated upstream flows
-    if(any(IM$UpstreamIsRunoff)) {
-      IM <- UpdateQsimUpstream(IM, RunOptions[[IM$id]]$IndPeriod_Run, OutputsModel)
+    # Update x[[id]]$Qupstream with simulated upstream flows
+    if(any(x[[id]]$UpstreamIsRunoff)) {
+      x[[id]] <- UpdateQsimUpstream(x[[id]], RunOptions[[id]]$IndPeriod_Run, OutputsModel)
     }
 
     # Run the model for the sub-basin
-    OutputsModel[[IM$id]] <- RunModel.InputsModel(
-      IM,
-      RunOptions = RunOptions[[IM$id]],
-      Param = Param[[IM$id]]
+    OutputsModel[[id]] <- RunModel.InputsModel(
+      x[[id]],
+      RunOptions = RunOptions[[id]],
+      Param = Param[[id]]
     )
-
   }
+  attr(OutputsModel, "Qm3s") <- OutputsModelQsim(x, OutputsModel, RunOptions[[1]]$IndPeriod_Run)
   return(OutputsModel)
 }
diff --git a/R/RunModel.Supervisor.R b/R/RunModel.Supervisor.R
index 285a3eb..22542c8 100644
--- a/R/RunModel.Supervisor.R
+++ b/R/RunModel.Supervisor.R
@@ -70,5 +70,6 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
   for(id in getSD_Ids(x$InputsModel)) {
     x$OutputsModel[[id]]$Qsim <- Qsim[[id]]
   }
+  attr(x$OutputsModel, "Qm3s") <- OutputsModelQsim(x$InputsModel, x$OutputsModel)
   return(x$OutputsModel)
 }
diff --git a/R/utils.R b/R/utils.R
index 49be122..52105ad 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -97,3 +97,55 @@ checkRunModelParameters <- function(InputsModel, RunOptions, Param) {
   if(!is.list(Param) || !all(names(InputsModel) %in% names(Param))) stop("Argument `Param` must be a list with names equal to nodes IDs")
   return()
 }
+
+
+#' Create a data.frame with simulated flows at each nodes of the [GRiwrm] object
+#'
+#' @details
+#' This function can only be called inside [RunModel.GRiwrmInputsmodel] or [RunModel.Supervisor]
+#' because it needs a `GRiwrmInputsModel` object internally modified by these functions
+#' (`Qupstream` updated with simulated flows).
+#'
+#' @param InputsModel a `GRiwrmInputsModel` object created by [CreateInputsModel.GRiwrm]
+#' @param OutputsModel a `GRiwrmOutputsModel` object created by [RunModel.GRiwrmInputsmodel] or [RunModel.Supervisor]
+#' @param IndPeriod_Run an [integer] vector (See [airGR::CreateRunOptions])
+#'
+#' @return a [data.frame] containing the simulated flows (in m3/time step) structured with the following columns:
+#' - 'DatesR' containing the timestamps of the time series
+#' - one column by node with the simulated flows
+#'
+OutputsModelQsim <- function(InputsModel, OutputsModel, IndPeriod_Run) {
+  griwrm <- attr(InputsModel, "GRiwrm")
+  # Get simulated flow for each node
+  # Flow for each node is available in InputsModel$Qupstream except for the downstream node
+  upperNodes <- griwrm$id[!is.na(griwrm$down)]
+  lQsim <- lapply(
+    upperNodes,
+    function(x, griwrm, IndPeriod_Run) {
+      node <- griwrm$down[griwrm$id == x]
+      InputsModel[[node]]$Qupstream[IndPeriod_Run, x]
+    },
+    griwrm = griwrm, IndPeriod_Run = IndPeriod_Run
+  )
+  names(lQsim) <- upperNodes
+  # Flow of the downstream node is only available in OutputsModel[[node]]$Qsim
+  downNode <- names(InputsModel)[length(InputsModel)]
+  lQsim[[downNode]] <- OutputsModel[[downNode]]$Qsim
+
+  # Conversion to m3/s
+  lQsim <- lapply(
+    names(lQsim),
+    function(x) {
+      i <- which(griwrm$id == x)
+      if(is.na(griwrm$area[i])) { # m3/time step => m3/s
+        return(lQsim[[x]] / attr(InputsModel, "TimeStep"))
+      } else { # mm/time step => m3/s
+        return(lQsim[[x]] * griwrm$area[i] * 1E3 / attr(InputsModel, "TimeStep"))
+      }
+    }
+  )
+  names(lQsim) <- c(upperNodes, downNode)
+  dfQsim <- cbind(data.frame(DatesR = as.POSIXct(InputsModel[[1]]$DatesR[IndPeriod_Run])),
+                  do.call(cbind,lQsim))
+  return(dfQsim)
+}
-- 
GitLab