From 8c715e0fd8c10d7cb466ad462acf906a875b3a6e Mon Sep 17 00:00:00 2001
From: Dorchies David <david.dorchies@irstea.fr>
Date: Mon, 25 May 2020 16:27:07 +0200
Subject: [PATCH] feat: CreateInputsModel for GRIWRM

- Works in vignette 01

Refs #4
---
 R/CreateInputsModel.R          | 135 +++++++++++++++++++++++++++++++++
 R/ginet.R                      |   2 +-
 R/gits.R                       |   8 ++
 vignettes/01_First_network.Rmd |  59 ++++++++++++--
 4 files changed, 196 insertions(+), 8 deletions(-)
 create mode 100644 R/CreateInputsModel.R

diff --git a/R/CreateInputsModel.R b/R/CreateInputsModel.R
new file mode 100644
index 0000000..72eb3fa
--- /dev/null
+++ b/R/CreateInputsModel.R
@@ -0,0 +1,135 @@
+#' Create InputsModel object for either airGR or GRIWRM
+#'
+#' @param x
+#' @param ...
+#'
+#' @return
+#' @export
+#'
+#' @examples
+CreateInputsModel <- function(x, ...) {
+  UseMethod("CreateInputsModel", x)
+}
+
+#' Wrapper for the airGR::CreateInputsModel function
+#'
+#' @param FUN_MOD
+#' @param DatesR
+#' @param Precip
+#' @param PrecipScale
+#' @param PotEvap
+#' @param TempMean
+#' @param TempMin
+#' @param TempMax
+#' @param ZInputs
+#' @param HypsoData
+#' @param NLayers
+#' @param QobsUpstr
+#' @param LengthHydro
+#' @param BasinAreas
+#' @param verbose
+#'
+#' @return
+#' @export
+#'
+#' @examples
+CreateInputsModel.default <- function(FUN_MOD,
+                              DatesR,
+                              Precip, PrecipScale = TRUE,
+                              PotEvap = NULL,
+                              TempMean = NULL, TempMin = NULL, TempMax = NULL,
+                              ZInputs = NULL, HypsoData = NULL, NLayers = 5,
+                              QobsUpstr = NULL, LengthHydro = NULL, BasinAreas = NULL,
+                              verbose = TRUE) {
+  airGR::CreateInputsModel(FUN_MOD, DatesR, Precip, PrecipScale, PotEvap,
+                           TempMean, TempMin, TempMax, ZInputs, HypsoData, NLayers,
+                           QobsUpstr, LengthHydro, BasinAreas, verbose)
+}
+
+
+#' Create InputsModel object for a GRIWRM network
+#'
+#' @param ginet
+#' @param girop
+#' @param gits
+#'
+#' @return
+#' @export
+#'
+#' @examples
+CreateInputsModel.Griwrm <- function(ginet, girop, gits, verbose = TRUE) {
+
+  InputsModel <- CreateEmptyGriwrmInputsModel()
+
+  for(id in getNodeRanking(ginet)) {
+    if(verbose) cat("CreateInputsModel.griwrm: Treating sub-basin", id, "...\n")
+    InputsModel[[id]] <- CreateOneGriwrmInputsModel(id, ginet, girop, gits)
+  }
+}
+
+
+#' Create an empty InputsModel object for GRIWRM
+#'
+#' @return
+#' @export
+#'
+#' @examples
+CreateEmptyGriwrmInputsModel <- function() {
+  InputsModel <- list()
+  class(InputsModel) <- append(class(InputsModel), "GriwrmInputsModel")
+  return(InputsModel)
+}
+
+
+#' Create one InputsModel for a GRIWRM node
+#'
+#' @param ginet
+#' @param girop
+#' @param gits
+#' @param id
+#'
+#' @return
+#' @export
+#'
+#' @examples
+CreateOneGriwrmInputsModel <- function(id, ginet, girop, gits) {
+  node <- ginet[ginet$id == id,]
+  FUN_MOD <- girop$model[girop$id == id]
+
+  # Set hydraulic parameters
+  UpstrNodes <- ginet$id[ginet$down == id & !is.na(ginet$down)]
+  QobsUpstr <- NULL
+  LengthHydro <- NULL
+  BasinAreas <- NULL
+
+  if(length(UpstrNodes) > 0) {
+    # Sub-basin with hydraulic routing
+    for(idUpstrNode in UpstrNodes) {
+      QobsUpstr1 <- matrix(gits[[idUpstrNode]]$Qobs, ncol = 1)
+      if(is.null(QobsUpstr)) {
+        QobsUpstr <- QobsUpstr1
+      } else {
+        QobsUpstr <- cbind(QobsUpstr, QobsUpstr1)
+      }
+    }
+    LengthHydro <- matrix(ginet$length[girop$id %in% UpstrNodes] , nrow = 1)
+    BasinAreas <- matrix(
+      c(
+        girop$area[girop$id %in% UpstrNodes],
+        girop$area[girop$id == id] - sum(girop$area[girop$id %in% UpstrNodes])
+      ),
+      nrow = 1
+    )
+  }
+
+  # Set model inputs
+  CreateInputsModel(
+    FUN_MOD,
+    DatesR = gits$date,
+    Precip = gits[[id]]$Precip,
+    PotEvap = gits[[id]]$PotEvap,
+    QobsUpstr = QobsUpstr,
+    LengthHydro = LengthHydro,
+    BasinAreas = BasinAreas
+  )
+}
diff --git a/R/ginet.R b/R/ginet.R
index 6e8d5b9..7a95be0 100644
--- a/R/ginet.R
+++ b/R/ginet.R
@@ -16,7 +16,7 @@ Ginet <- function(db, cols = list(id = "id", down = "down", length = "length", r
   if(!keep_all) {
     db <- dplyr::select(db, names(cols))
   }
-  class(db) <- append(class(db), "Ginet")
+  class(db) <- append(class(db), c("Griwrm", "Ginet"))
   db
 }
 
diff --git a/R/gits.R b/R/gits.R
index eaf1fd9..ad0f31a 100644
--- a/R/gits.R
+++ b/R/gits.R
@@ -13,6 +13,14 @@ Gits <- function(id, ts,
 
   cols <- as.list(cols)
   ts <- dplyr::rename(ts, unlist(cols))
+
+  if(any(is.na(ts$Qobs))) {
+    stop("Qobs should not contain any NA")
+  }
+  if(any(ts$Qobs < 0)) {
+    stop("Qobs should be strictly positive")
+  }
+
   gitsOut <- list(date = ts$date)
   cols$date <- NULL
   gitsOut[[id]] <- dplyr::select(ts, names(cols))
diff --git a/vignettes/01_First_network.Rmd b/vignettes/01_First_network.Rmd
index 0db5a50..1c935c7 100644
--- a/vignettes/01_First_network.Rmd
+++ b/vignettes/01_First_network.Rmd
@@ -26,30 +26,47 @@ seine_nodes <- readr::read_delim(
 seine_nodes
 ```
 
-Create the ginet object
+Create the ginet object which lists the nodes and describes the network diagram. It's a dataframe of class `Ginet` and `Griwrm` with specific column names:
+
+- `id`: the identifier of the node in the network
+- `down`: the identifier of the next node downstream
+- `length`: hydraulic distance to the next downstream node
+- `runoff`: does the node is a rainfall run-off model?
+
+`Ginet` function helps to rename the columns of the dataframe and assign the variable classes.
 
 ```{r}
 # Specify that all nodes are of run-off type
 seine_nodes$runoff <- TRUE
 # Convert distance in km as it the unit used by airGR
-seine_nodes$distance_aval <- seine_nodes$distance_aval / 1000
-# Generate the ginet object
-ginet <- Ginet(seine_nodes, list(id = "id_sgl", down = "id_aval", length = "distance_aval"))
+seine_nodes$length <- seine_nodes$distance_aval / 1000
+# Generate the ginet object 
+ginet <- Ginet(seine_nodes, list(id = "id_sgl", down = "id_aval"))
+ginet
 ```
 
-Create the girop object
+Each line of the `Ginet` object having the `runoff` columns switched to `TRUE` should have a corresponding line in the `Girop` object which contains the parameters of the rainfall run-off models. 
+
+The `Girop` object is a dataframe of class `Girop` with specific columns:
+
+- `id`: the identifier of the node in the network
+- `area`: the total area of the basin (including upstream sub-basins) at the location of the node (km<sup>2</sup>)
+- `model`: the name of the rainfall run-off model used (e.g. "RunModel_GR4J")
+- `params`: a list containing the calibration parameters of the model
 
 ```{r}
 # Specify which run-off model to use
 seine_nodes$model = "RunModel_GR4J"
 # Generate girop object
 girop <- Girop(seine_nodes, list(id = "id_sgl", area = "area"))
+girop
 ```
 
 
-## Load data
+## Observation time series 
+
+Loading hydrometeorological data on the Seine river basin from the ClimAware project: 
 
-Hydrometeorological data on the Seine river basin from the ClimAware project: 
 ```{r, warning=FALSE, message=FALSE}
 urls <- 
   file.path(
@@ -62,6 +79,13 @@ load_ts <- function(x) {
   ts <- readr::read_delim(file = x, 
              delim = ";", skip = 16, trim_ws = TRUE)
   ts$Date <- as.POSIXlt(lubridate::ymd(ts$Date))
+  # Interpolation of data gap in the discharge time serie
+  browser
+  ts[ts$Qnat < 0, "Qnat"] <- NA
+  if(is.na(ts$Qnat[nrow(ts)])) {
+    ts$Qnat[nrow(ts)] <- 0 # End of time series converge to zero
+  }
+  ts$Qnat <- zoo::na.approx(ts$Qnat)
   ts
 }
 
@@ -69,9 +93,17 @@ l <- lapply(urls, load_ts)
 
 ```
 
+`Gits` object is a list containing a item named `date` with a timestamp vector of the time series and items named by the identifier of each node. These items contain a dataframe with the observations.
+
+The Gits function creates a `Gits` object 
+
 ```{r}
 gits <- Gits(ginet$id[1], l[[ginet$id[1]]], cols = list(date = "Date", Precip = "Ptot", PotEvap = "ETP", Qobs = "Qnat"))
+```
 
+Copy the observations for each node the ginet network:
+
+```{r}
 for(id in ginet$id) {
   l[[id]]$Qnat
   l[[id]]$Qnat <- l[[id]]$Qnat * 86.4 / girop$area[girop$id == id]
@@ -80,6 +112,19 @@ for(id in ginet$id) {
 }
 ```
 
+## Generate the GRIWRM InputsModel object
+
+The GRIWRM InputsModel object is a list of airGR InputsModel. The identifier of the sub-basin is used as key in the list which is ordered from upstream to downstream.
+
+The airGR CreateInputsModel function is extended in order to handle the ginet object which describe the basin diagram:
+
+
+```{r}
+  InputsModel <- CreateInputsModel(ginet, girop, gits)
+```
+
+
+
 ## Save data for next vignettes
 
 ```{r}
-- 
GitLab