From db0de232ba243e1f2b22076181cae3d78dea7691 Mon Sep 17 00:00:00 2001
From: Dorchies David <david.dorchies@irstea.fr>
Date: Mon, 5 Oct 2020 12:09:56 +0200
Subject: [PATCH] v1.6.3.1 feat(SD): Change order of param for LAG
- The lag parameter is now in first position instead of the last one
Refs #34, #60
---
DESCRIPTION | 10 ++--
NEWS.md | 9 +++
R/Calibration_Michel.R | 16 +++---
R/CreateCalibOptions.R | 90 +++++++++++++++---------------
R/RunModel.R | 56 ++++---------------
R/RunModel_LAG.R | 43 ++++++++++++++
tests/testthat/test-RunModel_LAG.R | 12 ++--
7 files changed, 128 insertions(+), 108 deletions(-)
create mode 100644 R/RunModel_LAG.R
diff --git a/DESCRIPTION b/DESCRIPTION
index 44ca4908..5f014f72 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,23 +1,23 @@
Package: airGR
Type: Package
Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
-Version: 1.6.2.3
-Date: 2020-07-16
+Version: 1.6.3.1
+Date: 2020-10-05
Authors@R: c(
person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
person("Olivier", "Delaigue", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-7668-8468"), email = "airGR@inrae.fr"),
person("Guillaume", "Thirel", role = c("aut"), comment = c(ORCID = "0000-0002-1444-1830")),
person("Charles", "Perrin", role = c("aut", "ths"), comment = c(ORCID = "0000-0001-8552-1881")),
person("Claude", "Michel", role = c("aut", "ths")),
- person("Vazken", "Andréassian", role = c("ctb", "ths"), comment = c(ORCID = "0000-0001-7124-9303")),
- person("François", "Bourgin", role = c("ctb"), comment = c(ORCID = "0000-0002-2820-7260", vignette = "'Parameter estimation' vignettes")),
+ person("Vazken", "Andréassian", role = c("ctb", "ths"), comment = c(ORCID = "0000-0001-7124-9303")),
+ person("François", "Bourgin", role = c("ctb"), comment = c(ORCID = "0000-0002-2820-7260", vignette = "'Parameter estimation' vignettes")),
person("Pierre", "Brigode", role = c("ctb"), comment = c(ORCID = "0000-0001-8257-0741")),
person("Nicolas", "Le Moine", role = c("ctb")),
person("Thibaut", "Mathevet", role = c("ctb"), comment = c(ORCID = "0000-0002-4142-4454")),
person("Safouane", "Mouelhi", role = c("ctb")),
person("Ludovic", "Oudin", role = c("ctb"), comment = c(ORCID = "0000-0002-3712-0933")),
person("Raji", "Pushpalatha", role = c("ctb")),
- person("Audrey", "Valéry", role = c("ctb"))
+ person("Audrey", "Valéry", role = c("ctb"))
)
Depends: R (>= 3.0.1)
Suggests: knitr, rmarkdown, coda, DEoptim, dplyr, FME, ggmcmc, hydroPSO, Rmalschains, testthat
diff --git a/NEWS.md b/NEWS.md
index d3b32a32..842f8784 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,14 @@
## Release History of the airGR Package
+### 1.6.3.1 Release Notes (2020-10-05)
+
+#### New features
+
+- Change order of parameters: LAG is now the first parameter instead of the last
+
+____________________________________________________________________________________
+
+
### 1.6.2.2 Release Notes (2020-07-15)
#### Bug fixes
diff --git a/R/Calibration_Michel.R b/R/Calibration_Michel.R
index 9be0e426..11a50103 100644
--- a/R/Calibration_Michel.R
+++ b/R/Calibration_Michel.R
@@ -101,8 +101,8 @@ Calibration_Michel <- function(InputsModel,
}
ParamOut <- NA * ParamIn
NParam <- ncol(ParamIn)
- ParamOut[, 1:(NParam-1)] <- FUN1(ParamIn[, 1:(NParam-1)], Direction)
- ParamOut[, NParam ] <- FUN_LAG(as.matrix(ParamIn[, NParam ]), Direction)
+ ParamOut[, 2:NParam] <- FUN1(ParamIn[, 2:NParam], Direction)
+ ParamOut[, 1 ] <- FUN_LAG(as.matrix(ParamIn[, 1]), Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
}
@@ -157,9 +157,9 @@ Calibration_Michel <- function(InputsModel,
}
ParamOut <- NA * ParamIn
NParam <- ncol(ParamIn)
- ParamOut[, 1:(NParam-5)] <- FUN1( ParamIn[, 1:(NParam-5)], Direction)
- ParamOut[, (NParam-4):(NParam-1)] <- FUN2( ParamIn[, (NParam-4):(NParam-1)], Direction)
- ParamOut[, NParam ] <- FUN_LAG(as.matrix(ParamIn[, NParam ]), Direction)
+ ParamOut[, 2:(NParam-4)] <- FUN1( ParamIn[, 2:(NParam-4)], Direction)
+ ParamOut[, (NParam-3):NParam] <- FUN2( ParamIn[, (NParam-3):NParam], Direction)
+ ParamOut[, 1 ] <- FUN_LAG(as.matrix(ParamIn[, 1 ]), Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
}
@@ -174,9 +174,9 @@ Calibration_Michel <- function(InputsModel,
}
ParamOut <- NA * ParamIn
NParam <- ncol(ParamIn)
- ParamOut[, 1:(NParam-3)] <- FUN1( ParamIn[, 1:(NParam-3)], Direction)
- ParamOut[, (NParam-2):(NParam-1)] <- FUN2( ParamIn[, (NParam-2):(NParam-1)], Direction)
- ParamOut[, NParam ] <- FUN_LAG(as.matrix(ParamIn[, NParam ]), Direction)
+ ParamOut[, 2:(NParam-2)] <- FUN1( ParamIn[, 2:(NParam-2)], Direction)
+ ParamOut[, (NParam-1):NParam] <- FUN2( ParamIn[, (NParam-1):NParam], Direction)
+ ParamOut[, 1 ] <- FUN_LAG(as.matrix(ParamIn[, 1 ]), Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
}
diff --git a/R/CreateCalibOptions.R b/R/CreateCalibOptions.R
index 09ac45f8..422d2c18 100644
--- a/R/CreateCalibOptions.R
+++ b/R/CreateCalibOptions.R
@@ -7,9 +7,9 @@ CreateCalibOptions <- function(FUN_MOD,
SearchRanges = NULL,
StartParamList = NULL,
StartParamDistrib = NULL) {
-
+
ObjectClass <- NULL
-
+
FUN_MOD <- match.fun(FUN_MOD)
FUN_CALIB <- match.fun(FUN_CALIB)
if(!is.null(FUN_TRANSFO)) {
@@ -23,7 +23,7 @@ CreateCalibOptions <- function(FUN_MOD,
}
##check_FUN_MOD
BOOL <- FALSE
-
+
if (identical(FUN_MOD, RunModel_GR4H)) {
ObjectClass <- c(ObjectClass, "GR4H")
BOOL <- TRUE
@@ -59,7 +59,7 @@ CreateCalibOptions <- function(FUN_MOD,
if (identical(FUN_MOD, RunModel_CemaNeigeGR4H)) {
ObjectClass <- c(ObjectClass, "CemaNeigeGR4H")
BOOL <- TRUE
- }
+ }
if (identical(FUN_MOD, RunModel_CemaNeigeGR5H)) {
ObjectClass <- c(ObjectClass, "CemaNeigeGR5H")
BOOL <- TRUE
@@ -86,10 +86,10 @@ CreateCalibOptions <- function(FUN_MOD,
stop("incorrect 'FUN_MOD' for use in 'CreateCalibOptions'")
return(NULL)
}
-
+
##check_FUN_CALIB
BOOL <- FALSE
-
+
if (identical(FUN_CALIB, Calibration_Michel)) {
ObjectClass <- c(ObjectClass, "HBAN")
BOOL <- TRUE
@@ -97,9 +97,9 @@ CreateCalibOptions <- function(FUN_MOD,
if (!BOOL) {
stop("incorrect 'FUN_CALIB' for use in 'CreateCalibOptions'")
return(NULL)
-
+
}
-
+
##check_FUN_TRANSFO
if (is.null(FUN_TRANSFO)) {
##_set_FUN1
@@ -145,11 +145,11 @@ CreateCalibOptions <- function(FUN_MOD,
FUN2 <- TransfoParam_CemaNeigeHyst
} else {
FUN2 <- TransfoParam_CemaNeige
- }
+ }
##_set_FUN_LAG
if (IsSD) {
FUN_LAG <- TransfoParam_LAG
- }
+ }
##_set_FUN_TRANSFO
if (sum(ObjectClass %in% c("GR4H", "GR5H", "GR4J", "GR5J", "GR6J", "GR2M", "GR1A", "CemaNeige")) > 0) {
if (!IsSD) {
@@ -162,14 +162,14 @@ CreateCalibOptions <- function(FUN_MOD,
}
ParamOut <- NA * ParamIn
NParam <- ncol(ParamIn)
- ParamOut[, 1:(NParam-1)] <- FUN1(ParamIn[, 1:(NParam-1)], Direction)
- ParamOut[, NParam ] <- FUN_LAG(as.matrix(ParamIn[, NParam ]), Direction)
+ ParamOut[, 2:NParam] <- FUN1(ParamIn[, 2:NParam], Direction)
+ ParamOut[, 1 ] <- FUN_LAG(as.matrix(ParamIn[, 1]), Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
}
return(ParamOut)
}
- }
+ }
} else {
if (IsHyst & !IsSD) {
FUN_TRANSFO <- function(ParamIn, Direction) {
@@ -215,9 +215,9 @@ CreateCalibOptions <- function(FUN_MOD,
}
ParamOut <- NA * ParamIn
NParam <- ncol(ParamIn)
- ParamOut[, 1:(NParam - 5) ] <- FUN1( ParamIn[, 1:(NParam - 5) ], Direction)
- ParamOut[, (NParam - 4):(NParam - 1)] <- FUN2( ParamIn[, (NParam - 4):(NParam - 1)], Direction)
- ParamOut[, NParam ] <- FUN_LAG(as.matrix(ParamIn[, NParam ]), Direction)
+ ParamOut[, 2:(NParam - 4) ] <- FUN1( ParamIn[, 2:(NParam - 4) ], Direction)
+ ParamOut[, (NParam - 3):NParam] <- FUN2( ParamIn[, (NParam - 3):NParam], Direction)
+ ParamOut[, 1 ] <- FUN_LAG(as.matrix(ParamIn[, 1 ]), Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
}
@@ -233,12 +233,12 @@ CreateCalibOptions <- function(FUN_MOD,
ParamOut <- NA * ParamIn
NParam <- ncol(ParamIn)
if (NParam <= 3) {
- ParamOut[, 1:(NParam - 3)] <- FUN1(cbind(ParamIn[, 1:(NParam - 3)]), Direction)
+ ParamOut[, 2:(NParam - 2)] <- FUN1(cbind(ParamIn[, 2:(NParam - 2)]), Direction)
} else {
- ParamOut[, 1:(NParam - 3)] <- FUN1( ParamIn[, 1:(NParam - 3)], Direction)
+ ParamOut[, 2:(NParam - 2)] <- FUN1( ParamIn[, 2:(NParam - 2)], Direction)
}
- ParamOut[, (NParam - 2):(NParam - 1)] <- FUN2( ParamIn[, (NParam - 2):(NParam - 1)], Direction)
- ParamOut[, NParam ] <- FUN_LAG(as.matrix(ParamIn[, NParam]), Direction)
+ ParamOut[, (NParam - 1):NParam] <- FUN2( ParamIn[, (NParam - 1):NParam], Direction)
+ ParamOut[, 1 ] <- FUN_LAG(as.matrix(ParamIn[, 1]), Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
}
@@ -251,7 +251,7 @@ CreateCalibOptions <- function(FUN_MOD,
stop("'FUN_TRANSFO' was not found")
return(NULL)
}
-
+
##NParam
if ("GR4H" %in% ObjectClass) {
NParam <- 4
@@ -279,10 +279,10 @@ CreateCalibOptions <- function(FUN_MOD,
}
if ("CemaNeigeGR4H" %in% ObjectClass) {
NParam <- 6
- }
+ }
if ("CemaNeigeGR5H" %in% ObjectClass) {
NParam <- 7
- }
+ }
if ("CemaNeigeGR4J" %in% ObjectClass) {
NParam <- 6
}
@@ -298,7 +298,7 @@ CreateCalibOptions <- function(FUN_MOD,
if (IsSD) {
NParam <- NParam + 1
}
-
+
##check_FixedParam
if (is.null(FixedParam)) {
FixedParam <- rep(NA, NParam)
@@ -316,13 +316,13 @@ CreateCalibOptions <- function(FUN_MOD,
warning("You have not set any parameter in 'FixedParam'")
}
}
-
+
##check_SearchRanges
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 {
if (!is.matrix(SearchRanges)) {
stop("'SearchRanges' must be a matrix")
@@ -340,7 +340,7 @@ CreateCalibOptions <- function(FUN_MOD,
stop("Incompatibility between 'SearchRanges' ncol and 'FUN_MOD'")
}
}
-
+
##check_StartParamList_and_StartParamDistrib__default_values
if (("HBAN" %in% ObjectClass & is.null(StartParamList) & is.null(StartParamDistrib))) {
if ("GR4H" %in% ObjectClass) {
@@ -356,7 +356,7 @@ CreateCalibOptions <- function(FUN_MOD,
if (("GR5H" %in% ObjectClass) & !("interception" %in% ObjectClass)) {
ParamT <- matrix(c(+3.28, -0.39, +4.14, -9.54, -7.49,
+3.62, -0.19, +4.80, -9.00, -6.31,
- +4.01, -0.04, +5.43, -7.53, -5.33), ncol=5, byrow = TRUE);
+ +4.01, -0.04, +5.43, -7.53, -5.33), ncol=5, byrow = TRUE);
}
if ("GR4J" %in% ObjectClass) {
ParamT <- matrix(c(+5.13, -1.60, +3.03, -9.05,
@@ -367,7 +367,7 @@ CreateCalibOptions <- function(FUN_MOD,
ParamT <- matrix(c(+5.17, -1.13, +3.08, -9.37, -7.45,
+5.55, -0.46, +3.75, -9.09, -4.69,
+6.10, -0.11, +4.43, -8.60, -0.66), ncol = 5, byrow = TRUE)
-
+
}
if ("GR6J" %in% ObjectClass) {
ParamT <- matrix(c(+3.60, -1.00, +3.30, -9.10, -0.90, +3.00,
@@ -384,8 +384,8 @@ CreateCalibOptions <- function(FUN_MOD,
-0.38,
+1.39), ncol = 1, byrow = TRUE)
}
-
-
+
+
if ("CemaNeige" %in% ObjectClass) {
ParamT <- matrix(c(-9.96, +6.63,
-9.14, +6.90,
@@ -395,7 +395,7 @@ CreateCalibOptions <- function(FUN_MOD,
ParamT <- matrix(c(+5.12, -1.18, +4.34, -9.69, -9.96, +6.63,
+5.58, -0.85, +4.74, -9.47, -9.14, +6.90,
+6.01, -0.50, +5.14, -8.87, +4.10, +7.21), ncol = 6, byrow = TRUE)
- }
+ }
if (("CemaNeigeGR5H" %in% ObjectClass) & ("interception" %in% ObjectClass)) {
ParamT <- matrix(c(+3.46, -1.25, +4.04, -9.53, -9.34, -9.96, +6.63,
+3.74, -0.41, +4.78, -8.94, -3.33, -9.14, +6.90,
@@ -404,7 +404,7 @@ CreateCalibOptions <- function(FUN_MOD,
if (("CemaNeigeGR5H" %in% ObjectClass) & !("interception" %in% ObjectClass)) {
ParamT <- matrix(c(+3.28, -0.39, +4.14, -9.54, -7.49, -9.96, +6.63,
+3.62, -0.19, +4.80, -9.00, -6.31, -9.14, +6.90,
- +4.01, -0.04, +5.43, -7.53, -5.33, +4.10, +7.21), ncol = 7, byrow = TRUE);
+ +4.01, -0.04, +5.43, -7.53, -5.33, +4.10, +7.21), ncol = 7, byrow = TRUE);
}
if ("CemaNeigeGR4J" %in% ObjectClass) {
ParamT <- matrix(c(+5.13, -1.60, +3.03, -9.05, -9.96, +6.63,
@@ -421,7 +421,7 @@ CreateCalibOptions <- function(FUN_MOD,
+3.90, -0.50, +4.10, -8.70, +0.10, +4.00, -9.14, +6.90,
+4.50, +0.50, +5.00, -8.10, +1.10, +5.00, +4.10, +7.21), ncol = 8, byrow = TRUE)
}
-
+
if (IsHyst) {
ParamTHyst <- matrix(c(-7.00, -7.00,
-0.00, -0.00,
@@ -429,17 +429,17 @@ CreateCalibOptions <- function(FUN_MOD,
ParamT <- cbind(ParamT, ParamTHyst)
}
if (IsSD) {
- ParamTSD <- matrix(c(+1.25,
- +2.50,
+ ParamTSD <- matrix(c(+1.25,
+ +2.50,
+5.00), ncol = 1, byrow = TRUE)
- ParamT <- cbind(ParamT, ParamTSD)
+ ParamT <- cbind(ParamTSD, ParamT)
}
-
+
StartParamList <- NULL
StartParamDistrib <- TransfoParam(ParamIn = ParamT, Direction = "TR", FUN_TRANSFO = FUN_TRANSFO)
-
+
}
-
+
##check_StartParamList_and_StartParamDistrib__format
if ("HBAN" %in% ObjectClass & !is.null(StartParamList)) {
if (!is.matrix(StartParamList)) {
@@ -469,11 +469,11 @@ CreateCalibOptions <- function(FUN_MOD,
stop("Incompatibility between 'StartParamDistrib' ncol and 'FUN_MOD'")
}
}
-
-
+
+
##Create_CalibOptions
CalibOptions <- list(FixedParam = FixedParam, SearchRanges = SearchRanges)
-
+
if (!is.null(StartParamList)) {
CalibOptions <- c(CalibOptions, list(StartParamList = StartParamList))
}
@@ -481,7 +481,7 @@ CreateCalibOptions <- function(FUN_MOD,
CalibOptions <- c(CalibOptions, list(StartParamDistrib = StartParamDistrib))
}
class(CalibOptions) <- c("CalibOptions", ObjectClass)
-
+
return(CalibOptions)
-
+
}
\ No newline at end of file
diff --git a/R/RunModel.R b/R/RunModel.R
index 9ca647e8..c4eeebc9 100644
--- a/R/RunModel.R
+++ b/R/RunModel.R
@@ -3,51 +3,19 @@ RunModel <- function (InputsModel, RunOptions, Param, FUN_MOD) {
FUN_MOD <- match.fun(FUN_MOD)
if (inherits(InputsModel, "SD")) {
- OutputsModelDown <- FUN_MOD(InputsModel = InputsModel, RunOptions = RunOptions,
- Param = Param[-length(Param)])
- OutputsModelDown$QsimDown <- OutputsModelDown$Qsim
-
- if (inherits(InputsModel, "daily")) {
- TimeStep <- 60 * 60 * 24
- }
- if (inherits(InputsModel, "hourly")) {
- TimeStep <- 60 * 60
- }
-
- # propagation time from upstream meshes to outlet
- PT <- InputsModel$LengthHydro / Param[length(Param)] / TimeStep
- HUTRANS <- rbind(1 - (PT - floor(PT)), PT - floor(PT))
-
- NbUpBasins <- length(InputsModel$LengthHydro)
- LengthTs <- length(OutputsModelDown$QsimDown)
- OutputsModelDown$Qsim <- OutputsModelDown$QsimDown * InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1E3
-
- for (upstream_basin in seq_len(NbUpBasins)) {
- Qupstream <- InputsModel$Qupstream[RunOptions$IndPeriod_Run, upstream_basin]
- if(!is.na(InputsModel$BasinAreas[upstream_basin])) {
- # Upstream flow with area needs to be converted to m3 by time step
- Qupstream <- Qupstream * InputsModel$BasinAreas[upstream_basin] * 1E3
- }
- OutputsModelDown$Qsim <- OutputsModelDown$Qsim +
- c(rep(0, floor(PT[upstream_basin])),
- Qupstream[1:(LengthTs - floor(PT[upstream_basin]))]) *
- HUTRANS[1, upstream_basin] +
- c(rep(0, floor(PT[upstream_basin] + 1)),
- Qupstream[1:(LengthTs - floor(PT[upstream_basin]) - 1)]) *
- HUTRANS[2, upstream_basin]
- }
- # Warning for negative flows
- if(any(OutputsModelDown$Qsim < 0)) {
- warning(length(which(OutputsModelDown$Qsim < 0)), " time steps with negative flow, set to zero.")
- OutputsModelDown$Qsim[OutputsModelDown$Qsim < 0] <- 0
- }
- # Convert back Qsim to mm
- OutputsModelDown$Qsim <- OutputsModelDown$Qsim / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1E3
-
+ # LAG Model take one parameter at the beginning of the vector
+ iFirstParamRunOffModel <- 2
} else {
+ # All parameters
+ iFirstParamRunOffModel <- 1
+ }
- OutputsModelDown <- FUN_MOD(InputsModel = InputsModel, RunOptions = RunOptions, Param = Param)
-
+ OutputsModel <- FUN_MOD(InputsModel = InputsModel, RunOptions = RunOptions,
+ Param = Param[iFirstParamRunOffModel:length(Param)])
+
+ if (inherits(InputsModel, "SD")) {
+ InputsModel$OutputsModel <- OutputsModel
+ OutputsModel <- RunModel_LAG(InputsModel, RunOptions, Param[1])
}
- return(OutputsModelDown)
+ return(OutputsModel)
}
\ No newline at end of file
diff --git a/R/RunModel_LAG.R b/R/RunModel_LAG.R
new file mode 100644
index 00000000..3f06a19d
--- /dev/null
+++ b/R/RunModel_LAG.R
@@ -0,0 +1,43 @@
+RunModel_LAG <- function(InputsModel,RunOptions,Param) {
+
+ OutputsModel <- InputsModel$OutputsModel
+ OutputsModel$QsimDown <- OutputsModel$Qsim
+
+ if (inherits(InputsModel, "daily")) {
+ TimeStep <- 60 * 60 * 24
+ }
+ if (inherits(InputsModel, "hourly")) {
+ TimeStep <- 60 * 60
+ }
+
+ # propagation time from upstream meshes to outlet
+ PT <- InputsModel$LengthHydro / Param[length(Param)] / TimeStep
+ HUTRANS <- rbind(1 - (PT - floor(PT)), PT - floor(PT))
+
+ NbUpBasins <- length(InputsModel$LengthHydro)
+ LengthTs <- length(OutputsModel$QsimDown)
+ OutputsModel$Qsim <- OutputsModel$QsimDown * InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1E3
+
+ for (upstream_basin in seq_len(NbUpBasins)) {
+ Qupstream <- InputsModel$Qupstream[RunOptions$IndPeriod_Run, upstream_basin]
+ if(!is.na(InputsModel$BasinAreas[upstream_basin])) {
+ # Upstream flow with area needs to be converted to m3 by time step
+ Qupstream <- Qupstream * InputsModel$BasinAreas[upstream_basin] * 1E3
+ }
+ OutputsModel$Qsim <- OutputsModel$Qsim +
+ c(rep(0, floor(PT[upstream_basin])),
+ Qupstream[1:(LengthTs - floor(PT[upstream_basin]))]) *
+ HUTRANS[1, upstream_basin] +
+ c(rep(0, floor(PT[upstream_basin] + 1)),
+ Qupstream[1:(LengthTs - floor(PT[upstream_basin]) - 1)]) *
+ HUTRANS[2, upstream_basin]
+ }
+ # Warning for negative flows
+ if(any(OutputsModel$Qsim < 0)) {
+ warning(length(which(OutputsModel$Qsim < 0)), " time steps with negative flow, set to zero.")
+ OutputsModel$Qsim[OutputsModel$Qsim < 0] <- 0
+ }
+ # Convert back Qsim to mm
+ OutputsModel$Qsim <- OutputsModel$Qsim / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1E3
+ return(OutputsModel)
+}
\ No newline at end of file
diff --git a/tests/testthat/test-RunModel_LAG.R b/tests/testthat/test-RunModel_LAG.R
index a1fb28e6..f46338dc 100644
--- a/tests/testthat/test-RunModel_LAG.R
+++ b/tests/testthat/test-RunModel_LAG.R
@@ -67,7 +67,7 @@ test_that("Upstream basin with nil area should return same Qdown as GR4J alone",
OutputsSD <-
RunModel(InputsModel,
RunOptions,
- Param = c(Param, 1),
+ Param = c(1, Param),
FUN_MOD = RunModel_GR4J)
expect_equal(OutputsGR4JOnly$Qsim, OutputsSD$Qsim)
})
@@ -80,13 +80,13 @@ test_that(
OutputsSD <-
RunModel(InputsModel,
RunOptions,
- Param = c(Param, 1),
+ Param = c(1, Param),
FUN_MOD = RunModel_GR4J)
expect_equal(OutputsSD$Qsim, Qupstream[Ind_Run])
}
)
-ParamSD = c(Param, InputsModel$LengthHydro / (24 * 60 * 60)) # Speed corresponding to one time step delay
+ParamSD = c(InputsModel$LengthHydro / (24 * 60 * 60), Param) # Speed corresponding to one time step delay
QlsGR4Only <-
OutputsGR4JOnly$Qsim * InputsModel$BasinAreas[2] * 1E6 / 86400
@@ -103,7 +103,7 @@ test_that("1 input with lag of 1 time step delay out gives an output delayed of
test_that("1 input with lag of 0.5 time step delay out gives an output delayed of 0.5 time step", {
OutputsSD <-
- RunModel(InputsModel, RunOptions, Param = c(Param, InputsModel$LengthHydro / (12 * 3600)), FUN_MOD = RunModel_GR4J)
+ RunModel(InputsModel, RunOptions, Param = c(InputsModel$LengthHydro / (12 * 3600), Param), FUN_MOD = RunModel_GR4J)
QlsSdSim <-
OutputsSD$Qsim * sum(InputsModel$BasinAreas) * 1E6 / 86400
QlsUpstLagObs <-
@@ -135,8 +135,8 @@ test_that("Params from calibration with simulated data should be similar to init
CalibOptions = CalibOptions,
FUN_MOD = RunModel_GR4J
)
- expect_equal(OutputsCalib$ParamFinalR[1:4] / ParamSD[1:4], rep(1, 4), tolerance = 1E-2)
- expect_equal(OutputsCalib$ParamFinalR[5], ParamSD[5], tolerance = 2E-3)
+ expect_equal(OutputsCalib$ParamFinalR[2:5] / ParamSD[2:5], rep(1, 4), tolerance = 1E-2)
+ expect_equal(OutputsCalib$ParamFinalR[1], ParamSD[1], tolerance = 2E-3)
})
test_that("1 no area input with lag of 1 time step delay out gives an output delayed of one time step converted to mm", {
--
GitLab