feat: Overriding airGR CreateInputsModel, RunOptions and RunModel

Fix #4

NAMESPACE

 # Generated by roxygen2: do not edit by hand
 
+S3method(Calibration,default)
+S3method(Calibration,griwrm)
+S3method(CreateInputsModel,Griwrm)
+S3method(CreateInputsModel,default)
+S3method(CreateRunOptions,GriwrmInputsModel)
+S3method(CreateRunOptions,InputsModel)
+S3method(RunModel,GriwrmInputsModel)
+S3method(RunModel,InputsModel)
 S3method(merge,Gits)
+export(Calibration)
+export(CreateInputsModel)
+export(CreateRunOptions)
 export(Ginet)
 export(Girop)
 export(Gits)
-export(RunModelGriwrm)
-export(SetAirGrInputsAndOptions)
+export(RunModel)
 export(getNodeRanking)
 import(airGR)

R/RunModelGriwrm.R → R/CreateInputsModel.Griwrm.R

-#' Title
+#' Create InputsModel object for a GRIWRM network
 #'
 #' @param ginet
 #' @param girop
 #' @param gits
-#' @param IndPeriod_Run
-#' @param IndPeriod_WarmUp
 #'
 #' @return
-#' @import airGR
 #' @export
 #'
 #' @examples
-RunModelGriwrm <- function(ginet, girop, gits, IndPeriod_Run, IndPeriod_WarmUp = NULL, verbose = TRUE) {
+CreateInputsModel.Griwrm <- function(ginet, girop, gits, verbose = TRUE) {
 
-  OutputsModels <- list()
+  InputsModel <- CreateEmptyGriwrmInputsModel()
 
   for(id in getNodeRanking(ginet)) {
-    if(verbose) cat("*** Treating sub-basin", id, "... ***\n")
-
-    # Set InputsModel and RunOptions
-    lIO <- SetAirGrInputsAndOptions(
-      id, ginet, girop, gits, OutputsModels,
-      IndPeriod_Run, IndPeriod_WarmUp
-    )
-
-    # Prepare param for upstream sub-basin or basin with hydraulic routing
-    Param <- unlist(girop$params[girop$id == id])
+    if(verbose) cat("CreateInputsModel.griwrm: Treating sub-basin", id, "...\n")
+    InputsModel[[id]] <- CreateOneGriwrmInputsModel(id, ginet, girop, gits)
+  }
+  return(InputsModel)
+}
 
-    # Run the model for the sub-basin
-    OutputsModels[[id]] <- RunModel(
-      FUN_MOD = girop$model[girop$id == id],
-      InputsModel = lIO$InputsModel,
-      RunOptions = lIO$RunOptions,
-      Param = Param
-    )
 
-  }
-  OutputsModels
+#' Create an empty InputsModel object for GRIWRM nodes
+#'
+#' @return
+#'
+#' @examples
+CreateEmptyGriwrmInputsModel <- function() {
+  InputsModel <- list()
+  class(InputsModel) <- append(class(InputsModel), "GriwrmInputsModel")
+  return(InputsModel)
 }
 
-#' Title
+
+#' Create one InputsModel for a GRIWRM node
 #'
-#' @param id
 #' @param ginet
 #' @param girop
 #' @param gits
-#' @param IndPeriod_Run
-#' @param IndPeriod_WarmUp
+#' @param id
 #'
 #' @return
-#' @import airGR
-#' @export
 #'
 #' @examples
-SetAirGrInputsAndOptions <- function(id, ginet, girop, gits, OutputsModels, IndPeriod_Run, IndPeriod_WarmUp = NULL) {
+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)]
-  if(length(UpstrNodes) == 0) {
-    # Upstream sub-basin
-    QobsUpstr <- NULL
-    LengthHydro <- NULL
-    BasinAreas <- NULL
-  } else {
+  UpstreamNodes <- ginet$id[ginet$down == id & !is.na(ginet$down)]
+  QobsUpstr <- NULL
+  LengthHydro <- NULL
+  BasinAreas <- NULL
+
+  if(length(UpstreamNodes) > 0) {
     # Sub-basin with hydraulic routing
-    for(i in 1:length(UpstrNodes)) {
-      QobsUpstr1 <- matrix(
-        c(
-          rep(0, length(IndPeriod_WarmUp)),
-          OutputsModels[[UpstrNodes[i]]]$Qsim
-        ), ncol = 1
-      )
-      if(i == 1) {
+    for(idUpstrNode in UpstreamNodes) {
+      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)
+    LengthHydro <- matrix(ginet$length[girop$id %in% UpstreamNodes] , nrow = 1)
     BasinAreas <- matrix(
       c(
-        girop$area[girop$id %in% UpstrNodes],
-        girop$area[girop$id == id] - sum(girop$area[girop$id %in% UpstrNodes])
+        girop$area[girop$id %in% UpstreamNodes],
+        girop$area[girop$id == id] - sum(girop$area[girop$id %in% UpstreamNodes])
       ),
       nrow = 1
     )
   }
-  # Set model inputs
+
+  # Set model inputs with the airGR function
   InputsModel <- CreateInputsModel(
-    FUN_MOD = FUN_MOD,
+    FUN_MOD,
     DatesR = gits$date,
     Precip = gits[[id]]$Precip,
     PotEvap = gits[[id]]$PotEvap,
@@ -97,11 +82,14 @@ SetAirGrInputsAndOptions <- function(id, ginet, girop, gits, OutputsModels, IndP
     LengthHydro = LengthHydro,
     BasinAreas = BasinAreas
   )
-  # Set model options
-  RunOptions <- CreateRunOptions(
-    FUN_MOD = FUN_MOD,
-    InputsModel = InputsModel, IndPeriod_Run = IndPeriod_Run,
-    IniStates = NULL, IniResLevels = NULL, IndPeriod_WarmUp = IndPeriod_WarmUp
-  )
-  list(InputsModel = InputsModel, RunOptions = RunOptions)
+
+  # Add Identifiers of connected nodes in order to be able to update them with simulated flows
+  InputsModel$id <- id
+  if(length(UpstreamNodes) > 0) {
+    InputsModel$UpstreamNodes <- UpstreamNodes
+  }
+  # Add the model function
+  InputsModel$FUN_MOD <- FUN_MOD
+
+  return(InputsModel)
 }

R/CreateInputsModel.R

@@ -10,126 +10,3 @@
 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
-  )
-}

R/CreateInputsModel.default.R

+#' 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
+#' @import airGR
+#' @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)
+}

R/CreateRunOptions.GriwrmInputsModel.R

+#' Title
+#'
+#' @param InputsModel
+#' @param IndPeriod_WarmUp
+#' @param IndPeriod_Run
+#' @param IniStates
+#' @param IniResLevels
+#' @param Imax
+#' @param Outputs_Cal
+#' @param Outputs_Sim
+#' @param MeanAnSolidPrecip
+#' @param IsHyst
+#' @param warnings
+#' @param verbose
+#'
+#' @return
+#' @export
+#'
+#' @examples
+CreateRunOptions.GriwrmInputsModel <- function(InputsModel, IndPeriod_WarmUp = NULL, IndPeriod_Run,
+                                               IniStates = NULL, IniResLevels = NULL, Imax = NULL,
+                                               Outputs_Cal = NULL, Outputs_Sim = "all",
+                                               MeanAnSolidPrecip = NULL, IsHyst = FALSE,
+                                               warnings = TRUE, verbose = TRUE) {
+  RunOptions <- list()
+  class(RunOptions) <- append(class(RunOptions), "GriwrmRunOptions")
+  for(InputsModelBasin in InputsModel) {
+    RunOptions[[InputsModelBasin$id]] <- CreateRunOptions(
+      InputsModel = InputsModelBasin,
+      IndPeriod_WarmUp = IndPeriod_WarmUp,
+      IndPeriod_Run = IndPeriod_Run,
+      IniStates = IniStates,
+      IniResLevels = IniResLevels,
+      Imax = Imax,
+      Outputs_Cal = Outputs_Cal,
+      Outputs_Sim = Outputs_Sim,
+      MeanAnSolidPrecip = MeanAnSolidPrecip,
+      IsHyst = IsHyst,
+      warnings = warnings,
+      verbose = verbose
+    )
+  }
+  return(RunOptions)
+}

R/CreateRunOptions.InputsModel.R

+#' Title
+#'
+#' @param FUN_MOD
+#' @param InputsModel
+#' @param IndPeriod_WarmUp
+#' @param IndPeriod_Run
+#' @param IniStates
+#' @param IniResLevels
+#' @param Imax
+#' @param Outputs_Cal
+#' @param Outputs_Sim
+#' @param MeanAnSolidPrecip
+#' @param IsHyst
+#' @param warnings
+#' @param verbose
+#'
+#' @return
+#' @export
+#'
+#' @examples
+CreateRunOptions.InputsModel <- function(InputsModel,
+                                         IndPeriod_WarmUp = NULL, IndPeriod_Run,
+                                         IniStates = NULL, IniResLevels = NULL, Imax = NULL,
+                                         Outputs_Cal = NULL, Outputs_Sim = "all",
+                                         MeanAnSolidPrecip = NULL, IsHyst = FALSE,
+                                         warnings = TRUE, verbose = TRUE) {
+  airGR::CreateRunOptions(FUN_MOD = InputsModel$FUN_MOD,
+                          InputsModel = InputsModel,
+                          IndPeriod_WarmUp = IndPeriod_WarmUp,
+                          IndPeriod_Run = IndPeriod_Run,
+                          IniStates = IniStates,
+                          IniResLevels = IniResLevels,
+                          Imax = Imax,
+                          Outputs_Cal = Outputs_Cal,
+                          Outputs_Sim = Outputs_Sim,
+                          MeanAnSolidPrecip = MeanAnSolidPrecip,
+                          IsHyst = IsHyst,
+                          warnings = warnings,
+                          verbose = verbose)
+}

R/CreateRunOptions.R

+#' Title
+#'
+#' @param ...
+#' @param InputsModel
+#'
+#' @return
+#' @export
+#'
+#' @examples
+CreateRunOptions <- function(InputsModel, ...) {
+  UseMethod("CreateRunOptions", InputsModel)
+}

R/RunModel.GriwrmInputsModel.R

+#' Title
+#'
+#' @param ginet
+#' @param girop
+#' @param gits
+#' @param IndPeriod_Run
+#' @param IndPeriod_WarmUp
+#'
+#' @return
+#' @export
+#'
+#' @examples
+RunModel.GriwrmInputsModel <- function(InputsModel, RunOptions, girop, verbose = TRUE) {
+
+  OutputsModels <- list()
+
+  for(IM in InputsModel) {
+    if(verbose) cat("RunModel.GriwrmInputsModel: Treating sub-basin", IM$id, "...\n")
+
+    # Update InputsModel$QobsUpstr with simulated upstream flows
+    if(length(IM$UpstreamNodes) > 0) {
+      for(i in 1:length(IM$UpstreamNodes)) {
+        QobsUpstr1 <- matrix(
+          c(
+            rep(0, length(RunOptions[[IM$id]]$IndPeriod_WarmUp)),
+            OutputsModels[[IM$UpstreamNodes[i]]]$Qsim
+          ), ncol = 1
+        )
+        if(i == 1) {
+          IM$QobsUpstr <- QobsUpstr1
+        } else {
+          IM$QobsUpstr <- cbind(IM$QobsUpstr, QobsUpstr1)
+        }
+      }
+    }
+
+    # Run the model for the sub-basin
+    OutputsModels[[IM$id]] <- RunModel(
+      InputsModel = IM,
+      RunOptions = RunOptions[[IM$id]],
+      Param = unlist(girop$params[girop$id == IM$id])
+    )
+
+  }
+  return(OutputsModels)
+}

R/RunModel.InputsModel.R

+#' Wrapper for \code{\link[airGR]{RunModel}} which performs a single model run with the provided function over the selected period.
+#'
+#' @param InputsModel object of class \emph{InputsModel}, see \code{\link[airGR]{CreateInputsModel}} for details.
+#' @param RunOptions object of class \emph{RunOptions}, see \code{\link[airGR]{CreateRunOptions}} for details.
+#' @param Param numeric vector of model parameters.
+#' @param FUN_MOD hydrological model function (e.g. \code{\link[airGR]{RunModel_GR4J}}, \code{\link[airGR]{RunModel_CemaNeigeGR4J}}).
+#'
+#' @return
+#' @export
+#'
+#' @examples
+RunModel.InputsModel <- function(InputsModel, RunOptions, Param, FUN_MOD = NULL) {
+  if(is.null(FUN_MOD)) {
+    FUN_MOD <- InputsModel$FUN_MOD
+  }
+  airGR::RunModel(InputsModel, RunOptions, Param, FUN_MOD)
+}

R/RunModel.R

+#' RunModel function for both airGR and GriwrmInputsModel object
+#'
+#' @param InputsModel
+#' @param ...
+#'
+#' @return
+#' @export
+#'
+#' @examples
+RunModel <- function(InputsModel, ...) {
+  UseMethod("RunModel", InputsModel)
+}

vignettes/01_First_network.Rmd

@@ -80,10 +80,9 @@ load_ts <- function(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[nrow(ts)] <- 0 # No value at the end: time serie converge to zero
   }
   ts$Qnat <- zoo::na.approx(ts$Qnat)
   ts
@@ -120,7 +119,7 @@ The airGR CreateInputsModel function is extended in order to handle the ginet ob
 
 
 ```{r}
-  InputsModel <- CreateInputsModel(ginet, girop, gits)
+InputsModel <- CreateInputsModel(ginet, girop, gits)
 ```
 
 
@@ -129,6 +128,6 @@ The airGR CreateInputsModel function is extended in order to handle the ginet ob
 
 ```{r}
 dir.create("_cache", showWarnings = FALSE)
-save(ginet, girop, gits, file = "_cache/seine.RData")
+save(ginet, girop, gits, InputsModel, file = "_cache/seine.RData")
 ```
 

vignettes/02_First_run.Rmd

@@ -9,13 +9,25 @@ output: html_document
 knitr::opts_chunk$set(echo = TRUE)
 ```
 
-## Load parameters of GR4 run-off model
+## Load libraries
 
 ```{r}
 library(griwrm)
+```
+
+
+## Load parameters of GR4 run-off model
+
+### Loading network and time series data
+
+Run `vignette("01_First_network", package = "griwrm")` before this one in order to create the Rdata file loaded below:
+
+```{r}
 load("_cache/seine.RData")
 ```
 
+### Loading 
+
 Data comes from calibration of ClimAware project with naturalised flows.
 
 ```{r}
@@ -27,7 +39,7 @@ ClimAwareParams
 ```
 
 The lag $\tau_0$ and route $K_0$ parameters of TGR are expressed as time delay in hours corresponding to the delay time between the farest upstream inlet and the outlet of the sub-basin. 
-Almost all sub basin has only a lag parameter. The only exception is for La Marne à Noisiel (NOISI_17) that has a routing parameter which can be approximated to a single lag parameter equals to $\tau_0 + K_0$.
+Almost all sub basin has only a lag parameter. The only exception is for La Marne à Noisiel (NOISI_17) that has a routing parameter which can be approximated to a single lag parameter equals to $\tau_0 + K_0$.
 
 This lag parameter has to be converted in a speed in m/s used in the airGR lag model:
 
@@ -48,22 +60,52 @@ for(id in girop$id) {
   girop$params[girop$id == id] <- list(params[params$id == id, c("S", "IGF", "KR", "T", sLag)])
 }
 
+dplyr::select(params, id, S, IGF, KR, T, LAG)
 ```
 
-## Run the SD model for the whole basin
+## GriwmRunOptions object
+
+The `CreateRunOptions()` function allows to prepare the options required to the `RunModel()` function.
+
+The user must at least define the following arguments:
+
+* InputsModel: the associated input data
+* IndPeriod_Run: the period on which the model is run
+
 
 ```{r}
-# Time settings
 library(lubridate)
 IndPeriod_Run <- seq(
-  which(gits$date == (gits$date[1] + months(12))), # Set aside warm-up period
-  length(gits$date)) # Until the end of the time series
+  which(InputsModel[[1]]$DatesR == (InputsModel[[1]]$DatesR[1] + months(12))), # Set aside warm-up period
+  length(InputsModel[[1]]$DatesR) # Until the end of the time series
+) 
+```
+
+The warmup period could also be defined as is:
+
+```{r}
 IndPeriod_WarmUp = seq(1,IndPeriod_Run[1]-1)
+```
 
-OutputsModels <- RunModelGriwrm(
-  ginet, girop, gits, 
-  IndPeriod_Run = IndPeriod_Run, 
-  IndPeriod_WarmUp = IndPeriod_WarmUp
+
+
+```{r}
+RunOptions <- CreateRunOptions(
+  InputsModel = InputsModel, 
+  IndPeriod_WarmUp = IndPeriod_WarmUp, 
+  IndPeriod_Run = IndPeriod_Run
+)
+```
+
+
+
+## Run the SD model for the whole basin
+
+```{r}
+OutputsModels <- RunModel(
+  InputsModel = InputsModel,
+  RunOptions = RunOptions,
+  girop
 )
 ```