test-RunModel_LAG.R 8.54 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
                               InputsModel = InputsModel,
55
                               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
137
  expect_equal(QlsSdSim - QlsGR4Only, QlsUpstLagObs)
})


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

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", {
166
  Qm3GR4Only <- OutputsGR4JOnly$Qsim * BasinAreas[2L] * 1e3
167
  # Specify that upstream flow is not related to an area
168
  InputsModel$BasinAreas <- c(NA, BasinAreas[2L])
169
  # Convert upstream flow to m3/day
170
171
172
173
  InputsModel$Qupstream <- matrix(Qupstream, ncol = 1) * BasinAreas[1L] * 1e3
  
  OutputsSD <- RunModel(InputsModel, RunOptions, Param = ParamSD, FUN_MOD = RunModel_GR4J)
  
174
  expect_false(any(is.na(OutputsSD$Qsim)))
175
176
177
178
  
  Qm3SdSim <- OutputsSD$Qsim * sum(InputsModel$BasinAreas, na.rm = TRUE) * 1e3
  Qm3UpstLagObs <- c(0, InputsModel$Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]])
  
179
180
  expect_equal(Qm3SdSim - Qm3GR4Only, Qm3UpstLagObs)
})
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212

test_that("Warm start should give same result as warmed model", {
  InputsModel$BasinAreas <- rep(BasinInfo$BasinArea, 3)
  InputsModel$Qupstream <- cbind(InputsModel$Qupstream, InputsModel$Qupstream)
  InputsModel$LengthHydro <- c(1000, 1500)
  
  ParamSD[1] <- ParamSD[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-1991
  RunOptions <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
                                                  InputsModel = InputsModel, IndPeriod_Run = c(Ind_Run1, Ind_Run2)))
  OutputsModel <- RunModel(InputsModel = InputsModel,
                                RunOptions = RunOptions, Param = ParamSD, FUN_MOD = RunModel_GR4J)
  # 1990
  RunOptions1 <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
                                                   InputsModel = InputsModel, IndPeriod_Run = Ind_Run1))
  OutputsModel1 <- RunModel(InputsModel = InputsModel,
                                 RunOptions = RunOptions1, Param = ParamSD, FUN_MOD = RunModel_GR4J)
  # Warm start 1991
  RunOptions2 <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
                                  InputsModel = InputsModel, IndPeriod_Run = Ind_Run2, 
                                  IndPeriod_WarmUp = 0L,
                                  IniStates = OutputsModel1$StateEnd)
  OutputsModel2 <- RunModel(InputsModel = InputsModel,
                                 RunOptions = RunOptions2, Param = ParamSD, FUN_MOD = RunModel_GR4J)
  # Compare 1991 Qsim from warm started and from 1990-1991
  names(OutputsModel2$Qsim) <- NULL
  expect_equal(OutputsModel2$Qsim, OutputsModel$Qsim[366:730])
})