airGR issueshttps://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues2021-03-02T17:57:29+01:00https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/104RunModel_Lag: `StateEnd` is wrong if upstream flow is in mm / time step2021-03-02T17:57:29+01:00Dorchies DavidRunModel_Lag: `StateEnd` is wrong if upstream flow is in mm / time stepIn the lag model upstream flows related to a sub-basin are converted from mm / time step to m3 / time step. Initial states and upstream flows are merged in order to have the memory of the upstream flows before the simulation period but no conversion is applied to initial states.
So the initial states should either be converted if they are stored in mm/ts or stored in m3/ts.
There are pro and cons for the 2 solutions:
1. store states in mm / ts and converted in simulation: `Qupstream` is already in mm/ts so it seems consistent to store the corresponding initial state (which is `Qupstream` before the simulation period) in the same unity.
2. store states in m3/ts and use it directly in simulation which can avoid a computation burden due to a conversion at each run.
Concerning the second solution, maybe it could be a good idea to also convert `Qupstream` in m3/ts in `CreateInputsModel` in order to avoid an "on-the-fly" conversion at each run of `RunModel_Lag`. So, all upstream flows (States + Qupstream) would be in m3/ts and no computation burden would occur due to conversion.In the lag model upstream flows related to a sub-basin are converted from mm / time step to m3 / time step. Initial states and upstream flows are merged in order to have the memory of the upstream flows before the simulation period but no conversion is applied to initial states.
So the initial states should either be converted if they are stored in mm/ts or stored in m3/ts.
There are pro and cons for the 2 solutions:
1. store states in mm / ts and converted in simulation: `Qupstream` is already in mm/ts so it seems consistent to store the corresponding initial state (which is `Qupstream` before the simulation period) in the same unity.
2. store states in m3/ts and use it directly in simulation which can avoid a computation burden due to a conversion at each run.
Concerning the second solution, maybe it could be a good idea to also convert `Qupstream` in m3/ts in `CreateInputsModel` in order to avoid an "on-the-fly" conversion at each run of `RunModel_Lag`. So, all upstream flows (States + Qupstream) would be in m3/ts and no computation burden would occur due to conversion.https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/103RunModel_Lag: `StateEnd` is wrong in case of more than 1 upstream basin2021-03-02T08:41:05+01:00Dorchies DavidRunModel_Lag: `StateEnd` is wrong in case of more than 1 upstream basinThe code refers to the variable `Qupstream` which is the flow of the last upstream flow processed in the loop of lag calculation:
```r
OutputsModel$StateEnd$SD <- lapply(seq(NbUpBasins), function(x) {
Qupstream[(LengthTs - floor(PT[x])):LengthTs]
})
```
It should refer to `InputsModel$Qupstream` which is the matrix of all upstream flows.The code refers to the variable `Qupstream` which is the flow of the last upstream flow processed in the loop of lag calculation:
```r
OutputsModel$StateEnd$SD <- lapply(seq(NbUpBasins), function(x) {
Qupstream[(LengthTs - floor(PT[x])):LengthTs]
})
```
It should refer to `InputsModel$Qupstream` which is the matrix of all upstream flows.v1.7https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/102RunModel_Lag returns 2 values for a one time step run2021-03-02T17:57:29+01:00Dorchies DavidRunModel_Lag returns 2 values for a one time step runThe example below tries to run a daily SD model for the time step `1990-01-01` (adapted from `tests\testthat\test-RunModel_Lag.R`):
```r
library(airGR)
data(L0123001)
BasinAreas <- c(BasinInfo$BasinArea, BasinInfo$BasinArea)
# 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))
InputsModel <- CreateInputsModel(
FUN_MOD = RunModel_GR4J,
DatesR = BasinObs$DatesR,
Precip = BasinObs$P,
PotEvap = BasinObs$E,
Qupstream = matrix(Qupstream, ncol = 1),
LengthHydro = 1000,
BasinAreas = BasinAreas
)
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1990-01-01"),
which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1990-01-01"))
RunOptions <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel,
IndPeriod_Run = Ind_Run))
Param <- c(1., 257.237556, 1.012237, 88.234673, 2.207958)
OutputsModel <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions, Param = Param, FUN_MOD = RunModel_GR4J)
OutputsModel$Qsim
```
The output of this code is:
```r
> OutputsModel$Qsim
[1] 1.709998 1.715785
```
The error comes from the lines in `R\RunModel_Lag.R`:
```r
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]
```
Especially from the expression `(LengthTs - floor(PT[upstream_basin]) - 1)` which is equal to zero when `LengthTs = 1`. Initial states and Qupstream should be merge only if we need to address them.
I tried locally to solve the problem with this solution:
```r
Qupstream1 <- c(IniStates[[upstream_basin]][-length(IniStates[[upstream_basin]])], Qupstream[1:(LengthTs - floor(PT[upstream_basin]))])
Qupstream2 <- IniStates[[upstream_basin]]
if(LengthTs - floor(PT[upstream_basin]) - 1 > 0) Qupstream2 <- c(Qupstream2, Qupstream[1:(LengthTs - floor(PT[upstream_basin]) - 1)])
OutputsModel$Qsim <- OutputsModel$Qsim +
Qupstream1 * HUTRANS[1, upstream_basin] +
Qupstream2 * HUTRANS[2, upstream_basin]
```The example below tries to run a daily SD model for the time step `1990-01-01` (adapted from `tests\testthat\test-RunModel_Lag.R`):
```r
library(airGR)
data(L0123001)
BasinAreas <- c(BasinInfo$BasinArea, BasinInfo$BasinArea)
# 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))
InputsModel <- CreateInputsModel(
FUN_MOD = RunModel_GR4J,
DatesR = BasinObs$DatesR,
Precip = BasinObs$P,
PotEvap = BasinObs$E,
Qupstream = matrix(Qupstream, ncol = 1),
LengthHydro = 1000,
BasinAreas = BasinAreas
)
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1990-01-01"),
which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1990-01-01"))
RunOptions <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel,
IndPeriod_Run = Ind_Run))
Param <- c(1., 257.237556, 1.012237, 88.234673, 2.207958)
OutputsModel <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions, Param = Param, FUN_MOD = RunModel_GR4J)
OutputsModel$Qsim
```
The output of this code is:
```r
> OutputsModel$Qsim
[1] 1.709998 1.715785
```
The error comes from the lines in `R\RunModel_Lag.R`:
```r
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]
```
Especially from the expression `(LengthTs - floor(PT[upstream_basin]) - 1)` which is equal to zero when `LengthTs = 1`. Initial states and Qupstream should be merge only if we need to address them.
I tried locally to solve the problem with this solution:
```r
Qupstream1 <- c(IniStates[[upstream_basin]][-length(IniStates[[upstream_basin]])], Qupstream[1:(LengthTs - floor(PT[upstream_basin]))])
Qupstream2 <- IniStates[[upstream_basin]]
if(LengthTs - floor(PT[upstream_basin]) - 1 > 0) Qupstream2 <- c(Qupstream2, Qupstream[1:(LengthTs - floor(PT[upstream_basin]) - 1)])
OutputsModel$Qsim <- OutputsModel$Qsim +
Qupstream1 * HUTRANS[1, upstream_basin] +
Qupstream2 * HUTRANS[2, upstream_basin]
```v1.7