diff --git a/.regressionignore b/.regressionignore
index ab77dbfc6a0c0b933bdea02ba87b2836b4462036..def03cd47f4bfe1a3e6ef42a566d2e0d1b415d58 100644
--- a/.regressionignore
+++ b/.regressionignore
@@ -4,10 +4,7 @@
# ignored variable : [Topic]<SPACE>[Variable].
# Example for ignoring OutputsModel variable produced by example("RunModel_GR2M"): RunModel_GR2M OutputsModel
-RunModel_Lag LenghtHydro
-RunModel_Lag LengthHydro
RunModel_Lag Lag
+RunModel_Lag LenghtHydro
RunModel_Lag InputsModel
-RunModel_Lag OutputsModelDown
RunModel_Lag OutputsModel
-
diff --git a/R/CreateInputsModel.R b/R/CreateInputsModel.R
index d09031ca95dc9833571fe46cb383673f10b816a8..3bd0f14752105a78f3010a3d8cd744d646f6fa2e 100644
--- a/R/CreateInputsModel.R
+++ b/R/CreateInputsModel.R
@@ -210,7 +210,10 @@ CreateInputsModel <- function(FUN_MOD,
stop("'Qupstream' must have same number of rows as 'DatesR' length")
}
if(any(is.na(Qupstream))) {
- stop("'Qupstream' cannot contain any NA value")
+ warning("'Qupstream' contains NA values: model outputs will contain NAs")
+ }
+ if(any(LengthHydro > 1000)) {
+ warning("The unit of 'LengthHydro' has changed from m to km in v1.7 of airGR: values superior to 1000 km seem unrealistic")
}
}
diff --git a/R/RunModel_Lag.R b/R/RunModel_Lag.R
index 87f551f73f94e3ef91b93798b78ae35f2db5c0b6..f1ed6a2cd40bdd147d852b25b62a3a52f53eb01a 100644
--- a/R/RunModel_Lag.R
+++ b/R/RunModel_Lag.R
@@ -18,19 +18,13 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param) {
stop(paste("'Param' must be a vector of length", NParam, "and contain no NA"))
}
if (is.null(InputsModel$OutputsModel)) {
- stop(
- "'InputsModel' should contain an 'OutputsModel' key containing the output of the runoff of the downstream subcatchment"
- )
+ stop("'InputsModel' should contain an 'OutputsModel' key containing the output of the runoff of the downstream subcatchment")
}
if (is.null(InputsModel$OutputsModel$Qsim)) {
- stop(
- "'InputsModel$OutputsModel' should contain a key 'Qsim' containing the output of the runoff of the downstream subcatchment"
- )
+ stop("'InputsModel$OutputsModel' should contain a key 'Qsim' containing the output of the runoff of the downstream subcatchment")
}
if (sum(!is.na(InputsModel$OutputsModel$Qsim)) != length(RunOptions$IndPeriod_Run)) {
- stop(
- "'InputsModel$OutputsModel$Qim' should have the same lenght as 'RunOptions$IndPeriod_Run' and contain no NA"
- )
+ stop("'InputsModel$OutputsModel$Qim' should have the same lenght as 'RunOptions$IndPeriod_Run' and contain no NA")
}
OutputsModel <- InputsModel$OutputsModel
@@ -45,7 +39,7 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param) {
}
# propagation time from upstream meshes to outlet
- PT <- InputsModel$LengthHydro / Param[1L] / TimeStep
+ PT <- InputsModel$LengthHydro * 1e3 / Param[1L] / TimeStep
HUTRANS <- rbind(1 - (PT - floor(PT)), PT - floor(PT))
NbUpBasins <- length(InputsModel$LengthHydro)
@@ -75,33 +69,41 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param) {
rep(0, floor(PT[x] + 1))
})
}
+ # message("Initstates: ", paste(IniStates, collapse = ", "))
for (upstream_basin in seq_len(NbUpBasins)) {
- Qupstream <- InputsModel$Qupstream[RunOptions$IndPeriod_Run, upstream_basin]
+ Qupstream <- c(IniStates[[upstream_basin]],
+ 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
}
+ # message("Qupstream[", upstream_basin, "]: ", paste(Qupstream, collapse = ", "))
OutputsModel$Qsim <- OutputsModel$Qsim +
- c(IniStates[[upstream_basin]][-length(IniStates[[upstream_basin]])],
- Qupstream[1:(LengthTs - floor(PT[upstream_basin]))]) *
- HUTRANS[1, upstream_basin] +
- c(IniStates[[upstream_basin]],
- Qupstream[1:(LengthTs - floor(PT[upstream_basin]) - 1)]) *
- HUTRANS[2, upstream_basin]
+ Qupstream[2:(1 + LengthTs)] * HUTRANS[1, upstream_basin] +
+ Qupstream[1:LengthTs] * 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
+ # Warning for negative flows or NAs only in extended outputs
+ if(length(RunOptions$Outputs_Sim) > 2) {
+ if (any(OutputsModel$Qsim[!is.na(OutputsModel$Qsim)] < 0)) {
+ warning(length(which(OutputsModel$Qsim < 0)), " time steps with negative flow, set to zero.")
+ OutputsModel$Qsim[OutputsModel$Qsim < 0] <- 0
+ }
+ # Warning for NAs
+ if (any(is.na(OutputsModel$Qsim))) {
+ warning(length(which(is.na(OutputsModel$Qsim))), " time steps with NA values")
+ }
}
# Convert back Qsim to mm
OutputsModel$Qsim <- OutputsModel$Qsim / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
+ # message("Qsim: ", paste(OutputsModel$Qsim, collapse = ", "))
if ("StateEnd" %in% RunOptions$Outputs_Sim) {
OutputsModel$StateEnd$SD <- lapply(seq(NbUpBasins), function(x) {
- Qupstream[(LengthTs - floor(PT[x])):LengthTs]
+ lastTS <- RunOptions$IndPeriod_Run[length(RunOptions$IndPeriod_Run)]
+ InputsModel$Qupstream[(lastTS - floor(PT[x])):lastTS, x]
})
+ # message("StateEnd: ", paste(OutputsModel$StateEnd$SD, collapse = ", "))
}
return(OutputsModel)
diff --git a/man/CreateInputsModel.Rd b/man/CreateInputsModel.Rd
index 958f9e7eec4e2216912acb52c567a6636c443232..365955eb5d5af4700695578bcaae7e61ebfba4ae 100644
--- a/man/CreateInputsModel.Rd
+++ b/man/CreateInputsModel.Rd
@@ -52,7 +52,7 @@ CreateInputsModel(FUN_MOD, DatesR, Precip, PrecipScale = TRUE, PotEvap = NULL,
\item{Qupstream}{(optional) [numerical matrix] time series of upstream flows (catchment average) [mm/time step or m3/time step, see details], required to create the SD model inputs . See details}
-\item{LengthHydro}{(optional) [numeric] real giving the distance between the downstream outlet and each upstream inlet of the sub-catchment [m], required to create the SD model inputs . See details}
+\item{LengthHydro}{(optional) [numeric] real giving the distance between the downstream outlet and each upstream inlet of the sub-catchment [km], required to create the SD model inputs . See details}
\item{BasinAreas}{(optional) [numeric] real giving the area of each upstream sub-catchment [km2] and the area of the downstream sub-catchment in the last item, required to create the SD model inputs . See details}
@@ -83,7 +83,6 @@ Please note that if CemaNeige is used, and \code{ZInputs} is different than \cod
Users wanting to use a semi-distributed (SD) lag model should provide valid \code{Qupstream}, \code{LengthHydro}, and \code{BasinAreas} parameters. Each upstream sub-catchment is described by an upstream flow time series (one column in \code{Qupstream} matrix), a distance between the downstream outlet and the upstream inlet (one item in \code{LengthHydro}) and an area (one item in \code{BasinAreas}).
The order of the columns or of the items should be consistent for all these parameters. \code{BasinAreas} should contain one extra information (stored in the last item) which is the area of the downstream sub-catchment.
Upstream flows that are not related to a sub-catchment such as a release or withdraw spot are identified by an area equal to \code{NA} and an upstream flow expressed in m3/time step instead of mm/time step.
-Please note that the use of SD lag model require to use the \code{\link{RunModel}} function instead of \code{\link{RunModel_GR4J}} or the other \code{RunModel_*} functions.
}
\examples{
diff --git a/man/RunModel_Lag.Rd b/man/RunModel_Lag.Rd
index 1473ce9776390ee0f217d62d58a76a970ef069d4..b37487679c2734efdddda7f066a4ad9133e07cff 100644
--- a/man/RunModel_Lag.Rd
+++ b/man/RunModel_Lag.Rd
@@ -25,7 +25,7 @@ RunModel_Lag(InputsModel, RunOptions, Param)
\item{Param}{[numeric] vector of 1 parameter
\tabular{ll}{
- Lag \tab Mean flow velocity [m/s]
+ Velocity \tab Mean flow velocity [m/s]
}}
}
@@ -38,51 +38,62 @@ The list value contains an extra item named \code{QsimDown} which is a copy of \
\examples{
-library(airGR)
+#####################################################################
+## Simulation of a reservoir with a purpose of low-flow mitigation ##
+#####################################################################
+
+## ---- preparation of the InputsModel object
-## loading catchment data
+## loading package and catchment data
+library(airGR)
data(L0123001)
-## Simulating a reservoir
-# Withdrawing 1 m3/s with an instream flow of 1 m3/s
-Qupstream <- matrix(- unlist(lapply(BasinObs$Qls / 1000 - 1, function(x) {
- min(1, max(0,x, na.rm = TRUE))
-})), ncol = 1)
-# Except between July and November when releasing 3 m3/s
-month <- as.numeric(format(BasinObs$DatesR,"\%m"))
+## the reservoir withdraws 1 m3/s when it's possible considering the flow observed in the basin
+Qupstream <- matrix(-sapply(BasinObs$Qls / 1000 - 1, function(x) {
+ min(1, max(0, x, na.rm = TRUE))
+}), ncol = 1)
+
+## except between July and September when the reservoir releases 3 m3/s for low-flow mitigation
+month <- as.numeric(format(BasinObs$DatesR, "\%m"))
Qupstream[month >= 7 & month <= 9] <- 3
-# Conversion in m3/day
-Qupstream <- Qupstream * 86400
+Qupstream <- Qupstream * 86400 ## Conversion in m3/day
-## The reservoir is not an upstream subcachment: its areas is NA
+## the reservoir is not an upstream subcachment: its areas is NA
BasinAreas <- c(NA, BasinInfo$BasinArea)
-## Delay time between the reservoir and the catchment outlet is 2 days and the distance is 150 km
-LengthHydro <- 150000
-Lag <- LengthHydro / 2 / 86400 * 1000 # Conversion km/day -> m/s
+## delay time between the reservoir and the catchment outlet is 2 days and the distance is 150 km
+LengthHydro <- 150
-## preparation of the InputsModel object
InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
Precip = BasinObs$P, PotEvap = BasinObs$E,
Qupstream = Qupstream, LengthHydro = LengthHydro,
BasinAreas = BasinAreas)
+
+## ---- simulation of the basin with the reservoir influence
+
## 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
+## creation of the RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel, IndPeriod_Run = Ind_Run)
-## simulation of downstream subcatchment
+## simulation of the runoff of the catchment with a GR4J model
Param <- c(X1 = 257.238, X2 = 1.012, X3 = 88.235, X4 = 2.208)
OutputsModelDown <- RunModel_GR4J(InputsModel = InputsModel,
RunOptions = RunOptions, Param = Param)
+## with a delay of 2 days for 150 km, the flow velocity is 75 km per day
+Velocity <- (LengthHydro * 1e3 / 2) / (24 * 60 * 60) ## Conversion km/day -> m/s
+
+## add the output of the precipitation-runoff model in the InputsModel object
InputsModel$OutputsModel <- OutputsModelDown
+
+## run the lag model which routes precipitation-runoff model and upstream flows
OutputsModel <- RunModel_Lag(InputsModel = InputsModel,
- RunOptions = RunOptions, Param = Lag)
+ RunOptions = RunOptions, Param = Velocity)
## results preview of comparison between naturalised (observed) and influenced flow (simulated)
plot(OutputsModel, Qobs = OutputsModel$QsimDown)
diff --git a/tests/testthat/test-RunModel_Lag.R b/tests/testthat/test-RunModel_Lag.R
index bfef15d2833076d4dc30005615c46d9ed00eb2fe..0432fa091351418eaa4bfeee93a906dd805183f6 100644
--- a/tests/testthat/test-RunModel_Lag.R
+++ b/tests/testthat/test-RunModel_Lag.R
@@ -19,21 +19,6 @@ 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(
- FUN_MOD = RunModel_GR4J,
- DatesR = BasinObs$DatesR,
- Precip = BasinObs$P,
- PotEvap = BasinObs$E,
- Qupstream = matrix(BasinObs$Qmm, ncol = 1),
- LengthHydro = 1,
- BasinAreas = BasinAreas
- ),
- regexp = "'Qupstream' cannot contain any NA value"
- )
-})
-
# Qupstream = sinusoid synchronised on hydrological year from 0 mm to mean value of Qobs
Qupstream <- floor((sin((seq_along(BasinObs$Qmm)/365*2*3.14))+1) * mean(BasinObs$Qmm, na.rm = TRUE))
@@ -43,7 +28,7 @@ InputsModel <- CreateInputsModel(
Precip = BasinObs$P,
PotEvap = BasinObs$E,
Qupstream = matrix(Qupstream, ncol = 1),
- LengthHydro = 1000,
+ LengthHydro = 1,
BasinAreas = BasinAreas
)
@@ -85,7 +70,7 @@ test_that("'InputsModel$OutputsModel$Qim' should have the same lenght as 'RunOpt
)
})
-test_that("'InputsModel$OutputsModel$Qim' should contain no NA'", {
+test_that("'InputsModel$OutputsModel$Qsim' should contain no NA'", {
InputsModel$OutputsModel <- OutputsGR4JOnly
InputsModel$OutputsModel$Qsim[10L] <- NA
expect_error(
@@ -94,6 +79,37 @@ test_that("'InputsModel$OutputsModel$Qim' should contain no NA'", {
)
})
+test_that("'Qupstream' contain NA values", {
+ expect_warning(
+ InputsModel <- CreateInputsModel(
+ FUN_MOD = RunModel_GR4J,
+ DatesR = BasinObs$DatesR,
+ Precip = BasinObs$P,
+ PotEvap = BasinObs$E,
+ Qupstream = matrix(BasinObs$Qmm, ncol = 1),
+ LengthHydro = 1,
+ BasinAreas = BasinAreas
+ ),
+ regexp = "'Qupstream' contains NA values: model outputs will contain NAs"
+ )
+
+ RunOptions <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
+ InputsModel = InputsModel,
+ IndPeriod_Run = Ind_Run))
+ InputsModel$OutputsModel <- OutputsGR4JOnly
+ # Warning with RunModel
+ expect_warning(
+ RunModel_Lag(InputsModel = InputsModel, RunOptions = RunOptions, Param = 1),
+ regexp = "time steps with NA values"
+ )
+ # No warning during calibration
+ RunOptions$Outputs_Sim <- RunOptions$Outputs_Cal
+ expect_warning(
+ RunModel_Lag(InputsModel = InputsModel, RunOptions = RunOptions, Param = 1),
+ regexp = NA
+ )
+})
+
test_that("Upstream basin with nil area should return same Qdown as GR4J alone", {
UpstBasinArea <- InputsModel$BasinAreas[1L]
InputsModel$BasinAreas[1L] <- 0
@@ -114,7 +130,7 @@ test_that("Downstream basin with nil area and nul upstream length should return
expect_equal(OutputsSD$Qsim, Qupstream[Ind_Run])
})
-ParamSD <- c(InputsModel$LengthHydro / (24 * 60 * 60), Param) # Speed corresponding to one time step delay
+ParamSD <- c(InputsModel$LengthHydro * 1e3 / (24 * 60 * 60), Param) # Speed corresponding to one time step delay
QlsGR4Only <- OutputsGR4JOnly$Qsim * InputsModel$BasinAreas[2L] * 1e6 / 86400
@@ -127,7 +143,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(InputsModel$LengthHydro / (12 * 3600), Param),
+ Param = c(InputsModel$LengthHydro * 1e3 / (12 * 3600), Param),
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[1L] * 1e6 / 86400
@@ -180,7 +196,7 @@ test_that("1 no area input with lag of 1 time step delay out gives an output del
IM <- InputsModel
IM$BasinAreas <- rep(BasinInfo$BasinArea, 3)
IM$Qupstream <- cbind(IM$Qupstream, IM$Qupstream)
-IM$LengthHydro <- c(1000, 1500)
+IM$LengthHydro <- c(1, 1.5)
PSDini <- ParamSD
PSDini[1] <- PSDini[1] / 2 # 2 time step delay
diff --git a/vignettes/V05_sd_model.Rmd b/vignettes/V05_sd_model.Rmd
index 383d8a93cc9e35ade7d23c98a962c8991aa687db..632819c13cdd8c32d751f23f4933467a030de2d8 100644
--- a/vignettes/V05_sd_model.Rmd
+++ b/vignettes/V05_sd_model.Rmd
@@ -12,7 +12,6 @@ vignette: >
```{r, include=FALSE, fig.keep='none', results='hide'}
library(airGR)
options(digits = 3)
-library(imputeTS)
```
# Introduction
@@ -56,7 +55,9 @@ For the observed flow at the downstream outlet, we generate it with the assumpti
```{r}
QObsDown <- (BasinObs$Qmm + c(0, 0, BasinObs$Qmm[1:(length(BasinObs$Qmm)-2)])) / 2
+options(digits = 5)
summary(cbind(QObsUp = BasinObs$Qmm, QObsDown))
+options(digits = 3)
```
# Calibration of the upstream subcatchment
@@ -91,20 +92,14 @@ OutputsModelUp <- RunModel_GR4J(InputsModel = InputsModelUp, RunOptions = RunOpt
# Calibration of the downstream subcatchment with upstream flow observations
-Observed flow data contain `NA` values and a complete time series is mandatory for running the Lag model. We propose to complete the observed upstream flow with linear interpolation:
-
-```{r}
-QObsUp <- imputeTS::na_interpolation(BasinObs$Qmm)
-```
-
we need to create the `InputsModel` object completed with upstream information:
```{r}
InputsModelDown1 <- CreateInputsModel(
FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
Precip = BasinObs$P, PotEvap = BasinObs$E,
- Qupstream = matrix(QObsUp, ncol = 1), # upstream observed flow
- LengthHydro = 1e2 * 1e3, # distance between upstream catchment outlet & the downstream one [m]
+ Qupstream = matrix(BasinObs$Qmm, ncol = 1), # upstream observed flow
+ LengthHydro = 100, # distance between upstream catchment outlet & the downstream one [km]
BasinAreas = c(180, 180) # upstream and downstream areas [km²]
)
```
@@ -168,27 +163,27 @@ ParamDown2 <- OutputsCalibDown2$ParamFinalR
# Discussion
-## Identification of Lag parameter
+## Identification of Velocity parameter
-The theoretical Lag parameter should be equal to:
+The theoretical Velocity parameter should be equal to:
```{r}
-Lag <- InputsModelDown1$LengthHydro / (2 * 86400)
-paste(format(Lag), "m/s")
+Velocity <- InputsModelDown1$LengthHydro * 1e3 / (2 * 86400)
+paste(format(Velocity), "m/s")
```
Both calibrations overestimate this parameter:
```{r}
-mLag <- matrix(c(Lag,
- OutputsCalibDown1$ParamFinalR[1],
- OutputsCalibDown2$ParamFinalR[1]),
- ncol = 1,
- dimnames = list(c("theoretical",
- "calibrated with observed upstream flow",
- "calibrated with simulated upstream flow"),
- c("Lag parameter")))
-knitr::kable(mLag)
+mVelocity <- matrix(c(Velocity,
+ OutputsCalibDown1$ParamFinalR[1],
+ OutputsCalibDown2$ParamFinalR[1]),
+ ncol = 1,
+ dimnames = list(c("theoretical",
+ "calibrated with observed upstream flow",
+ "calibrated with simulated upstream flow"),
+ c("Velocity parameter")))
+knitr::kable(mVelocity)
```
## Value of the performance criteria with theoretical calibration
@@ -196,7 +191,7 @@ knitr::kable(mLag)
Theoretically, the parameters of the downstream GR4J model should be the same as the upstream one and we know the lag time. So this set of parameter should give a better performance criteria:
```{r}
-ParamDownTheo <- c(Lag, OutputsCalibUp$ParamFinalR)
+ParamDownTheo <- c(Velocity, OutputsCalibUp$ParamFinalR)
OutputsModelDownTheo <- RunModel(InputsModel = InputsModelDown2,
RunOptions = RunOptionsDown,
Param = ParamDownTheo,
@@ -219,7 +214,7 @@ comp <- cbind(comp, c(OutputsCalibUp$CritFinal,
CritDown1$CritValue,
OutputsCalibDown2$CritFinal,
CritDownTheo$CritValue))
-colnames(comp) <- c("Lag", paste0("X", 1:4), "NSE")
+colnames(comp) <- c("Velocity", paste0("X", 1:4), "NSE")
rownames(comp) <- c("Calibration of the upstream subcatchment",
"Calibration 1 with observed upstream flow",
"Validation 1 with simulated upstream flow",