v1.2.11.0 NEW: add an IsHyst argument in CreatCalibOptions to use hysteresis #5252

DESCRIPTION

 Package: airGR
 Type: Package
 Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
-Version: 1.2.10.0
+Version: 1.2.11.0
 Date: 2019-03-22
 Authors@R: c(
   person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),

NEWS.rmd

@@ -13,7 +13,7 @@ output:
 
 
 
-### 1.2.10.0 Release Notes (2019-03-22) 
+### 1.2.11.0 Release Notes (2019-03-22) 
 
 
 
@@ -50,6 +50,9 @@ output:
 - <code>CreateRunOptions()</code> now presents a <code>warnings</code> argument to replace the verbose action (the <code>verbose</code> argument is kept to print messages).
 
 
+- <code>CreateCalibOptions()</code> now presents a <code>IsHyst</code> argument to give the possibility to use the hysteresis with CemaNeige.
+
+
 
 #### Major user-visible changes
 

R/CreateCalibOptions.R

 CreateCalibOptions <- function(FUN_MOD,
                                FUN_CALIB = Calibration_Michel,
                                FUN_TRANSFO = NULL,
+                               IsHyst = FALSE,
                                FixedParam = NULL,
                                SearchRanges = NULL,
                                StartParamList = NULL,
@@ -13,7 +14,9 @@ CreateCalibOptions <- function(FUN_MOD,
     if(!is.null(FUN_TRANSFO)) {
       FUN_TRANSFO <- match.fun(FUN_TRANSFO)
     }
-    
+    if (!is.logical(IsHyst) | length(IsHyst) != 1L) {
+      stop("'IsHyst' must be a 'logical' of length 1")
+    }
     ##check_FUN_MOD
     BOOL <- FALSE
     
@@ -57,6 +60,9 @@ CreateCalibOptions <- function(FUN_MOD,
       ObjectClass <- c(ObjectClass, "CemaNeigeGR6J")
       BOOL <- TRUE
     }
+	if (IsHyst) {
+	  ObjectClass <- c(ObjectClass, "hysteresis")
+	}
     if (!BOOL) {
       stop("incorrect FUN_MOD for use in CreateCalibOptions")
       return(NULL)
@@ -100,7 +106,7 @@ CreateCalibOptions <- function(FUN_MOD,
         FUN1 <- TransfoParam_GR1A
       }
 	  if (identical(FUN_MOD, RunModel_CemaNeige)) {
-        if (inherits(FUN_MOD, "hysteresis")) {
+        if (IsHyst) {
           FUN1 <- TransfoParam_CemaNeigeHyst
         } else {
           FUN1 <- TransfoParam_CemaNeige
@@ -111,30 +117,50 @@ CreateCalibOptions <- function(FUN_MOD,
         return(NULL)
       }
       ##_set_FUN2
-      FUN2 <- TransfoParam_CemaNeige
-      
+      if (IsHyst) {
+        FUN2 <- TransfoParam_CemaNeigeHyst
+      } else {
+        FUN2 <- TransfoParam_CemaNeige
+      }		
       ##_set_FUN_TRANSFO
       if (sum(ObjectClass %in% c("GR4H", "GR4J", "GR5J", "GR6J", "GR2M", "GR1A", "CemaNeige")) > 0) {
         FUN_TRANSFO <- FUN1
       } else {
-        FUN_TRANSFO <- function(ParamIn, Direction) {
-          Bool <- is.matrix(ParamIn)
-          if (Bool == FALSE) {
-            ParamIn <- rbind(ParamIn)
+	     if (IsHyst) {
+           FUN_TRANSFO <- function(ParamIn, Direction) {
+             Bool <- is.matrix(ParamIn)
+             if (Bool == FALSE) {
+               ParamIn <- rbind(ParamIn)
+             }
+             ParamOut <- NA * ParamIn
+             NParam   <- ncol(ParamIn)
+             ParamOut[, 1:(NParam - 4)] <- FUN1(ParamIn[, 1:(NParam - 4)], Direction)
+             ParamOut[, (NParam - 3):NParam] <- FUN2(ParamIn[, (NParam - 3):NParam], Direction)
+             if (Bool == FALSE) {
+               ParamOut <- ParamOut[1, ]
+             }
+             return(ParamOut)
           }
-          ParamOut <- NA * ParamIn
-          NParam   <- ncol(ParamIn)
-          if (NParam <= 3) {
-            ParamOut[, 1:(NParam - 2)] <- FUN1(cbind(ParamIn[, 1:(NParam - 2)]), Direction)
-          } else {
-            ParamOut[, 1:(NParam - 2)] <- FUN1(ParamIn[, 1:(NParam - 2)], Direction)
-          }
-          ParamOut[, (NParam - 1):NParam] <- FUN2(ParamIn[, (NParam - 1):NParam], Direction)
-          if (Bool == FALSE) {
-            ParamOut <- ParamOut[1, ]
-          }
-          return(ParamOut)
-        }
+		} else {
+          FUN_TRANSFO <- function(ParamIn, Direction) {
+            Bool <- is.matrix(ParamIn)
+            if (Bool == FALSE) {
+              ParamIn <- rbind(ParamIn)
+            }
+            ParamOut <- NA * ParamIn
+            NParam   <- ncol(ParamIn)
+            if (NParam <= 3) {
+              ParamOut[, 1:(NParam - 2)] <- FUN1(cbind(ParamIn[, 1:(NParam - 2)]), Direction)
+            } else {
+              ParamOut[, 1:(NParam - 2)] <- FUN1(ParamIn[, 1:(NParam - 2)], Direction)
+            }
+            ParamOut[, (NParam - 1):NParam] <- FUN2(ParamIn[, (NParam - 1):NParam], Direction)
+            if (Bool == FALSE) {
+              ParamOut <- ParamOut[1, ]
+            }
+            return(ParamOut)
+		  }
+		}
       }
     }
     if (is.null(FUN_TRANSFO)) {
@@ -173,10 +199,9 @@ CreateCalibOptions <- function(FUN_MOD,
     if ("CemaNeigeGR6J" %in% ObjectClass) {
       NParam <- 8
     }
-	if (inherits(FUN_MOD, "hysteresis")) {
+	if (IsHyst) {
 	  NParam <- NParam + 2
 	}
-
     
     ##check_FixedParam
     if (is.null(FixedParam)) {
@@ -200,7 +225,6 @@ CreateCalibOptions <- function(FUN_MOD,
     if (is.null(SearchRanges)) {
       ParamT <-  matrix(c(rep(-9.99, NParam), rep(+9.99, NParam)),
                         ncol = NParam, byrow = TRUE)
-      
       SearchRanges <- TransfoParam(ParamIn = ParamT, Direction = "TR", FUN_TRANSFO = FUN_TRANSFO)
       
     } else {
@@ -297,7 +321,7 @@ CreateCalibOptions <- function(FUN_MOD,
                              +4.50, +0.50, +5.00, -8.10, +1.10, +5.00, +4.10, +7.21), ncol = 8, byrow = TRUE)
         }
 	  # }
-      if (inherits(FUN_MOD, "hysteresis")) {
+      if (IsHyst) {
         ParamTHyst <- matrix(c(-9.08, -6.99,
                                -8.00, -3.20,
                                -6.40, +9.99), ncol = 2, byrow = TRUE)

man/CreateCalibOptions.Rd

@@ -10,7 +10,7 @@
 
 \usage{
 CreateCalibOptions(FUN_MOD, FUN_CALIB = Calibration_Michel,
-  FUN_TRANSFO = NULL, FixedParam = NULL,
+  FUN_TRANSFO = NULL, IsHyst = FALSE, FixedParam = NULL,
   SearchRanges = NULL, StartParamList = NULL,
   StartParamDistrib = NULL)
 }
@@ -23,6 +23,8 @@ CreateCalibOptions(FUN_MOD, FUN_CALIB = Calibration_Michel,
 
 \item{FUN_TRANSFO}{(optional) [function] model parameters transformation function, if the \code{FUN_MOD} used is native in the package, \code{FUN_TRANSFO} is automatically defined}
 
+\item{IsHyst}{[boolean] boolean indicating if the hysteresis version of CemaNeige is used. See details}
+
 \item{FixedParam}{(optional) [numeric] vector giving the values set for the non-optimised parameter values (NParam columns, 1 line)
 \cr Example:
 \tabular{llllll}{
@@ -75,6 +77,12 @@ Creation of the \emph{CalibOptions} object required by the \code{Calibration*} f
 \details{
 Users wanting to use \code{FUN_MOD}, \code{FUN_CALIB} or \code{FUN_TRANSFO} functions that are not included in 
 the package must create their own \code{CalibOptions} object accordingly.
+
+##### CemaNeige version #####
+
+If \code{IsHyst = FALSE}, the original CemaNeige version from Valéry et al. (2014) is used.  \cr
+If \code{IsHyst = TRUE}, the CemaNeige version from Riboust et al. (2019) is used. Compared to the original version, this version of CemaNeige needs two more parameters and it includes a representation of the hysteretic relationship between the Snow Cover Area (SCA) and the Snow Water Equivalent (SWE) in the catchment. The hysteresis included in airGR is the Modified Linear hysteresis (LH*); it is represented on panel b) of Fig. 3 in Riboust et al. (2019). Riboust et al. (2019) advise to use the LH* version of CemaNeige with parameters calibrated using an objective function combining 75 \% of KGE calculated on discharge simulated from a rainfall-runoff model compared to observed discharge and 5 \% of KGE calculated on SCA on 5 CemaNeige elevation bands compared to satellite (e.g. MODIS) SCA (see Eq. (18), Table 3 and Fig. 6). Riboust et al. (2019)'s tests were realized with GR4J as the chosen rainfall-runoff model. \cr \cr
+
 }