Resolve "Handle CemaNeige compatibility"

R/CreateInputsModel.GRiwrm.R

@@ -5,7 +5,18 @@
 #' @param Precip Matrix or data frame of numeric containing precipitation in mm. Column names correspond to node IDs
 #' @param PotEvap Matrix or data frame of numeric containing potential evaporation in mm. Column names correspond to node IDs
 #' @param Qobs Matrix or data frame of numeric containing potential observed flow in mm. Column names correspond to node IDs
-#' @param ... further arguments passed to [airGR::CreateInputsModel]
+#' @param PrecipScale (optional) named [vector] of [logical] indicating if the mean of the precipitation interpolated on the elevation layers must be kept or not, required to create CemaNeige module inputs, default `TRUE` (the mean of the precipitation is kept to the original value)
+#' @param TempMean (optional) [matrix] or [data.frame] of time series of mean air temperature \[°C\], required to create the CemaNeige module inputs
+#' @param TempMin (optional) [matrix] or [data.frame] of time series of minimum air temperature \[°C\], possibly used to create the CemaNeige module inputs
+#' @param TempMax (optional) [matrix] or [data.frame] of time series of maximum air temperature \[°C\], possibly used to create the CemaNeige module inputs
+#' @param ZInputs  (optional) named [vector] of [numeric] giving the mean elevation of the Precip and Temp series (before extrapolation) \[m\], possibly used to create the CemaNeige module input
+#' @param HypsoData	(optional) [matrix] or [data.frame] containing 101 [numeric] rows: min, q01 to q99 and max of catchment elevation distribution \[m\], if not defined a single elevation is used for CemaNeige
+#' @param NLayers (optional) named [vector] of [numeric] integer giving the number of elevation layers requested [-], required to create CemaNeige module inputs, default=5
+#' @param ... used for compatibility with S3 methods
+#'
+#' @details Meteorological data are needed for the nodes of the network that represent a catchment simulated by a rainfall-runoff model. Instead of [airGR::CreateInputsModel] that has [numeric] [vector] as time series inputs, this function uses [matrix] or [data.frame] with the id of the sub-catchment as column names. For single values (`ZInputs` or `NLayers`), the function requires named [vector] with the id of the sub-catchment as name item. If an argument is optional, only the column or the named item has to be provided.
+#'
+#' See [airGR::CreateInputsModel] documentation for details concerning each input.
 #'
 #' @return GRiwrmInputsModel object equivalent to **airGR** InputsModel object for a semi-distributed model (See [airGR::CreateInputsModel])
 #' @export
@@ -53,16 +64,60 @@
 #'                             Qobs = Qobs)
 #' str(InputsModels)
 #'
-CreateInputsModel.GRiwrm <- function(x, DatesR, Precip, PotEvap, Qobs, ...) {
+CreateInputsModel.GRiwrm <- function(x, DatesR,
+                                     Precip,
+                                     PotEvap = NULL,
+                                     Qobs,
+                                     PrecipScale = TRUE,
+                                     TempMean = NULL, TempMin = NULL,
+                                     TempMax = NULL, ZInputs = NULL,
+                                     HypsoData = NULL, NLayers = 5, ...) {
+
+  # Check and format inputs
+  varNames <- c("Precip", "PotEvap", "TempMean",
+                "TempMin", "TempMax", "ZInputs", "HypsoData", "NLayers")
+  names(varNames) <- varNames
+  lapply(varNames, function(varName) {
+    v <- get(varName)
+    if(!is.null(v)) {
+      if(is.matrix(v) || is.data.frame(v)) {
+        if(is.null(colnames(v))) {
+          stop(sprintf(
+            "'%s' must have column names",
+            varName
+          ))
+        } else if(!all(colnames(v) %in% x$id)) {
+          stop(sprintf(
+            "'%s' column names must be included in 'id's of the GRiwrm object",
+            varName
+          ))
+        }
+      } else if (!varName %in% c("ZInputs", "NLayers")) {
+        stop(sprintf("'%s' must be a matrix or a data.frame", varName))
+      }
+    }
+  })
 
   InputsModel <- CreateEmptyGRiwrmInputsModel(x)
   Qobs[is.na(Qobs)] <- -99 # airGR::CreateInputsModel doesn't accept NA values
 
   for(id in getNodeRanking(x)) {
     message("CreateInputsModel.GRiwrm: Treating sub-basin ", id, "...")
-    InputsModel[[id]] <- CreateOneGRiwrmInputsModel(
-      id, x, DatesR,Precip[,id], PotEvap[,id], Qobs, ...
-    )
+    InputsModel[[id]] <-
+      CreateOneGRiwrmInputsModel(id = id,
+                                 griwrm = x,
+                                 DatesR = DatesR,
+                                 Precip = getInputBV(Precip, id),
+                                 PrecipScale,
+                                 PotEvap = getInputBV(PotEvap, id),
+                                 TempMean = getInputBV(TempMean, id),
+                                 TempMin = getInputBV(TempMin, id),
+                                 TempMax = getInputBV(TempMax, id),
+                                 ZInputs = getInputBV(ZInputs, id),
+                                 HypsoData = getInputBV(HypsoData, id),
+                                 NLayers = getInputBV(NLayers, id, 5),
+                                 Qobs = Qobs
+                                 )
   }
   attr(InputsModel, "TimeStep") <- getModelTimeStep(InputsModel)
   return(InputsModel)
@@ -74,6 +129,7 @@ CreateInputsModel.GRiwrm <- function(x, DatesR, Precip, PotEvap, Qobs, ...) {
 #' @param griwrm a `GRiwrm` object (See [CreateGRiwrm])
 #'
 #' @return \emph{GRiwrmInputsModel} empty object
+#' @noRd
 CreateEmptyGRiwrmInputsModel <- function(griwrm) {
   InputsModel <- list()
   class(InputsModel) <- c("GRiwrmInputsModel", class(InputsModel))
@@ -90,9 +146,10 @@ CreateEmptyGRiwrmInputsModel <- function(griwrm) {
 #' @param Precip time series of potential evapotranspiration (catchment average) (mm/time step)
 #' @param PotEvap time series of potential evapotranspiration (catchment average) (mm/time step)
 #' @param Qobs Matrix or data frame of numeric containing observed flow (mm/time step). Column names correspond to node IDs
-##'
+#'
 #' @return \emph{InputsModel} object for one.
-CreateOneGRiwrmInputsModel <- function(id, griwrm, DatesR, Precip, PotEvap, Qobs) {
+#' @noRd
+CreateOneGRiwrmInputsModel <- function(id, griwrm, ..., Qobs) {
   node <- griwrm[griwrm$id == id,]
   FUN_MOD <- griwrm$model[griwrm$id == id]
 
@@ -111,15 +168,19 @@ CreateOneGRiwrmInputsModel <- function(id, griwrm, DatesR, Precip, PotEvap, Qobs
         griwrm$area[griwrm$id %in% UpstreamNodes],
         node$area - sum(griwrm$area[griwrm$id %in% UpstreamNodes], na.rm = TRUE)
     )
+    if (BasinAreas[length(BasinAreas)] < 0) {
+      stop(sprintf(
+        "Area of the catchment %s must be greater than the sum of the areas of its upstream catchments",
+        id
+      ))
+    }
     names(BasinAreas) <- c(UpstreamNodes, id)
   }
 
   # Set model inputs with the **airGR** function
   InputsModel <- CreateInputsModel(
     FUN_MOD,
-    DatesR = DatesR,
-    Precip = Precip,
-    PotEvap = PotEvap,
+    ...,
     Qupstream = Qupstream,
     LengthHydro = LengthHydro,
     BasinAreas = BasinAreas
@@ -149,7 +210,7 @@ CreateOneGRiwrmInputsModel <- function(id, griwrm, DatesR, Precip, PotEvap, Qobs
 #' @param InputsModel a `GRiwrmInputsModel`
 #'
 #' @return A [numeric] representing the time step in seconds
-#'
+#' @noRd
 getModelTimeStep <- function(InputsModel) {
   TS <- sapply(InputsModel, function(x) {
     if (inherits(x, "hourly")) {
@@ -165,3 +226,32 @@ getModelTimeStep <- function(InputsModel) {
   }
   return(unique(TS))
 }
+
+#' Select the node input for input arguments of [airGR::CreateInputsModel]
+#'
+#' @param x [matrix] [data.frame] or named [vector] the input argument
+#' @param id [character] the id of the node
+#' @param unset default value if the id is not found in `x`
+#'
+#' @return the selected column or value in respect to `id`
+#' @noRd
+getInputBV <- function(x, id, unset = NULL) {
+  if(is.null(x)) {
+    return(unset)
+  }
+  if (is.matrix(x) || is.data.frame(x)) {
+    if (!id %in% colnames(x)) {
+      return(unset)
+    }
+  } else {
+    # vector (for ZInputs and NLayers)
+    if (length(x) == 1 && is.null(names(x))) {
+      return(x)
+    } else if(!id %in% names(x)) {
+      return(unset)
+    } else {
+      return(x[id])
+    }
+  }
+  return(x[, id])
+}

man/CreateEmptyGRiwrmInputsModel.Rd

-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/CreateInputsModel.GRiwrm.R
-\name{CreateEmptyGRiwrmInputsModel}
-\alias{CreateEmptyGRiwrmInputsModel}
-\title{Create an empty InputsModel object for \strong{airGRiwrm} nodes}
-\usage{
-CreateEmptyGRiwrmInputsModel(griwrm)
-}
-\arguments{
-\item{griwrm}{a \code{GRiwrm} object (See \link{CreateGRiwrm})}
-}
-\value{
-\emph{GRiwrmInputsModel} empty object
-}
-\description{
-Create an empty InputsModel object for \strong{airGRiwrm} nodes
-}

man/CreateInputsModel.GRiwrm.Rd

@@ -4,7 +4,21 @@
 \alias{CreateInputsModel.GRiwrm}
 \title{Create InputsModel object for a \strong{airGRiwrm} network}
 \usage{
-\method{CreateInputsModel}{GRiwrm}(x, DatesR, Precip, PotEvap, Qobs, ...)
+\method{CreateInputsModel}{GRiwrm}(
+  x,
+  DatesR,
+  Precip,
+  PotEvap = NULL,
+  Qobs,
+  PrecipScale = TRUE,
+  TempMean = NULL,
+  TempMin = NULL,
+  TempMax = NULL,
+  ZInputs = NULL,
+  HypsoData = NULL,
+  NLayers = 5,
+  ...
+)
 }
 \arguments{
 \item{x}{GRiwrm object describing the diagram of the semi-distributed model (See \link{CreateGRiwrm})}
@@ -17,7 +31,21 @@
 
 \item{Qobs}{Matrix or data frame of numeric containing potential observed flow in mm. Column names correspond to node IDs}
 
-\item{...}{further arguments passed to \link[airGR:CreateInputsModel]{airGR::CreateInputsModel}}
+\item{PrecipScale}{(optional) named \link{vector} of \link{logical} indicating if the mean of the precipitation interpolated on the elevation layers must be kept or not, required to create CemaNeige module inputs, default \code{TRUE} (the mean of the precipitation is kept to the original value)}
+
+\item{TempMean}{(optional) \link{matrix} or \link{data.frame} of time series of mean air temperature [°C], required to create the CemaNeige module inputs}
+
+\item{TempMin}{(optional) \link{matrix} or \link{data.frame} of time series of minimum air temperature [°C], possibly used to create the CemaNeige module inputs}
+
+\item{TempMax}{(optional) \link{matrix} or \link{data.frame} of time series of maximum air temperature [°C], possibly used to create the CemaNeige module inputs}
+
+\item{ZInputs}{(optional) named \link{vector} of \link{numeric} giving the mean elevation of the Precip and Temp series (before extrapolation) [m], possibly used to create the CemaNeige module input}
+
+\item{HypsoData}{(optional) \link{matrix} or \link{data.frame} containing 101 \link{numeric} rows: min, q01 to q99 and max of catchment elevation distribution [m], if not defined a single elevation is used for CemaNeige}
+
+\item{NLayers}{(optional) named \link{vector} of \link{numeric} integer giving the number of elevation layers requested \link{-}, required to create CemaNeige module inputs, default=5}
+
+\item{...}{used for compatibility with S3 methods}
 }
 \value{
 GRiwrmInputsModel object equivalent to \strong{airGR} InputsModel object for a semi-distributed model (See \link[airGR:CreateInputsModel]{airGR::CreateInputsModel})
@@ -25,6 +53,11 @@ GRiwrmInputsModel object equivalent to \strong{airGR} InputsModel object for a s
 \description{
 Create InputsModel object for a \strong{airGRiwrm} network
 }
+\details{
+Meteorological data are needed for the nodes of the network that represent a catchment simulated by a rainfall-runoff model. Instead of \link[airGR:CreateInputsModel]{airGR::CreateInputsModel} that has \link{numeric} \link{vector} as time series inputs, this function uses \link{matrix} or \link{data.frame} with the id of the sub-catchment as column names. For single values (\code{ZInputs} or \code{NLayers}), the function requires named \link{vector} with the id of the sub-catchment as name item. If an argument is optional, only the column or the named item has to be provided.
+
+See \link[airGR:CreateInputsModel]{airGR::CreateInputsModel} documentation for details concerning each input.
+}
 \examples{
 #################################################################
 # Run the `airGRRunModel_Lag` example in the GRiwrm fashion way #

man/CreateOneGRiwrmInputsModel.Rd

-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/CreateInputsModel.GRiwrm.R
-\name{CreateOneGRiwrmInputsModel}
-\alias{CreateOneGRiwrmInputsModel}
-\title{Create one InputsModel for a \strong{airGRiwrm} node}
-\usage{
-CreateOneGRiwrmInputsModel(id, griwrm, DatesR, Precip, PotEvap, Qobs)
-}
-\arguments{
-\item{id}{string of the node identifier}
-
-\item{griwrm}{See \link{CreateGRiwrm})}
-
-\item{DatesR}{vector of dates required to create the GR model and CemaNeige module inputs}
-
-\item{Precip}{time series of potential evapotranspiration (catchment average) (mm/time step)}
-
-\item{PotEvap}{time series of potential evapotranspiration (catchment average) (mm/time step)}
-
-\item{Qobs}{Matrix or data frame of numeric containing observed flow (mm/time step). Column names correspond to node IDs}
-}
-\value{
-\emph{InputsModel} object for one.
-}
-\description{
-Create one InputsModel for a \strong{airGRiwrm} node
-}

man/getModelTimeStep.Rd

-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/CreateInputsModel.GRiwrm.R
-\name{getModelTimeStep}
-\alias{getModelTimeStep}
-\title{Check time steps of the model of all the nodes and return the time step in seconds}
-\usage{
-getModelTimeStep(InputsModel)
-}
-\arguments{
-\item{InputsModel}{a \code{GRiwrmInputsModel}}
-}
-\value{
-A \link{numeric} representing the time step in seconds
-}
-\description{
-This function is called inside \link{CreateInputsModel.GRiwrm} for defining the time step of the big model.
-}

tests/testthat/helper_cemaneige.R

+#' Set up a data sample for Cemaneige tests
+#'
+#' A basin of 3 nodes with 2 identical upstream basins and a nul area downstream basin with a lag of 0 time step
+#'
+#' @return [list] with:
+#' - `BasinInfo`
+#' - `griwrm`
+#' - `Precip`
+#' - `TempMean`
+#' - `PotEvap`
+#' - `Qobs`
+#' - `HypsoData`
+#' @noRd
+#'
+#' @examples
+setUpCemaNeigeData <- function() {
+
+  library(airGR)
+  data(L0123001)
+  detach("package:airGR")
+
+  # Formatting observations for the hydrological models
+  # Each input data should be a matrix or a data.frame with the good id in the name of the column
+  ids <- c("Up1", "Up2", "Down")
+  l <- lapply(c("P", "T", "E", "Qmm"), function(x) {
+    m <- matrix(data = rep(BasinObs[[x]], 3), ncol = 3)
+    colnames(m) <- ids
+    return(m)
+  })
+  l$DatesR <- BasinObs$DatesR
+  names(l) <- c("Precip", "TempMean", "PotEvap", "Qobs")
+  l$HypsoData <- matrix(data = rep(BasinInfo$HypsoData, 3), ncol = 3)
+  colnames(l$HypsoData) <- ids
+  l$ZInputs <- rep(BasinInfo$HypsoData[51], 3)
+  names(l$ZInputs) <- ids
+
+  db <- data.frame(id = ids,
+                   length = c(0, 0, NA),
+                   down = c("Down", "Down", NA),
+                   area = c(rep(BasinInfo$BasinArea, 2), BasinInfo$BasinArea * 2),
+                   model = rep("RunModel_CemaNeigeGR4J", 3),
+                   stringsAsFactors = FALSE)
+
+  # Create GRiwrm object from the data.frame
+  l$griwrm <- CreateGRiwrm(db)
+
+  l$DatesR <- BasinObs$DatesR
+
+  return(l)
+}

tests/testthat/test-CreateInputsModel.R

+context("CreateInputsModel")
+
+l <- setUpCemaNeigeData()
+
+test_that("CemaNeige data should be in InputsModel", {
+  InputsModels <- CreateInputsModel(l$griwrm,
+                                    DatesR = l$DatesR,
+                                    Precip = l$Precip,
+                                    PotEvap = l$PotEvap,
+                                    TempMean = l$TempMean,
+                                    ZInputs = l$ZInputs,
+                                    HypsoData = l$HypsoData,
+                                    Qobs = l$Qobs)
+  l$DatesR <- as.data.frame(l$DatesR)
+  lapply(InputsModels, function(IM) {
+    lapply(c("DatesR", "Precip", "PotEvap"), function(varName) {
+      expect_equal(IM[[varName]], l[[varName]][,1])
+    })
+    expect_named(IM$LayerPrecip, paste0("L", seq(1,5)))
+    expect_named(IM$LayerTempMean, paste0("L", seq(1,5)))
+    expect_named(IM$LayerFracSolidPrecip, paste0("L", seq(1,5)))
+  })
+})
+
+test_that("downstream sub-catchment area should be positive", {
+  l$griwrm$area[3] <- 360
+  expect_error(CreateInputsModel(l$griwrm,
+                                 DatesR = l$DatesR,
+                                 Precip = l$Precip,
+                                 PotEvap = l$PotEvap,
+                                 TempMean = l$TempMean,
+                                 ZInputs = l$ZInputs,
+                                 HypsoData = l$HypsoData,
+                                 Qobs = l$Qobs))
+})
+
+test_that("handles mix of with and without CemaNeige nodes", {
+  l$griwrm[l$griwrm$id == "Down", "model"] <- "RunModel_GR4J"
+  l$TempMean <- l$TempMean[,1:2]
+  l$ZInputs <- l$ZInputs[1:2]
+  l$TempMean <- l$TempMean[,1:2]
+  l$HypsoData <- l$HypsoData[,1:2]
+  InputsModels <- CreateInputsModel(l$griwrm,
+                                    DatesR = l$DatesR,
+                                    Precip = l$Precip,
+                                    PotEvap = l$PotEvap,
+                                    TempMean = l$TempMean,
+                                    ZInputs = l$ZInputs,
+                                    HypsoData = l$HypsoData,
+                                    Qobs = l$Qobs)
+  expect_false(inherits(InputsModels$Down, "CemaNeige"))
+  expect_null(InputsModels$Down$LayerPrecip)
+})
+
+test_that("throws error on wrong column name", {
+  colnames(l$Precip)[1] <- "Up0"
+  expect_error(CreateInputsModel(l$griwrm,
+                                 DatesR = l$DatesR,
+                                 Precip = l$Precip,
+                                 PotEvap = l$PotEvap,
+                                 TempMean = l$TempMean,
+                                 ZInputs = l$ZInputs,
+                                 HypsoData = l$HypsoData,
+                                 Qobs = l$Qobs))
+  colnames(l$Precip) <- NULL
+  expect_error(CreateInputsModel(l$griwrm,
+                                 DatesR = l$DatesR,
+                                 Precip = l$Precip,
+                                 PotEvap = l$PotEvap,
+                                 TempMean = l$TempMean,
+                                 ZInputs = l$ZInputs,
+                                 HypsoData = l$HypsoData,
+                                 Qobs = l$Qobs))
+})
+
+test_that("throws error when missing CemaNeige data", {
+  expect_error(CreateInputsModel(l$griwrm,
+                                 DatesR = l$DatesR,
+                                 Precip = l$Precip,
+                                 PotEvap = l$PotEvap,
+                                 Qobs = l$Qobs))
+})

tests/testthat/test-RunModel.R

@@ -12,8 +12,9 @@ griwrm <- CreateGRiwrm(nodes, list(id = "gauge_id", down = "downstream_id", leng
 
 # InputsModel
 DatesR <- Severn$BasinsObs[[1]]$DatesR
-PrecipTot <- cbind(sapply(Severn$BasinsObs, function(x) {x$precipitation}))
-PotEvapTot <- cbind(sapply(Severn$BasinsObs, function(x) {x$peti}))
+BasinsObs <- Severn$BasinsObs[griwrm$id]
+PrecipTot <- cbind(sapply(BasinsObs, function(x) {x$precipitation}))
+PotEvapTot <- cbind(sapply(BasinsObs, function(x) {x$peti}))
 Precip <- ConvertMeteoSD(griwrm, PrecipTot)
 PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
 Qobs <- cbind(sapply(Severn$BasinsObs, function(x) {x$discharge_spec}))
@@ -97,3 +98,41 @@ test_that("RunModelSupervisor with multi time steps controller, two regulations
   )
   expect_equal(OM_Supervisor[["54057"]]$Qsim, OM_GriwrmInputs[["54057"]]$Qsim)
 })
+
+test_that("RunModel.GRiwrm handles CemaNeige", {
+  l <- setUpCemaNeigeData()
+  l$griwrm[l$griwrm$id == "Down", "model"] <- "RunModel_GR4J"
+  l$TempMean <- l$TempMean[,1:2]
+  l$ZInputs <- l$ZInputs[1:2]
+  l$TempMean <- l$TempMean[,1:2]
+  l$HypsoData <- l$HypsoData[,1:2]
+  InputsModels <- CreateInputsModel(l$griwrm,
+                                    DatesR = l$DatesR,
+                                    Precip = l$Precip,
+                                    PotEvap = l$PotEvap,
+                                    TempMean = l$TempMean,
+                                    ZInputs = l$ZInputs,
+                                    HypsoData = l$HypsoData,
+                                    Qobs = l$Qobs)
+  ## run period selection
+  Ind_Run <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1990-01-01"),
+                 which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1999-12-31"))
+  ## preparation of the RunOptions object
+  RunOptions <- suppressWarnings(CreateRunOptions(InputsModels,
+                                                  IndPeriod_Run = Ind_Run))
+  ids <- l$griwrm$id
+  names(ids) <- ids
+  Params <- lapply(ids, function(x) {
+    c(X1 = 408.774, X2 = 2.646, X3 = 131.264, X4 = 1.174,
+      CNX1 = 0.962, CNX2 = 2.249)
+  })
+  Params$Down <- c(1, Params$Down[1:4])
+  OutputsModel <- RunModel(
+    InputsModels,
+    RunOptions = RunOptions,
+    Param = Params
+  )
+  expect_named(OutputsModel, l$griwrm$id)
+  Qm3s <- attr(OutputsModel, "Qm3s")
+  expect_equal(Qm3s[,4], rowSums(Qm3s[,2:3]))
+})