server.R 25.3 KB
Newer Older
1
2
# server.R

3
shinyServer(function(input, output, session) {
4
5
  
  
6
7
8
9
10
11
12
13
14
15
16
17
  ## --------------- List of input names
  
  getInputs <- reactive({
    inputList <- sort(names(reactiveValuesToList(input)))
    inputList <- inputList[!grepl("^dy", inputList)]
    inputList <- inputList[!grepl("^CalButton", inputList)]
    inputList <- c(inputList, "DownloadTab", "DownloadPlot")
    return(inputList)
  })
  
  
  
18
19
  ## --------------- Data preparation
  
20
  getPrep <- reactive({
21
    
22
    TMGR  <- .TypeModelGR(input$HydroModel)
23
    PARAM <- c(input$X1, input$X2, input$X3, input$X4, input$X5, input$X6)[seq_len(TMGR$NbParam)]
24
    
25
    if (input$SnowModel == "CemaNeige") {
26
27
      PARAM <- c(PARAM, input$C1, input$C2)
    }
28
    
29
30
    OBS <- ObsGR(ObsBV = get(input$Dataset), HydroModel = input$HydroModel,
                 CemaNeige = input$SnowModel == "CemaNeige",
31
32
                 Precip = .ShinyGR.args$Precip, PotEvap = .ShinyGR.args$PotEvap,
                 Qobs = get(input$Dataset), TempMean = .ShinyGR.args$TempMean, 
33
34
                 ZInputs = .ShinyGR.args$ZInputs, HypsoData = .ShinyGR.args$HypsoData,
                 NLayers = .ShinyGR.args$NLayers)
35
    
36
37
38
39
40
    WUPPER <- c(OBS$InputsModel$DatesR[1L], input$Period[1]-.TypeModelGR(OBS)$TimeLag)
    if (WUPPER[2] < WUPPER[1]) {
      WUPPER[2] <- WUPPER[1]
    }
    
41
42
43
44
45
46
47
    ## Enable or disable automatic calibration (if there is Qobs or not)
    isQobs <- !all(is.na(OBS$Qobs[OBS$InputsModel$Dates >= input$Period[1] & OBS$InputsModel$Dates <= input$Period[2]]))
    if (isQobs) {
      shinyjs::enable("CalButton")
    } else {
      shinyjs::disable("CalButton")
    }
48
    return(list(TMGR = TMGR, OBS = OBS, WUPPER = WUPPER))
49
    
50
  })
51
  
52
  
53
54
55
  
  ## --------------- Calibration
  
56
  ## If the user calibrate the model
57
  CAL_click <- reactiveValues(valueButton = 0)
58
  
unknown's avatar
unknown committed
59
  
60
  ## Automatic calibration
61
  observeEvent(input$CalButton, {
62
63
64
65
66
67
    
    ## Desable all inputs during automatic calibration
    lapply(getInputs(), shinyjs::disable)
    shinyjs::disable("CalButton")
    
    ## Model calibration
68
69
70
    CAL_opt <- list(Crit    = gsub(" .*", "", input$TypeCrit),
                    Transfo = gsub("(\\D{3} \\[)(\\w{0,4})(\\W*Q\\W*\\])", "\\2", input$TypeCrit))
    
71
    CAL     <- CalGR(ObsGR = getPrep()$OBS, CalCrit = CAL_opt$Crit, transfo = CAL_opt$Transfo,
72
73
                     WupPer = substr(getPrep()$WUPPER, 1, 10),
                     CalPer = substr(c(input$Period[1], input$Period[2]), 1, 10), verbose = FALSE)
74
    PARAM   <- CAL$OutputsCalib$ParamFinalR
75
    
76
77
78
79
    updateSliderInput(session, inputId = "X1", value = PARAM[1L])
    updateSliderInput(session, inputId = "X2", value = PARAM[2L])
    updateSliderInput(session, inputId = "X3", value = PARAM[3L])
    updateSliderInput(session, inputId = "X4", value = PARAM[4L])
80
    if (getPrep()$TMGR$NbParam >= 5) {
81
82
      updateSliderInput(session, inputId = "X5", value = PARAM[5L])
    }
83
    if (getPrep()$TMGR$NbParam >= 6) {
84
85
      updateSliderInput(session, inputId = "X6", value = PARAM[6L])
    }
86
    if (input$SnowModel == "CemaNeige") {
87
88
89
      updateSliderInput(session, inputId = "C1", value = PARAM[length(PARAM)-1])
      updateSliderInput(session, inputId = "C2", value = PARAM[length(PARAM)])
    }
90
91
    updateActionButton(session, inputId = "CalButton", label = "Model calibrated", icon = icon("check"))
    CAL_click$valueButton <- 1
92
93
94
95
    
    ## Enable caliration
    shinyjs::enable("CalButton")
  }, priority = +20)
96
  
unknown's avatar
unknown committed
97
  
98
  ## Manual calibration
99
  observeEvent({input$Dataset ; input$HydroModel ; input$SnowModel ;
100
101
102
103
104
105
106
107
    input$X1 ; input$X2 ; input$X3 ; input$X4 ; input$X5 ; input$X6 ;
    input$TypeCrit ; input$Period}, {
      CAL_click$valueButton <-  CAL_click$valueButton - 1
      CAL_click$valueButton <- ifelse(CAL_click$valueButton < -1, -1, CAL_click$valueButton)
      if (CAL_click$valueButton < 0) {
        updateActionButton(session, inputId = "CalButton", label = "Run", icon = icon("refresh"))
      }
      
108
109
      ## Enable all inputs except automatic calibration
      lapply(getInputs(), shinyjs::enable)
110
111
112
    })
  
  
113
  
114
  ## --------------- Simulation
115
116
  
  getRES <- reactive({
117
    PARAM <- c(input$X1, input$X2, input$X3, input$X4, input$X5, input$X6)[seq_len(getPrep()$TMGR$NbParam)]
118
    if (input$SnowModel == "CemaNeige") {
119
      PARAM <- c(PARAM, input$C1, input$C2)
120
    }
121
    
122
123
124
125
126
    ## Simulated flows computation
    SIM_opt <- list(Crit    = c("ErrorCrit_NSE", "ErrorCrit_KGE"),
                    Transfo = c("NO", "sqrt", "log"))
    SIM <- lapply(SIM_opt$Crit, function(iCRIT) {
      SIM_transfo <- lapply(SIM_opt$Transfo, function(iTRSF) {
127
        iTRSF <- gsub("NO", "", iTRSF)
128
129
        iSIM  <- SimGR(ObsGR = getPrep()$OBS, Param = PARAM,
                       WupPer = substr(getPrep()$WUPPER, 1, 10),
130
131
132
                       SimPer = substr(c(input$Period[1], input$Period[2]), 1, 10),
                       transfo = iTRSF, verbose = FALSE)
        iCRIT <- ErrorCrit(InputsCrit = iSIM$OptionsCrit, OutputsModel = iSIM$OutputsModel, FUN_CRIT = get(iCRIT), verbose = FALSE)
133
134
135
        iCRIT <- iCRIT[c("CritName", "CritValue")]
        return(list(SIM = iSIM, CRIT = iCRIT))
      })
136
      names(SIM_transfo) <- SIM_opt$Transfo
137
      return(SIM_transfo)
138
    })
139
    names(SIM) <- SIM_opt$Crit
140
    
141
    ## Criteria computation
142
    CRIT <- lapply(SIM, function(iCRIT) {
143
      lapply(SIM_opt$Transfo, function(iTRSF) {
144
145
        iCRIT[[iTRSF]][["CRIT"]]
      })
146
    })
147
    CRIT <- as.data.frame(matrix(na.omit(unlist(CRIT)), ncol = 2, byrow = TRUE), stringsAsFactors = FALSE)
148
    colnames(CRIT) <- c("Criterion", "Value")
149
    rownames(CRIT) <- NULL    
150
151
    CRIT$Value     <- as.numeric(CRIT$Value)
    CRIT$Criterion <- gsub("\\[", " [", CRIT$Criterion)
152
153
154
155
    
    .GlobalEnv$.ShinyGR.hist[[length(.GlobalEnv$.ShinyGR.hist)+1]] <- list(Qsim = SIM$ErrorCrit_KGE$NO$SIM$OutputsModel$Qsim,
                                                                           Param = PARAM,
                                                                           TypeModel = SIM$ErrorCrit_KGE$NO$SIM$TypeModel)
156
    .GlobalEnv$.ShinyGR.hist <- .GlobalEnv$.ShinyGR.hist[!(duplicated(sapply(.GlobalEnv$.ShinyGR.hist, function(x) sum(x$Param)), fromLast = TRUE) & 
157
                                                             duplicated(sapply(.GlobalEnv$.ShinyGR.hist, function(x) x$TypeModel ), fromLast = TRUE))]
158
    .GlobalEnv$.ShinyGR.hist <- tail(.GlobalEnv$.ShinyGR.hist, n = 2)
159
    if (length(.GlobalEnv$.ShinyGR.hist) == 2 &  is.null(names(.GlobalEnv$.ShinyGR.hist[[1]]))) {
160
      .GlobalEnv$.ShinyGR.hist[[1]] <- NULL
161
    }
162
163
    if (length(.GlobalEnv$.ShinyGR.hist) == 2 & !is.null(names(.GlobalEnv$.ShinyGR.hist[[1]]))) {
      if (length(.GlobalEnv$.ShinyGR.hist[[1]]$Qsim) != length(.GlobalEnv$.ShinyGR.hist[[2]]$Qsim)) {
164
165
        OBSold <- getPrep()$OBS
        OBSold$TypeModel <- .GlobalEnv$.ShinyGR.hist[[1]]$TypeModel
166
167
        if (.TypeModelGR(OBSold)$CemaNeige & !.TypeModelGR(getPrep()$OBS)$CemaNeige) { # present: No CemaNeige ; old: CemaNeige
          OBSold <- ObsGR(ObsBV = get(input$Dataset), HydroModel = .TypeModelGR(OBSold)$NameModel,
168
169
170
171
172
                          CemaNeige = .TypeModelGR(OBSold)$CemaNeige,
                          Precip = .ShinyGR.args$Precip, PotEvap = .ShinyGR.args$PotEvap,
                          Qobs = get(input$Dataset), TempMean = .ShinyGR.args$TempMean, 
                          ZInputs = .ShinyGR.args$ZInputs, HypsoData = .ShinyGR.args$HypsoData,
                          NLayers = .ShinyGR.args$NLayers)
173
        }
174
        SIMold <- SimGR(ObsGR = OBSold,
175
176
177
178
179
                        Param = .GlobalEnv$.ShinyGR.hist[[1]]$Param,
                        WupPer = substr(getPrep()$WUPPER, 1, 10),
                        SimPer = substr(c(input$Period[1], input$Period[2]), 1, 10),
                        verbose = FALSE)
        .GlobalEnv$.ShinyGR.hist[[1]]$Qsim <- SIMold$OutputsModel$Qsim
180
      }
181
182
    }
    
183
    return(list(PARAM = PARAM, SIM = SIM$ErrorCrit_KGE$NO$SIM, SIMold = .GlobalEnv$.ShinyGR.hist, Crit = CRIT))
184
185
    
  })
186
  
187
  
188
189
190
191
  
  ## --------------- Plot
  
  ## Choice
192
193
194
195
196
197
198
  getPlotType <- reactive({
    switch(input$PlotType,
           "Model performance" = 1,
           "Flow time series"  = 2,
           "State variables"   = 3,
           "Model diagram"     = 4)
  })
199
  
200
201
202
  
  ## Models available considering the plot type
  observe({
203
    if (getPlotType() == 4) {
204
      updateSelectInput(session, inputId = "HydroModel", choice = c("GR4J", "GR5J"), selected = input$HydroModel)
205
      updateSelectInput(session, inputId = "SnowModel" , choice = c("None"))
206
    } else {
207
      updateSelectInput(session, inputId = "HydroModel", choice = c("GR4J", "GR5J", "GR6J"), selected = input$HydroModel)
208
      updateSelectInput(session, inputId = "SnowModel" , choice = c("None", "CemaNeige")   , selected = input$SnowModel)
209
    }
210
211
  })
  
unknown's avatar
unknown committed
212
  
213
214
  ## Plots available considering the model type
  observe({
215
    if (input$HydroModel == "GR6J") {
216
217
218
219
220
221
222
      updateSelectInput(session, inputId = "PlotType",
                        choice = c("Flow time series", "Model performance", "State variables"),
                        selected = input$PlotType)
    } else {
      updateSelectInput(session, inputId = "PlotType",
                        choice = c("Flow time series", "Model performance", "State variables", "Model diagram"),
                        selected = input$PlotType)
223
    }
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
  })
  
  
  ## Period slider responds to changes in the selected/zoomed dateWindow 
  observeEvent({input$dyPlotTS_date_window ; input$dyPlotSVs_date_window ; input$dyPlotMDp_date_window}, {
    if (!is.null(input$dyPlotTS_date_window) && getPlotType() == 2) {
      dateWindow <- as.POSIXct(strftime(input$dyPlotTS_date_window, "%Y-%m-%d %H:%M:%S"), tz = "UTC")
    }
    if (!is.null(input$dyPlotSVs_date_window) && getPlotType() == 3) {
      dateWindow <- as.POSIXct(strftime(input$dyPlotSVs_date_window, "%Y-%m-%d %H:%M:%S"), tz = "UTC")
    }
    if (!is.null(input$dyPlotMDp_date_window) && getPlotType() == 4) {
      dateWindow <- as.POSIXct(strftime(input$dyPlotMDp_date_window, "%Y-%m-%d %H:%M:%S"), tz = "UTC")
    }
    if (exists("dateWindow")) {
      updateSliderInput(session, inputId = "Period",
240
                        value = dateWindow + .TypeModelGR(input$HydroModel)$TimeLag)
241
    }
unknown's avatar
unknown committed
242
243
  })
  
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
  
  ## Reset period slider responds to dygraphs to mouse clicks
  observeEvent({input$dyPlotTS_click}, {
    updateSliderInput(session, inputId = "Period",
                      value = as.POSIXct(.ShinyGR.args$SimPer, tz = "UTC"))
  }, priority = +10)
  observeEvent({input$dyPlotSVs_click}, {
    updateSliderInput(session, inputId = "Period",
                      value = as.POSIXct(.ShinyGR.args$SimPer, tz = "UTC"))
  }, priority = +10)
  observeEvent({input$dyPlotSVq_click}, {
    updateSliderInput(session, inputId = "Period",
                      value = as.POSIXct(.ShinyGR.args$SimPer, tz = "UTC"))
  }, priority = +10)
  observeEvent({input$dyPlotMDp_click}, {
    updateSliderInput(session, inputId = "Period",
                      value = as.POSIXct(.ShinyGR.args$SimPer, tz = "UTC"))
  }, priority = +10)
  observeEvent({input$dyPlotMDe_click}, {
    updateSliderInput(session, inputId = "Period",
                      value = as.POSIXct(.ShinyGR.args$SimPer, tz = "UTC"))
  }, priority = +10)
  observeEvent({input$dyPlotMDq_click}, {
    updateSliderInput(session, inputId = "Period",
                      value = as.POSIXct(.ShinyGR.args$SimPer, tz = "UTC"))
  }, priority = +10)
  
  
unknown's avatar
unknown committed
272
273
274
  ## Target date slider
  observe({
    updateSliderInput(session, inputId = "Event",
275
                      min = input$Period[1L] + .TypeModelGR(input$HydroModel)$TimeLag,
unknown's avatar
unknown committed
276
                      max = input$Period[2L])
277
  }) 
unknown's avatar
unknown committed
278
279
  
  
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
  ## Graphical parameters
  getPlotPar <- reactive({
    if (.GlobalEnv$.ShinyGR.args$theme == "Cyborg") {
      col_bg <- "black"
      col_fg <- "white"
      par(bg = col_bg, fg = col_fg, col.axis = col_fg, col.lab = col_fg)
    } else if (.GlobalEnv$.ShinyGR.args$theme == "Flatly") {
      col_bg <- "#2C3E50"
      col_fg <- "white"
      par(bg = col_bg, fg = col_fg, col.axis = col_bg, col.lab = col_bg)
    } else {
      col_bg <- "white"
      col_fg <- "black"
      par(bg = col_bg , fg = col_fg)
    }
    return(list(col_bg = col_bg, col_fg = col_fg, par = par(no.readonly = TRUE)))
  })
297
  
298
299
300
301
  
  ## Plot model performance
  output$stPlotMP <- renderPlot({
    OutputsModel <- getRES()$SIM$OutputsModel
302
    IndPlot <- which(OutputsModel$DatesR >= input$Period[1L] & OutputsModel$DatesR <= input$Period[2L])
303
304
    par(getPlotPar()$par)
    par(cex.axis = 1.2)
305
306
307
    if (input$SnowModel != "CemaNeige") {
      par(oma = c(20, 0, 0, 0))
    }
308
    plot(OutputsModel, Qobs = getRES()$SIM$Qobs, IndPeriod_Plot = IndPlot, cex.lab = 1.2, cex.axis = 1.4, cex.leg = 1.4)
309
  }, bg = "transparent")
310
  
unknown's avatar
unknown committed
311
  
312
  ## Plot flow time series
313
  output$dyPlotTS <- dygraphs::renderDygraph({
314
315
316
317
318
319
    if (length(.GlobalEnv$.ShinyGR.hist) == 2 & input$ShowOldQsim == "Yes") {
      QsimOld <- getRES()$SIMold[[1]]$Qsim
    } else {
      QsimOld <- NULL
    }
    # Qold <- ifelse(format(OutputsModel2$DatesR, "%Y%m%d") > format(input$Event, "%Y%m%d"), NA, QsimOld[IndPlot])
320
    op <- getPlotPar()$par
321
322
    dg1 <- dyplot(getRES()$SIM, Qsup = QsimOld, Qsup.name = "Qold", RangeSelector = FALSE, LegendShow = "auto",
                  col.Q = c(op$fg, "orangered", "grey"), col.Precip = "#428BCA")
323
324
    dg1 <- dygraphs::dyOptions(dg1, axisLineColor = op$fg, axisLabelColor = op$fg, retainDateWindow = FALSE)
    dg1 <- dygraphs::dyLegend(dg1, show = "follow", width = 325)
325
326
  })
  
unknown's avatar
unknown committed
327
  
328
  ## Plot state variables stores
329
  output$dyPlotSVs <- dygraphs::renderDygraph({
330
331
332
333
    OutputsModel <- getRES()$SIM$OutputsModel
    data <- data.frame(DatesR = OutputsModel$DatesR,
                       prod.  = OutputsModel$Prod,
                       rout.  = OutputsModel$Rout)
334
    data.xts <- xts::xts(data[, -1L], order.by = data$DatesR)
335
336
    
    op <- getPlotPar()$par
337
338
339
340
341
342
    dg2 <- dygraphs::dygraph(data.xts, group = "state_var", ylab = "store [mm]")
    dg2 <- dygraphs::dyOptions(dg2, colors = c("#00008B", "#008B8B"),
                               fillGraph = TRUE, fillAlpha = 0.3,
                               drawXAxis = FALSE, axisLineColor = op$fg, axisLabelColor = op$fg, retainDateWindow = FALSE)
    dg2 <- dygraphs::dyLegend(dg2, show = "always", width = 325)
    dg2 <- dygraphs::dyCrosshair(dg2, direction = "vertical")
343
344
  })
  
unknown's avatar
unknown committed
345
  
346
  ## Plot state variables Q
347
  output$dyPlotSVq <- dygraphs::renderDygraph({
348
349
350
351
    OutputsModel <- getRES()$SIM$OutputsModel
    IndPlot <- which(OutputsModel$DatesR >= input$Period[1L] & OutputsModel$DatesR <= input$Period[2L])
    OutputsModel2 <- sapply(OutputsModel[seq_len(which(names(OutputsModel) == "Qsim"))], function(x) x[IndPlot])
    OutputsModel2 <- c(OutputsModel2, Qobs = list(getRES()$SIM$Qobs[IndPlot]))
352
    
353
354
    data <- data.frame(DatesR = OutputsModel2$DatesR,
                       Qr     = OutputsModel2$QR,
355
                       Qd     = OutputsModel2$QD,
356
357
                       Qsim   = OutputsModel2$Qsim,
                       Qobs   = OutputsModel2$Qobs)
358
    if (input$HydroModel == "GR6J") {
359
      data$QrExp <- OutputsModel2$QRExp
360
    } else {
361
      data$QrExp <- NA
362
    }
363
    data.xts <- xts::xts(data[, -1L], order.by = data$DatesR)
364
365
    
    op <- getPlotPar()$par
366
367
368
369
370
371
372
373
374
375
    dg3 <- dygraphs::dygraph(data.xts, group = "state_var", ylab = "flow [mm/d]", main = " ")
    dg3 <- dygraphs::dyOptions(dg3, fillAlpha = 1.0,
                               axisLineColor = op$fg, axisLabelColor = op$fg,
                               titleHeight = 10, retainDateWindow = FALSE)
    dg3 <- dygraphs::dyStackedRibbonGroup(dg3, name = c("Qd", "Qr", "QrExp"),
                                          color = c("#FFD700", "#EE6300", "brown"), strokeBorderColor = "black")
    dg3 <- dygraphs::dySeries(dg3, name = "Qobs", fillGraph = FALSE, drawPoints = TRUE, color = op$fg)
    dg3 <- dygraphs::dySeries(dg3, name = "Qsim", fillGraph = FALSE, color = "orangered")
    dg3 <- dygraphs::dyCrosshair(dg3, direction = "vertical")
    dg3 <- dygraphs::dyLegend(dg3, show = "always", width = 325)
376
  })
377
  
unknown's avatar
unknown committed
378
  
379
  ## Plot model diagram precipitation
380
  output$dyPlotMDp <- dygraphs::renderDygraph({
381
382
    # barChartPrecip <- scan(file = system.file("plugins/barChartPrecip.js", package = "airGRteaching"),
    #                        what = "character", quiet = TRUE)
383
384
    data <- data.frame(DatesR  = getRES()$SIM$OutputsModel$DatesR,
                       precip. = getRES()$SIM$OutputsModel$Precip)
385
    data.xts <- xts::xts(data[, -1L, drop = FALSE], order.by = data$DatesR)
386
    
387
388
389
390
391
392
393
    dg4 <- dygraphs::dygraph(data.xts, group = "mod_diag", ylab = "precip. [mm/d]")
    dg4 <- dygraphs::dyOptions(dg4, colors = "#428BCA", drawXAxis = FALSE, retainDateWindow = FALSE)
    dg4 <- dygraphs::dyBarSeries(dg4, name = "precip.")
    dg4 <- dygraphs::dyAxis(dg4, name = "y", valueRange = c(max(data.xts[, "precip."], na.rm = TRUE), -1e-3))
    dg4 <- dygraphs::dyEvent(dg4, input$Event, color = "orangered")
    dg4 <- dygraphs::dyLegend(dg4, show = "onmouseover", width = 225)
    dg4 <- dygraphs::dyCrosshair(dg4, direction = "vertical")
394
  })
395
  
unknown's avatar
unknown committed
396
  
397
  ## Plot model diagram ETP
398
  output$dyPlotMDe <- dygraphs::renderDygraph({
399
400
    op <- getPlotPar()$par
    data <- data.frame(DatesR = getRES()$SIM$OutputsModel$DatesR,
401
                       PET    = getRES()$SIM$OutputsModel$PotEvap)
402
    data.xts <- xts::xts(data[, -1L, drop = FALSE], order.by = data$DatesR)
403
    
404
405
406
407
408
409
410
411
    dg5 <- dygraphs::dygraph(data.xts, group = "mod_diag", ylab = "PET [mm/d]", main = " ")
    dg5 <- dygraphs::dyOptions(dg5, colors = "#A4C400", drawPoints = TRUE,
                               strokeWidth = 0, pointSize = 2, drawXAxis = FALSE,
                               axisLineColor = op$fg, axisLabelColor = op$fg,
                               titleHeight = 10, retainDateWindow = FALSE)
    dg5 <- dygraphs::dyEvent(dg5, input$Event, color = "orangered")
    dg5 <- dygraphs::dyLegend(dg5, show = "onmouseover", width = 225)
    dg5 <- dygraphs::dyCrosshair(dg5, direction = "vertical")
412
  })
413
  
unknown's avatar
unknown committed
414
  
415
  ## Plot model diagram flow
416
  output$dyPlotMDq <- dygraphs::renderDygraph({
417
418
    if (length(.GlobalEnv$.ShinyGR.hist) == 2 & input$ShowOldQsim == "Yes") {
      QsimOld <- getRES()$SIMold[[1]]$Qsim
419
420
421
    } else {
      QsimOld <- NA
    }
422
423
424
425
    OutputsModel <- getRES()$SIM$OutputsModel
    IndPlot <- which(OutputsModel$DatesR >= input$Period[1L] & OutputsModel$DatesR <= input$Period[2L])
    OutputsModel2 <- sapply(OutputsModel[seq_len(which(names(OutputsModel) == "Qsim"))], function(x) x[IndPlot])
    OutputsModel2 <- c(OutputsModel2, Qobs = list(getRES()$SIM$Qobs[IndPlot]))
unknown's avatar
unknown committed
426
    OutputsModel2$Qsim <- ifelse(format(OutputsModel2$DatesR, "%Y%m%d") > format(input$Event, "%Y%m%d"), NA, OutputsModel2$Qsim)
427
    OutputsModel2$Qold <- ifelse(format(OutputsModel2$DatesR, "%Y%m%d") > format(input$Event, "%Y%m%d"), NA, QsimOld[IndPlot])
428
    
429
430
431
432
    data <- data.frame(DatesR  = OutputsModel2$DatesR,
                       Qobs    = OutputsModel2$Qobs,
                       Qsim    = OutputsModel2$Qsim,
                       QsimOld = OutputsModel2$Qold)
433
    data.xts <- xts::xts(data[, -1L, drop = FALSE], order.by = data$DatesR)
434
435
    
    op <- getPlotPar()$par
436
437
438
439
440
441
442
443
444
    dg6 <- dygraphs::dygraph(data.xts, group = "mod_diag", ylab = "flow [mm/d]", main = " ")
    dg6 <- dygraphs::dyOptions(dg6, colors = c(op$fg, "grey", "orangered"), drawPoints = TRUE,
                               axisLineColor = op$fg, axisLabelColor = op$fg,
                               titleHeight = 10, retainDateWindow = FALSE)
    dg6 <- dygraphs::dySeries(dg6, name = "QsimOld", drawPoints = FALSE, strokePattern = "dashed")
    dg6 <- dygraphs::dySeries(dg6, name = "Qsim"   , drawPoints = FALSE)
    dg6 <- dygraphs::dyEvent(dg6, input$Event, color = "orangered")
    dg6 <- dygraphs::dyLegend(dg6, show = "onmouseover", width = 225)
    dg6 <- dygraphs::dyCrosshair(dg6, direction = "vertical")
445
  })
446
  
unknown's avatar
unknown committed
447
  
448
449
450
451
452
453
  ## Plot model diagram chart
  output$stPlotMD <- renderPlot({
    OutputsModel <- getRES()$SIM$OutputsModel
    IndPlot <- which(OutputsModel$DatesR >= input$Period[1L] & OutputsModel$DatesR <= input$Period[2L])
    OutputsModel2 <- sapply(OutputsModel[seq_len(which(names(OutputsModel) == "Qsim"))], function(x) x[IndPlot])
    OutputsModel2 <- c(OutputsModel2, Qobs = list(getRES()$SIM$Qobs[IndPlot]))
454
    
455
456
457
    par(getPlotPar()$par)
    airGRteaching:::DiagramGR(OutputsModel = OutputsModel2, Param = getRES()$PARAM,
                              SimPer = input$Period, EventDate = input$Event,
458
                              HydroModel = input$HydroModel)
459
  }, bg = "transparent")
460
461
  
  
462
  
463
464
  ## --------------- Criteria table
  
465
  output$Criteria <- renderTable({
466
    ## Table created in order to choose order the criteria in the table output
467
    tabCrit_gauge <- data.frame(Criterion = c("NSE [Q]", "NSE [sqrt(Q)]", "NSE [log(Q)]",
468
                                              "KGE [Q]", "KGE [sqrt(Q)]", "KGE [log(Q)]"),
469
                                ID        = 1:6, stringsAsFactors = FALSE)
470
    
471
    tabCrit_out <- merge(tabCrit_gauge, getRES()$Crit, by = "Criterion", all.x = TRUE)
472
473
    tabCrit_out <- tabCrit_out[order(tabCrit_out$ID), ]
    tabCrit_out <- tabCrit_out[, !colnames(tabCrit_out) %in% "ID"]
474
475
476
    
    ## Color the cell of the crietaia uses during the calibration
    if (CAL_click$valueButton >= 0) {
477
      CellCol <- '<div style="color: #FFFFFF; background-color: #A4C400; border: 5px solid #A4C400; position:relative; top: 0px; left: 5px; padding: 0px; margin: -5px -0px -8px -10px;">
478
<span>9999</span></div>'
479
      CellCol_id  <- which(tabCrit_out$Criterion == input$TypeCrit)
480
481
482
483
484
      tabCrit_out[CellCol_id, 1] <- gsub("9999", tabCrit_out[CellCol_id, 1], CellCol)
    }
    
    return(tabCrit_out)
  }, sanitize.text.function = function(x) x)
485
486
487
488
489
  
  
  
  ## --------------- Download buttons
  
490
  ## Download simulation table
491
492
493
  output$DownloadTab <- downloadHandler(
    filename = function() {
      filename <- "TabSim"
494
      filename <- sprintf("airGR_%s_%s.csv", filename, gsub("(.*)( )(\\d{2})(:)(\\d{2})(:)(\\d{2})", "\\1_\\3h\\5m\\7s", Sys.time()))
495
496
497
498
499
500
501
    },
    content = function(file) {
      OBS <- getPrep()$OBS
      SIM <- getRES()$SIM
      if (input$SnowModel != "CemaNeige") {
        PrecipSim <- NA
        FracSolid <- NA
502
        TempMean  <- NA
503
504
505
506
507
508
509
      } else {
        PrecipSol <- rowMeans(as.data.frame(OBS$InputsModel$LayerPrecip) * as.data.frame(OBS$InputsModel$LayerFracSolidPrecip), na.rm = TRUE)
        PrecipSim <- rowMeans(as.data.frame(OBS$InputsModel$LayerPrecip), na.rm = TRUE)
        FracSolid <- PrecipSol / PrecipSim
        FracSolid <- ifelse(is.na(FracSolid)  & PrecipSol == 0 & PrecipSim == 0, 0, FracSolid)
        PrecipSim <- PrecipSim[SIM$OptionsSimul$IndPeriod_Run]
        FracSolid <- FracSolid[SIM$OptionsSimul$IndPeriod_Run]
510
511
512
        FracSolid <- round(FracSolid, digits = 3)
        TempMean  <- rowMeans(as.data.frame(OBS$InputsModel$LayerTempMean), na.rm = TRUE)
        TempMean  <- TempMean[SIM$OptionsSimul$IndPeriod_Run]
513
      }
514
515
516
517
518
519
520
521
522
523
524
525
      TabSim <- data.frame(Dates                     = SIM$OutputsModel$DatesR,
                           PotEvap                   = SIM$OutputsModel$PotEvap,
                           PrecipObs                 = SIM$OutputsModel$Precip,
                           PrecipSim_CemaNeige       = PrecipSim,
                           PrecipFracSolid_CemaNeige = FracSolid,
                           TempMeanSim_CemaNeige     = TempMean,
                           Qobs                      = SIM$OptionsCrit$Qobs,
                           Qsim                      = SIM$OutputsModel$Qsim)
      colnames(TabSim) <- sprintf("%s [%s]", colnames(TabSim), c("-", rep("mm",3), "-", "°C", rep("mm", 2)))
      colnames(TabSim) <- ifelse(grepl("mm", colnames(TabSim)),
                                 gsub("mm", paste0("mm/", .TypeModelGR(OBS)$TimeUnit),colnames(TabSim)),
                                 colnames(TabSim))
526
527
528
529
      write.table(TabSim, file = file, row.names = FALSE, sep = ";")
    }
  )
  
530
  ## Download plots
531
532
533
  output$DownloadPlot <- downloadHandler(
    filename = function() {
      filename <- switch(input$PlotType,
534
535
536
537
                         "Model performance" = "PlotModelPerf",
                         "Flow time series"  = "PlotFlowTimeSeries",
                         "State variables"   = "PlotStateVar",
                         "Model diagram"     = "PlotModelDiag")
538
      filename <- sprintf("airGR_%s_%s.png", filename, gsub("(.*)( )(\\d{2})(:)(\\d{2})(:)(\\d{2})", "\\1_\\3h\\5m\\7s", Sys.time()))
539
540
541
    },
    content = function(file) {
      k <- 1.75
542
543
544
545
546
547
548
549
550
      ParamTitle <- c("X1", "X2", "X3", "X4", "X5", "X6")[seq_len(getPrep()$TMGR$NbParam)]
      if (input$SnowModel == "CemaNeige") {
        ParamTitle <- c(ParamTitle, "C1", "C2")
      }
      ParamTitle <- paste(ParamTitle, getRES()$PARAM, sep = " = ", collapse = ", ")
      PngTitle <- sprintf("%s - %s/%s\n%s\n%s", input$Dataset,
                          input$HydroModel, ifelse(input$SnowModel == "CemaNeige", "CemaNeige", "No snow model"),
                          paste0(input$Period, collapse = " - "),
                          ParamTitle)
551
552
      if (getPlotType() == 1) {
        png(filename = file, width = 1000*k,  height = ifelse(input$SnowModel != "CemaNeige", 700*k, 1100*k), pointsize = 14, res = 150)
553
        par(oma = c(0, 0, 4, 0))
554
        plot(getRES()$SIM)
555
        mtext(text = PngTitle, side = 3, outer = TRUE, cex = 0.8, line = 1.2)
556
557
558
559
        dev.off()
      }
      if (getPlotType() == 2) {
        png(filename = file, width = 1000*k, height = 600*k, pointsize = 14, res = 150)
560
        par(oma = c(0, 0, 4, 0))
561
        plot(getRES()$SIM, which = c( "Precip", "Flows"))
562
        mtext(text = PngTitle, side = 3, outer = TRUE, cex = 0.8, line = 1.2)
563
564
565
566
567
568
        dev.off()
      }
    }
  )
  
  
569
570
})