test-RunModel_LAG.R 9.16 KB
Newer Older
1
context("RunModel_Lag")
2
3
4

data(L0123001)

5
test_that("'BasinAreas' must have one more element than 'LengthHydro'", {
6
7
8
9
10
11
  expect_error(
    InputsModel <- CreateInputsModel(
      FUN_MOD = RunModel_GR4J,
      DatesR = BasinObs$DatesR,
      Precip = BasinObs$P,
      PotEvap = BasinObs$E,
12
      Qupstream = matrix(BasinObs$Qmm, ncol = 1),
13
14
      LengthHydro = 1,
      BasinAreas = 1
15
    ),
16
    regexp = "'BasinAreas' must have one more element than 'LengthHydro'"
17
  )
18
19
})

20
21
BasinAreas <- c(BasinInfo$BasinArea, BasinInfo$BasinArea)

22
test_that("'Qupstream' cannot contain any NA value", {
23
24
25
26
27
28
  expect_error(
    InputsModel <- CreateInputsModel(
      FUN_MOD = RunModel_GR4J,
      DatesR = BasinObs$DatesR,
      Precip = BasinObs$P,
      PotEvap = BasinObs$E,
29
      Qupstream = matrix(BasinObs$Qmm, ncol = 1),
30
      LengthHydro = 1,
31
      BasinAreas = BasinAreas
32
    ),
33
    regexp = "'Qupstream' cannot contain any NA value"
34
  )
35
36
})

37
# Qupstream = sinusoid synchronised on hydrological year from 0 mm to mean value of Qobs
38
Qupstream <- floor((sin((seq_along(BasinObs$Qmm)/365*2*3.14))+1) * mean(BasinObs$Qmm, na.rm = TRUE))
39
40

InputsModel <- CreateInputsModel(
41
42
43
44
  FUN_MOD = RunModel_GR4J,
  DatesR = BasinObs$DatesR,
  Precip = BasinObs$P,
  PotEvap = BasinObs$E,
45
  Qupstream = matrix(Qupstream, ncol = 1),
46
  LengthHydro = 1000,
47
  BasinAreas = BasinAreas
48
49
)

50
51
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"))
52

53
RunOptions <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
54
55
                                                InputsModel = InputsModel,
                                                IndPeriod_Run = Ind_Run))
56

57
58
test_that("InputsModel parameter should contain an OutputsModel key", {
  expect_error(
59
    RunModel_Lag(InputsModel = InputsModel, RunOptions = RunOptions, Param = 1),
60
61
62
63
    regexp = "'InputsModel' should contain an 'OutputsModel' key"
  )
})

64
Param <- c(257.237556, 1.012237, 88.234673, 2.207958) #  From vignettes/V01_get_started
65

66
67
68
OutputsGR4JOnly <- RunModel_GR4J(InputsModel = InputsModel,
                                 RunOptions = RunOptions,
                                 Param = Param)
69

70
71
72
73
test_that("InputsModel$OutputsModel should contain a Qsim key", {
  InputsModel$OutputsModel <- OutputsGR4JOnly
  InputsModel$OutputsModel$Qsim <- NULL
  expect_error(
74
    RunModel_Lag(InputsModel = InputsModel, RunOptions = RunOptions, Param = 1),
75
76
77
78
79
80
81
82
    regexp = "should contain a key 'Qsim'"
  )
})

test_that("'InputsModel$OutputsModel$Qim' should have the same lenght as 'RunOptions$IndPeriod_Run'", {
  InputsModel$OutputsModel <- OutputsGR4JOnly
  InputsModel$OutputsModel$Qsim <- c(InputsModel$OutputsModel$Qsim, 0)
  expect_error(
83
    RunModel_Lag(InputsModel = InputsModel, RunOptions = RunOptions, Param = 1),
84
85
86
87
88
89
    regexp = "should have the same lenght as"
  )
})

test_that("'InputsModel$OutputsModel$Qim' should contain no NA'", {
  InputsModel$OutputsModel <- OutputsGR4JOnly
90
  InputsModel$OutputsModel$Qsim[10L] <- NA
91
  expect_error(
92
    RunModel_Lag(InputsModel = InputsModel, RunOptions = RunOptions, Param = 1),
93
94
95
96
    regexp = "contain no NA"
  )
})

97
test_that("Upstream basin with nil area should return same Qdown as GR4J alone", {
98
99
100
101
102
103
  UpstBasinArea <- InputsModel$BasinAreas[1L]
  InputsModel$BasinAreas[1L] <- 0
  OutputsSD <- RunModel(InputsModel,
                        RunOptions,
                        Param = c(1, Param),
                        FUN_MOD = RunModel_GR4J)
104
105
106
  expect_equal(OutputsGR4JOnly$Qsim, OutputsSD$Qsim)
})

107
108
109
110
111
112
113
114
115
test_that("Downstream basin with nil area and nul upstream length should return same Qdown as Qupstream alone", {
  InputsModel$LengthHydro <- 0
  InputsModel$BasinAreas <- c(BasinInfo$BasinArea, 0)
  OutputsSD <- RunModel(InputsModel,
                        RunOptions,
                        Param = c(1, Param),
                        FUN_MOD = RunModel_GR4J)
  expect_equal(OutputsSD$Qsim, Qupstream[Ind_Run])
})
116

117
ParamSD <- c(InputsModel$LengthHydro / (24 * 60 * 60), Param) # Speed corresponding to one time step delay
118

119
QlsGR4Only <- OutputsGR4JOnly$Qsim * InputsModel$BasinAreas[2L] * 1e6 / 86400
120

121
test_that("1 input with lag of 1 time step delay out gives an output delayed of one time step", {
122
123
124
  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)]]) * InputsModel$BasinAreas[1L] * 1e6 / 86400
125
126
127
  expect_equal(QlsSdSim - QlsGR4Only, QlsUpstLagObs)
})

128
test_that("1 input with lag of 0.5 time step delay out gives an output delayed of 0.5 time step", {
129
130
131
132
133
  OutputsSD <- RunModel(InputsModel, RunOptions,
                        Param = c(InputsModel$LengthHydro / (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
134
135
136
  expect_equal(QlsSdSim - QlsGR4Only, QlsUpstLagObs)
})

137
138
139
140
141
142
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",
143
144
    Obs = (c(0, Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]]) * BasinAreas[1L] +
             BasinObs$Qmm[Ind_Run] * BasinAreas[2L]) / sum(BasinAreas)
145
146
  )
  CalibOptions <- CreateCalibOptions(
147
148
    FUN_MOD = RunModel_GR4J,
    FUN_CALIB = Calibration_Michel,
149
150
151
    IsSD = TRUE
  )
  OutputsCalib <- Calibration_Michel(
152
153
154
155
156
157
    InputsModel = InputsModel,
    RunOptions = RunOptions,
    InputsCrit = InputsCrit,
    CalibOptions = CalibOptions,
    FUN_MOD = RunModel_GR4J
  )
158
159
  expect_equal(OutputsCalib$ParamFinalR[2:5] / ParamSD[2:5], rep(1, 4), tolerance = 1e-2)
  expect_equal(OutputsCalib$ParamFinalR[1L], ParamSD[1L], tolerance = 2e-3)
160
})
161
162

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", {
163
  Qm3GR4Only <- OutputsGR4JOnly$Qsim * BasinAreas[2L] * 1e3
164
  # Specify that upstream flow is not related to an area
165
  InputsModel$BasinAreas <- c(NA, BasinAreas[2L])
166
  # Convert upstream flow to m3/day
167
  InputsModel$Qupstream <- matrix(Qupstream, ncol = 1) * BasinAreas[1L] * 1e3
168

169
  OutputsSD <- RunModel(InputsModel, RunOptions, Param = ParamSD, FUN_MOD = RunModel_GR4J)
170

171
  expect_false(any(is.na(OutputsSD$Qsim)))
172

173
174
  Qm3SdSim <- OutputsSD$Qsim * sum(InputsModel$BasinAreas, na.rm = TRUE) * 1e3
  Qm3UpstLagObs <- c(0, InputsModel$Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]])
175

176
177
  expect_equal(Qm3SdSim - Qm3GR4Only, Qm3UpstLagObs)
})
178

179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
# *** IniStates tests ***
IM <- InputsModel
IM$BasinAreas <- rep(BasinInfo$BasinArea, 3)
IM$Qupstream <- cbind(IM$Qupstream, IM$Qupstream)
IM$LengthHydro <- c(1000, 1500)

PSDini <- ParamSD
PSDini[1] <- PSDini[1] / 2 # 2 time step delay
Ind_Run1 <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1990-01-01"),
                which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1990-12-31"))
Ind_Run2 <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1991-01-01"),
                which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1991-12-31"))

# 1990
RunOptions1 <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
194
195
                                                 InputsModel = IM,
                                                 IndPeriod_Run = Ind_Run1))
196
OutputsModel1 <- RunModel(InputsModel = IM,
197
198
                          RunOptions = RunOptions1, Param = PSDini,
                          FUN_MOD = RunModel_GR4J)
199
200
# 1990-1991
RunOptions <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
201
202
                                                InputsModel = IM,
                                                IndPeriod_Run = c(Ind_Run1, Ind_Run2)))
203
OutputsModel <- RunModel(InputsModel = IM,
204
205
206
                         RunOptions = RunOptions,
                         Param = PSDini,
                         FUN_MOD = RunModel_GR4J)
207

208
209
210
test_that("Warm start should give same result as warmed model", {
  # Warm start 1991
  RunOptions2 <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
211
                                  InputsModel = IM, IndPeriod_Run = Ind_Run2,
212
213
                                  IndPeriod_WarmUp = 0L,
                                  IniStates = OutputsModel1$StateEnd)
214
  OutputsModel2 <- RunModel(InputsModel = IM,
215
216
217
                            RunOptions = RunOptions2,
                            Param = PSDini,
                            FUN_MOD = RunModel_GR4J)
218
219
220
  # Compare 1991 Qsim from warm started and from 1990-1991
  names(OutputsModel2$Qsim) <- NULL
  expect_equal(OutputsModel2$Qsim, OutputsModel$Qsim[366:730])
221
222
223
224
225
226
227
228
229
230
231
})

test_that("Error on Wrong length of iniState$SD", {
  OutputsModel1$StateEnd$SD[[1]] <- c(1,1)
  RunOptions2 <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
                                  InputsModel = IM, IndPeriod_Run = Ind_Run2,
                                  IndPeriod_WarmUp = 0L,
                                  IniStates = OutputsModel1$StateEnd)
  expect_error(RunModel(InputsModel = IM, RunOptions = RunOptions2, Param = PSDini, FUN_MOD = RunModel_GR4J)
  )
})