diff --git a/DESCRIPTION b/DESCRIPTION
index a517127e0f90c09bc051d3cd465221930d762036..44ca4908f104de6441c1fe698e680cdc700b98d0 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
Package: airGR
Type: Package
Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
-Version: 1.6.2.2
-Date: 2020-06-05
+Version: 1.6.2.3
+Date: 2020-07-16
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"),
diff --git a/NEWS.md b/NEWS.md
index 173a46a9284cfde05987613917c5c4da7443a314..d3b32a32a338b7d1308f5cbe16322c669b80b097 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,14 @@
## Release History of the airGR Package
+### 1.6.2.2 Release Notes (2020-07-15)
+
+#### Bug fixes
+
+- Major error in lag model implementation
+
+____________________________________________________________________________________
+
+
### 1.6.2.1 Release Notes (2020-06-05)
#### New features
diff --git a/R/RunModel.R b/R/RunModel.R
index dabe11d2aae2dda1651dede10b31eb2037b48a59..9ca647e892c6cf21ea7fad4d77e99ddef7f240c6 100644
--- a/R/RunModel.R
+++ b/R/RunModel.R
@@ -30,10 +30,10 @@ RunModel <- function (InputsModel, RunOptions, Param, FUN_MOD) {
}
OutputsModelDown$Qsim <- OutputsModelDown$Qsim +
c(rep(0, floor(PT[upstream_basin])),
- Qupstream[(1 + floor(PT[upstream_basin])):LengthTs]) *
+ Qupstream[1:(LengthTs - floor(PT[upstream_basin]))]) *
HUTRANS[1, upstream_basin] +
c(rep(0, floor(PT[upstream_basin] + 1)),
- Qupstream[(2 + floor(PT[upstream_basin])):LengthTs]) *
+ Qupstream[1:(LengthTs - floor(PT[upstream_basin]) - 1)]) *
HUTRANS[2, upstream_basin]
}
# Warning for negative flows
diff --git a/tests/testthat/test-RunModel_LAG.R b/tests/testthat/test-RunModel_LAG.R
index 42a18ceb240f9c51f0be93a844fdee3636efdd83..a1fb28e624c25652f44e7ce76da3aa84b662490e 100644
--- a/tests/testthat/test-RunModel_LAG.R
+++ b/tests/testthat/test-RunModel_LAG.R
@@ -2,7 +2,6 @@ context("RunModel_LAG")
data(L0123001)
-
test_that("'BasinAreas' must have one more element than 'LengthHydro'", {
expect_error(
InputsModel <- CreateInputsModel(
@@ -18,6 +17,8 @@ test_that("'BasinAreas' must have one more element than 'LengthHydro'", {
)
})
+BasinAreas <- c(BasinInfo$BasinArea, BasinInfo$BasinArea)
+
test_that("'Qupstream' cannot contain any NA value", {
expect_error(
InputsModel <- CreateInputsModel(
@@ -27,14 +28,14 @@ test_that("'Qupstream' cannot contain any NA value", {
PotEvap = BasinObs$E,
Qupstream = matrix(BasinObs$Qmm, ncol = 1),
LengthHydro = 1,
- BasinAreas = c(1, 2)
+ BasinAreas = BasinAreas
),
regexp = "'Qupstream' cannot contain any NA value"
)
})
-Qupstream = BasinObs$Qmm
-Qupstream[is.na(Qupstream)] = mean(Qupstream, na.rm = TRUE)
+# Qupstream = sinusoid synchronised on hydrological year from 0 mm to mean value of Qobs
+Qupstream = floor((sin((1:length(BasinObs$Qmm)/365*2*3.14))+1)*mean(BasinObs$Qmm, na.rm = T))
InputsModel <- CreateInputsModel(
FUN_MOD = RunModel_GR4J,
@@ -43,7 +44,7 @@ InputsModel <- CreateInputsModel(
PotEvap = BasinObs$E,
Qupstream = matrix(Qupstream, ncol = 1),
LengthHydro = 1000,
- BasinAreas = c(BasinInfo$BasinArea * 2, BasinInfo$BasinArea)
+ BasinAreas = BasinAreas
)
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1990-01-01"),
@@ -60,7 +61,6 @@ OutputsGR4JOnly <-
RunOptions = RunOptions,
Param = Param)
-
test_that("Upstream basin with nil area should return same Qdown as GR4J alone", {
UpstBasinArea = InputsModel$BasinAreas[1]
InputsModel$BasinAreas[1] <- 0
@@ -88,25 +88,40 @@ test_that(
ParamSD = c(Param, InputsModel$LengthHydro / (24 * 60 * 60)) # Speed corresponding to one time step delay
+QlsGR4Only <-
+ OutputsGR4JOnly$Qsim * InputsModel$BasinAreas[2] * 1E6 / 86400
+
test_that("1 input with lag of 1 time step delay out gives an output delayed of one time step", {
- QlsGR4Only <-
- OutputsGR4JOnly$Qsim * InputsModel$BasinAreas[2] * 1E6 / 86400
OutputsSD <-
RunModel(InputsModel, RunOptions, Param = ParamSD, FUN_MOD = RunModel_GR4J)
QlsSdSim <-
OutputsSD$Qsim * sum(InputsModel$BasinAreas) * 1E6 / 86400
QlsUpstLagObs <-
- c(0, Qupstream[Ind_Run[1:(length(Ind_Run) - 1)] + 1]) * InputsModel$BasinAreas[1] * 1E6 / 86400
+ c(0, Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]]) * InputsModel$BasinAreas[1] * 1E6 / 86400
expect_equal(QlsSdSim - QlsGR4Only, QlsUpstLagObs)
})
+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)
+ QlsSdSim <-
+ OutputsSD$Qsim * sum(InputsModel$BasinAreas) * 1E6 / 86400
+ QlsUpstLagObs <-
+ (Qupstream[Ind_Run] + c(0, Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]]))/2 * InputsModel$BasinAreas[1] * 1E6 / 86400
+ expect_equal(QlsSdSim - QlsGR4Only, QlsUpstLagObs)
+})
+
+
test_that("Params from calibration with simulated data should be similar to initial params", {
InputsCrit <- CreateInputsCrit(
FUN_CRIT = ErrorCrit_NSE,
InputsModel = InputsModel,
RunOptions = RunOptions,
VarObs = "Q",
- Obs = BasinObs$Qmm[Ind_Run]
+ Obs = (
+ c(0, Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]]) * BasinAreas[1] +
+ BasinObs$Qmm[Ind_Run] * BasinAreas[2]
+ ) / sum(BasinAreas)
)
CalibOptions <- CreateCalibOptions(
FUN_MOD = RunModel_GR4J,
@@ -114,31 +129,32 @@ test_that("Params from calibration with simulated data should be similar to init
IsSD = TRUE
)
OutputsCalib <- Calibration_Michel(
- InputsModel = InputsModel,
- RunOptions = RunOptions,
- InputsCrit = InputsCrit,
- CalibOptions = CalibOptions,
- FUN_MOD = RunModel_GR4J
- )
- expect_equal(OutputsCalib$ParamFinalR, ParamSD, tolerance = 1E-3)
+ InputsModel = InputsModel,
+ RunOptions = RunOptions,
+ InputsCrit = InputsCrit,
+ 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)
})
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", {
Qm3GR4Only <-
- OutputsGR4JOnly$Qsim * InputsModel$BasinAreas[2] * 1E3
+ OutputsGR4JOnly$Qsim * BasinAreas[2] * 1E3
# Specify that upstream flow is not related to an area
- InputsModel$BasinAreas = c(NA, BasinInfo$BasinArea)
- # Convert upstream flow to Liter/day
- InputsModel$Qupstream <- matrix(Qupstream, ncol = 1) * BasinInfo$BasinArea * 2
+ InputsModel$BasinAreas = c(NA, BasinAreas[2])
+ # Convert upstream flow to m3/day
+ InputsModel$Qupstream <- matrix(Qupstream, ncol = 1) * BasinAreas[1] * 1E3
+
OutputsSD <-
RunModel(InputsModel, RunOptions, Param = ParamSD, FUN_MOD = RunModel_GR4J)
-
+
expect_false(any(is.na(OutputsSD$Qsim)))
-
- Qm3SdSim <-
- OutputsSD$Qsim * sum(InputsModel$BasinAreas, na.rm = TRUE) * 1E3
+
+ Qm3SdSim <- OutputsSD$Qsim * sum(InputsModel$BasinAreas, na.rm = TRUE) * 1E3
Qm3UpstLagObs <-
- c(0, InputsModel$Qupstream[Ind_Run[1:(length(Ind_Run) - 1)] + 1])
-
+ c(0, InputsModel$Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]])
+
expect_equal(Qm3SdSim - Qm3GR4Only, Qm3UpstLagObs)
})