diff --git a/R/CreateCalibOptions.R b/R/CreateCalibOptions.R
index fc71940be53f63ee7ff9d7b93ccf398e1b434ce3..d944d9dc7954eec2567409b17a2e4cbebf6ec3b2 100644
--- a/R/CreateCalibOptions.R
+++ b/R/CreateCalibOptions.R
@@ -23,17 +23,21 @@ CreateCalibOptions <- function(x, ...) {
#' @export
CreateCalibOptions.InputsModel <- function(x,
...) {
- if (!exists("FUN_MOD") && !is.null(x$FUN_MOD)) {
- airGR::CreateCalibOptions(
- FUN_MOD = x$FUN_MOD,
- IsSD = !is.null(x$Qupstream) & x$FUN_MOD != "RunModel_Lag",
- ...
- )
- } else {
- airGR::CreateCalibOptions(
- ...
- )
+ dots <- list(...)
+ # Add FUN_MOD in parameters if carried by InputsModel
+ if (!"FUN_MOD" %in% names(dots)) {
+ if(!is.null(x$FUN_MOD)) {
+ dots$FUN_MOD <- x$FUN_MOD
+ } else {
+ stop(" The parameter `FUN_MOD` must be defined")
+ }
}
+ # Automatically define IsSD for intermediate basin GR models
+ dots$IsSD = !is.null(x$Qupstream) & dots$FUN_MOD != "RunModel_Lag"
+ # Add IsHyst in parameters if carried by InputsModel
+ if (!is.null(x$model$IsHyst)) dots$IsHyst <- x$model$IsHyst
+ # Call airGR function
+ do.call(airGR::CreateCalibOptions, dots)
}
#' @rdname CreateCalibOptions
diff --git a/R/CreateInputsModel.GRiwrm.R b/R/CreateInputsModel.GRiwrm.R
index 50edf1fbd53f26d05818c1278dba3c69eaa90001..115916133134eee00fa8783bdb8234d88c67e920 100644
--- a/R/CreateInputsModel.GRiwrm.R
+++ b/R/CreateInputsModel.GRiwrm.R
@@ -35,6 +35,8 @@
#' @param NLayers (optional) named [vector] of [numeric] integer giving the number
#' of elevation layers requested [-], required to create CemaNeige module
#' inputs, default=5
+#' @param IsHyst [logical] boolean indicating if the hysteresis version of
+#' CemaNeige is used. See details of [airGR::CreateRunOptions].
#' @param ... used for compatibility with S3 methods
#'
#' @details Meteorological data are needed for the nodes of the network that
@@ -70,7 +72,8 @@ CreateInputsModel.GRiwrm <- function(x, DatesR,
PrecipScale = TRUE,
TempMean = NULL, TempMin = NULL,
TempMax = NULL, ZInputs = NULL,
- HypsoData = NULL, NLayers = 5, ...) {
+ HypsoData = NULL, NLayers = 5,
+ IsHyst = FALSE, ...) {
# Check and format inputs
varNames <- c("Precip", "PotEvap", "TempMean", "Qobs", "Qmin",
@@ -212,7 +215,8 @@ CreateInputsModel.GRiwrm <- function(x, DatesR,
HypsoData = getInputBV(HypsoData, id),
NLayers = getInputBV(NLayers, id, 5),
Qobs = Qobs,
- Qmin = getInputBV(Qmin, id)
+ Qmin = getInputBV(Qmin, id),
+ IsHyst = IsHyst
)
}
attr(InputsModel, "TimeStep") <- getModelTimeStep(InputsModel)
@@ -247,7 +251,7 @@ CreateEmptyGRiwrmInputsModel <- function(griwrm) {
#'
#' @return \emph{InputsModel} object for one.
#' @noRd
-CreateOneGRiwrmInputsModel <- function(id, griwrm, ..., Qobs, Qmin) {
+CreateOneGRiwrmInputsModel <- function(id, griwrm, ..., Qobs, Qmin, IsHyst) {
np <- getNodeProperties(id, griwrm)
if (np$Diversion) {
@@ -326,7 +330,7 @@ CreateOneGRiwrmInputsModel <- function(id, griwrm, ..., Qobs, Qmin) {
# Add the model function
InputsModel$FUN_MOD <- FUN_MOD
- featModel <- .GetFeatModel(InputsModel)
+ featModel <- .GetFeatModel(InputsModel, IsHyst)
InputsModel$isUngauged <- griwrm$model[griwrm$id == id] == "Ungauged"
InputsModel$gaugedId <- griwrm$donor[griwrm$id == id]
InputsModel$hasUngaugedNodes <- hasUngaugedNodes(id, griwrm)
@@ -334,8 +338,9 @@ CreateOneGRiwrmInputsModel <- function(id, griwrm, ..., Qobs, Qmin) {
list(
indexParamUngauged = ifelse(inherits(InputsModel, "SD"), 0, 1) +
seq.int(featModel$NbParam),
- hasX4 = grepl("RunModel_GR[456][HJ]", FUN_MOD),
- iX4 = ifelse(inherits(InputsModel, "SD"), 5, 4)
+ hasX4 = grepl("RunModel_(CemaNeige|)GR[456][HJ]", FUN_MOD),
+ iX4 = ifelse(inherits(InputsModel, "SD"), 5, 4),
+ IsHyst = featModel$IsHyst
)
InputsModel$hasDiversion <- np$Diversion
InputsModel$isReservoir <- np$Reservoir
@@ -446,7 +451,7 @@ hasUngaugedNodes <- function(id, griwrm) {
#' function to extract model features partially copied from airGR:::.GetFeatModel
#' @importFrom utils tail
#' @noRd
-.GetFeatModel <- function(InputsModel) {
+.GetFeatModel <- function(InputsModel, IsHyst) {
path <- system.file("modelsFeatures/FeatModelsGR.csv", package = "airGR")
FeatMod <- read.table(path, header = TRUE, sep = ";", stringsAsFactors = FALSE)
NameFunMod <- ifelse(test = FeatMod$Pkg %in% "airGR",
@@ -461,6 +466,13 @@ hasUngaugedNodes <- function(id, griwrm) {
if (FeatMod$IsSD) {
FeatMod$NbParam <- FeatMod$NbParam + 1
}
+ FeatMod$IsHyst <- FALSE
+ if (grepl("CemaNeige", FeatMod$NameMod)) {
+ FeatMod$IsHyst <- IsHyst
+ if (IsHyst) {
+ FeatMod$NbParam <- FeatMod$NbParam + 2
+ }
+ }
return(FeatMod)
}
diff --git a/R/CreateRunOptions.R b/R/CreateRunOptions.R
index e6762391d008c9c850e558160ca951d6ca8f56c7..c74cc68539bd87d676e1d2428b9cc1ddd63c239a 100644
--- a/R/CreateRunOptions.R
+++ b/R/CreateRunOptions.R
@@ -43,6 +43,8 @@ CreateRunOptions.InputsModel <- function(x, ...) {
stop(" The parameter `FUN_MOD` must be defined")
}
}
+ # Add IsHyst in parameters if carried by InputsModel
+ if (!is.null(x$model$IsHyst)) dots$IsHyst <- x$model$IsHyst
# Temporary fix waiting for resolution of HYCAR-Hydro/airgr#167
if (identical(match.fun(dots$FUN_MOD), RunModel_Lag)) {
diff --git a/man/CreateInputsModel.GRiwrm.Rd b/man/CreateInputsModel.GRiwrm.Rd
index f77b5cebec78cbbdab61ae6f8839cfb17ede6077..26b9a66bc318a276b31f98ce25216ad896e2e9db 100644
--- a/man/CreateInputsModel.GRiwrm.Rd
+++ b/man/CreateInputsModel.GRiwrm.Rd
@@ -18,6 +18,7 @@
ZInputs = NULL,
HypsoData = NULL,
NLayers = 5,
+ IsHyst = FALSE,
...
)
}
@@ -70,6 +71,9 @@ if not defined a single elevation is used for CemaNeige}
of elevation layers requested \link{-}, required to create CemaNeige module
inputs, default=5}
+\item{IsHyst}{\link{logical} boolean indicating if the hysteresis version of
+CemaNeige is used. See details of \link[airGR:CreateRunOptions]{airGR::CreateRunOptions}.}
+
\item{...}{used for compatibility with S3 methods}
}
\value{
diff --git a/tests/testthat/helper_1_RunModel.R b/tests/testthat/helper_1_RunModel.R
index 565d840b3609b927d2f389ba54df9d02106b9e62..3e9acbede60e83b4f6ad62b24b4d4f1a59ce50fb 100644
--- a/tests/testthat/helper_1_RunModel.R
+++ b/tests/testthat/helper_1_RunModel.R
@@ -15,7 +15,8 @@ setupRunModel <-
runRunOptions = TRUE,
runRunModel = TRUE,
griwrm = NULL,
- Qobs2 = NULL) {
+ Qobs2 = NULL,
+ IsHyst = FALSE) {
data(Severn)
@@ -43,12 +44,20 @@ setupRunModel <-
# Convert meteo data to SD (remove upstream areas)
Precip <- ConvertMeteoSD(griwrm, PrecipTot)
PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
+ if (IsHyst) {
+ TempMean <- PotEvap+5 # Fake temperatures
+ } else {
+ TempMean <- NULL
+ }
# set up inputs
if (!runInputsModel)
return(environment())
InputsModel <-
- suppressWarnings(CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs = Qobs2))
+ suppressWarnings(CreateInputsModel(griwrm, DatesR, Precip, PotEvap,
+ TempMean = TempMean,
+ Qobs = Qobs2,
+ IsHyst = IsHyst))
# RunOptions
if (!runRunOptions)
diff --git a/tests/testthat/helper_Calibration.R b/tests/testthat/helper_Calibration.R
index c2928118c67bc5c677964329fa5947a110a2e6bb..08dc63672e9cff31494c3ebff5f344520423ed00 100644
--- a/tests/testthat/helper_Calibration.R
+++ b/tests/testthat/helper_Calibration.R
@@ -1,14 +1,19 @@
runCalibration <- function(nodes = NULL,
Qobs2 = NULL,
InputsCrit = NULL,
+ CalibOptions = NULL,
FUN_CRIT = ErrorCrit_KGE2,
- runRunModel = FALSE) {
+ runRunModel = FALSE,
+ IsHyst = FALSE) {
if (is.null(nodes)) {
griwrm <- NULL
} else {
griwrm <- CreateGRiwrm(nodes)
}
- e <- setupRunModel(griwrm = griwrm, runRunModel = runRunModel, Qobs2 = Qobs2)
+ e <- setupRunModel(griwrm = griwrm,
+ runRunModel = runRunModel,
+ Qobs2 = Qobs2,
+ IsHyst = IsHyst)
for(x in ls(e)) assign(x, get(x, e))
rm(e)
np <- getAllNodesProperties(griwrm)
@@ -22,7 +27,8 @@ runCalibration <- function(nodes = NULL,
)
}
- CalibOptions <- CreateCalibOptions(InputsModel)
+ if (is.null(CalibOptions)) CalibOptions <- CreateCalibOptions(InputsModel)
+
OutputsCalib <- Calibration(InputsModel, RunOptions, InputsCrit, CalibOptions)
Param <- sapply(OutputsCalib, "[[", "ParamFinalR")
return(environment())
diff --git a/tests/testthat/test-RunModel_Ungauged.R b/tests/testthat/test-RunModel_Ungauged.R
index 1a03a4fd329df7e13dfe3432a35fe881dad652f0..5b6dd7bc0ea885da64c547ceb9fad707a7bec414 100644
--- a/tests/testthat/test-RunModel_Ungauged.R
+++ b/tests/testthat/test-RunModel_Ungauged.R
@@ -233,3 +233,59 @@ test_that("Donor node with diversion should work", {
for(x in ls(e)) assign(x, get(x, e))
expect_equal(OC_ref$`54032`$CritFinal, OutputsCalib$`54032`$CritFinal, tolerance = 1E-3)
})
+
+test_that("Cemaneige with hysteresis works", {
+ nodes <- loadSevernNodes()
+ nodes <- nodes[nodes$id %in% c("54057", "54032", "54001"), ]
+ nodes$model <- "RunModel_CemaNeigeGR4J"
+ nodes$model[nodes$id != 54057] <- "Ungauged"
+ griwrm <- CreateGRiwrm(nodes)
+
+ # # The custom ErrorCrit function !!!
+ ErrorCrit_KGE3 <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
+ OutputsCritQ <- suppressMessages(
+ ErrorCrit_KGE2(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE)
+ )
+ InputsCrit$Obs <- InputsCrit$SCA #a adapter
+ OutputsModel$Qsim <- OutputsModel$CemaNeigeLayers[[1]]$Gratio #a adapter
+ OutputsCritSCA <- suppressMessages(
+ ErrorCrit_KGE2(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE)
+ )
+ OutputsCritQ$CritValue <-
+ (OutputsCritQ$CritValue + OutputsCritSCA$CritValue) / 2
+ OutputsCritQ$CritName <- "(0.5 * KGE2[Q] + 0.5 * KGE2[SCA])"
+ return(OutputsCritQ)
+ }
+ class(ErrorCrit_KGE3) <- c("FUN_CRIT", class(ErrorCrit_KGE3))
+
+ e <- suppressWarnings(
+ setupRunModel(griwrm = griwrm, runRunModel = FALSE, IsHyst = TRUE)
+ )
+ for(x in ls(e)) assign(x, get(x, e))
+
+ expect_true(all(sapply(InputsModel, function(x) x$model$hasX4)))
+
+ np <- getAllNodesProperties(griwrm)
+ InputsCrit <- CreateInputsCrit(
+ InputsModel,
+ FUN_CRIT = ErrorCrit_KGE3,
+ RunOptions = RunOptions,
+ Obs = Qobs[IndPeriod_Run, np$id[np$RunOff & np$calibration == "Gauged"], drop = FALSE],
+ )
+ InputsCrit$`54057`$SCA <- runif(length(IndPeriod_Run)) # Fake SCA
+ CalibOptions <- CreateCalibOptions(InputsModel)
+ CO <- lapply(CalibOptions, function(x) {
+ x$StartParamList <- matrix(
+ c(0.605, 320.596, -0.042, 37.991, 2.221, 0.705, 6.764, 85.000, 0.850),
+ nrow = 1)
+ x$StartParamDistrib <- NULL
+ x
+ })
+ class(CO) <- class(CalibOptions)
+ e <- suppressWarnings(
+ runCalibration(nodes, InputsCrit = InputsCrit, CalibOptions = CO, IsHyst = TRUE)
+ )
+ for(x in ls(e)) assign(x, get(x, e))
+ expect_equal(sapply(Param, length),
+ c("54057" = 9, "54032" = 9, "54001" = 8))
+})