server.R 39.4 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
18
  ## --------------- 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)
  })
  
  
  
19
20
  ## --------------- Data preparation
  
21
  getPrep <- reactive({
22
    
23
    TMGR  <- .TypeModelGR(input$HydroModel)
24
    PARAM <- c(input$X1, input$X2, input$X3, input$X4, input$X5, input$X6)[seq_len(TMGR$NbParam)]
25
    
26
    if (input$SnowModel == "CemaNeige") {
27
28
      PARAM <- c(PARAM, input$C1, input$C2)
    }
29
    if (input$Dataset == "Unnamed watershed") {
30
      ObsDF <- NULL
31
    } else {
32
33
      # ObsDF <- get(input$Dataset)
      ObsDF <- .ShinyGR.args$ObsDF[[input$Dataset]]
34
    }
35
    PREP <- PrepGR(ObsDF = ObsDF,
36
37
38
39
40
41
42
43
44
45
                   DatesR = .ShinyGR.args$DatesR,
                   Precip = .ShinyGR.args$Precip, PotEvap = .ShinyGR.args$PotEvap,
                   Qobs = .ShinyGR.args$Qobs, TempMean = .ShinyGR.args$TempMean, 
                   ZInputs = .ShinyGR.args$ZInputs[[input$Dataset]],
                   HypsoData = .ShinyGR.args$HypsoData[[input$Dataset]],
                   NLayers = .ShinyGR.args$NLayers[[input$Dataset]],
                   HydroModel = input$HydroModel,
                   CemaNeige = input$SnowModel == "CemaNeige")
    
    WUPPER <- c(PREP$InputsModel$DatesR[1L], input$Period[1]-.TypeModelGR(PREP)$TimeLag)
46
47
48
49
    if (WUPPER[2] < WUPPER[1]) {
      WUPPER[2] <- WUPPER[1]
    }
    
50
    ## Enable or disable automatic calibration (if there is Qobs or not)
51
    isQobs <- !all(is.na(PREP$Qobs[PREP$InputsModel$Dates >= input$Period[1] & PREP$InputsModel$Dates <= input$Period[2]]))
52
53
54
    # if ( isQobs | input$Period[1L] != input$Period[2L]) {
    #   shinyjs::enable("CalButton")
    # }
55
    if (!isQobs | input$Period[1L] == input$Period[2L]) {
56
57
      shinyjs::disable("CalButton")
    }
58
    
59
    return(list(TMGR = TMGR, PREP = PREP, WUPPER = WUPPER))
60
    
61
  })
62
  
63
  
64
65
66
  
  ## --------------- Calibration
  
67
  ## If the user calibrate the model
68
  CAL_click <- reactiveValues(valueButton = 0)
69
  
unknown's avatar
unknown committed
70
  
71
  ## Automatic calibration
72
  observeEvent(input$CalButton, {
73
74
75
76
77
78
    
    ## Desable all inputs during automatic calibration
    lapply(getInputs(), shinyjs::disable)
    shinyjs::disable("CalButton")
    
    ## Model calibration
79
    CAL_opt <- list(Crit    = gsub(" .*", "", input$TypeCrit),
80
                    Transfo = gsub("1", "inv", gsub("(\\D{3} \\[)(\\w{0,4})(\\W*Q\\W*\\])", "\\2", input$TypeCrit)))
81
    CAL     <- CalGR(PrepGR = getPrep()$PREP, CalCrit = CAL_opt$Crit, transfo = CAL_opt$Transfo,
82
83
                     WupPer = substr(getPrep()$WUPPER, 1, 10),
                     CalPer = substr(c(input$Period[1], input$Period[2]), 1, 10), verbose = FALSE)
84
    PARAM   <- CAL$OutputsCalib$ParamFinalR
85
    
86
87
88
89
    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])
90
    if (getPrep()$TMGR$NbParam >= 5) {
91
92
      updateSliderInput(session, inputId = "X5", value = PARAM[5L])
    }
93
    if (getPrep()$TMGR$NbParam >= 6) {
94
95
      updateSliderInput(session, inputId = "X6", value = PARAM[6L])
    }
96
    if (input$SnowModel == "CemaNeige") {
97
98
99
      updateSliderInput(session, inputId = "C1", value = PARAM[length(PARAM)-1])
      updateSliderInput(session, inputId = "C2", value = PARAM[length(PARAM)])
    }
100
101
    updateActionButton(session, inputId = "CalButton", label = "Model calibrated", icon = icon("check"))
    CAL_click$valueButton <- 1
102
103
104
105
106
107
    shinyjs::disable("CalButton")
    ## Enable calibration
    # if (input$Period[1L] != input$Period[2L]) {
    #   lapply(getInputs(), shinyjs::enable)
    #   shinyjs::enable("CalButton")
    # }
108
  }, priority = +20)
109

unknown's avatar
unknown committed
110
  
111
  ## Manual calibration
112
  observeEvent({input$Dataset ; input$HydroModel ; input$SnowModel ;
113
    input$X1 ; input$X2 ; input$X3 ; input$X4 ; input$X5 ; input$X6 ; input$C1 ; input$C2 ; 
114
    input$TypeCrit ; input$Period}, {
115
      CAL_click$valueButton <- CAL_click$valueButton - 1
116
117
118
      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"))
119
        shinyjs::enable("CalButton")
120
121
      }
      
122
      ## Enable all inputs except automatic calibration
123
124
125
      if (input$Period[1L] != input$Period[2L]) {
        lapply(getInputs(), shinyjs::enable)
      }
126
    })
127

128
  
129
  
130
  ## --------------- Simulation
131
  
132
  getSim <- reactive({
133
    PARAM <- c(input$X1, input$X2, input$X3, input$X4, input$X5, input$X6)[seq_len(getPrep()$TMGR$NbParam)]
134
    if (input$SnowModel == "CemaNeige") {
135
      PARAM <- c(PARAM, input$C1, input$C2)
136
    }
137
    
138
    ## Simulated flows computation
139
    SIM <- SimGR(PrepGR = getPrep()$PREP, Param = PARAM,
140
                 WupPer = substr(getPrep()$WUPPER, 1, 10),
141
                 SimPer = substr(c(input$Period[1], input$Period[2]), 1, 10),
142
                 verbose = FALSE)
143
    
144
    ## Criteria computation
145
146
    CRIT_opt <- list(Crit    = c("ErrorCrit_NSE", "ErrorCrit_KGE"),
                     Transfo = c("NO", "sqrt", "inv"))
147
148
    CRIT <- lapply(CRIT_opt$Crit, function(iCRIT) {
      Qtransfo <- lapply(CRIT_opt$Transfo, function(iTRSF) {
149
        iInputsCrit <- SIM$OptionsCrit
150
151
        iTRSF <- gsub("NO", "", iTRSF)
        iInputsCrit$transfo <- iTRSF
152
153
154
        iCRIT <- ErrorCrit(InputsCrit = iInputsCrit, OutputsModel = SIM$OutputsModel, FUN_CRIT = get(iCRIT), verbose = FALSE)
        iCRIT <- iCRIT[c("CritName", "CritValue")]
        return(iCRIT)
155
      })
156
      return(Qtransfo)
157
    })
158
    CRIT <- as.data.frame(matrix(na.omit(unlist(CRIT)), ncol = 2, byrow = TRUE), stringsAsFactors = FALSE)
159
    colnames(CRIT) <- c("Criterion", "Value")
160
    rownames(CRIT) <- NULL    
161
162
    CRIT$Value     <- as.numeric(CRIT$Value)
    CRIT$Criterion <- gsub("\\[", " [", CRIT$Criterion)
163
    
164
    ## Recording past simulations
165
    .GlobalEnv$.ShinyGR.hist[[length(.GlobalEnv$.ShinyGR.hist)+1]] <- list(Qsim      = SIM$OutputsModel$Qsim,
166
                                                                           Param     = PARAM,
167
                                                                           TypeModel = SIM$TypeModel,
168
169
                                                                           Crit      = CRIT,
                                                                           Dataset   = input$Dataset)
170
    
171
    .GlobalEnv$.ShinyGR.hist <- .GlobalEnv$.ShinyGR.hist[!(duplicated(sapply(.GlobalEnv$.ShinyGR.hist, function(x) sum(x$Param)), fromLast = TRUE) & 
172
                                                             duplicated(sapply(.GlobalEnv$.ShinyGR.hist, function(x) x$TypeModel), fromLast = TRUE))]
173
    .GlobalEnv$.ShinyGR.hist <- tail(.GlobalEnv$.ShinyGR.hist, n = 2)
174

175
    if (length(.GlobalEnv$.ShinyGR.hist) == 2 &  is.null(names(.GlobalEnv$.ShinyGR.hist[[1]]))) {
176
      .GlobalEnv$.ShinyGR.hist[[1]] <- NULL
177
    }
178
179
180
181
182
    if (length(.GlobalEnv$.ShinyGR.hist) == 2) {
      if (.GlobalEnv$.ShinyGR.hist[[1]]$Dataset != .GlobalEnv$.ShinyGR.hist[[2]]$Dataset) { # reset Qold when new dataset
        .GlobalEnv$.ShinyGR.hist[[1]] <- NULL
      }
    }
183
    if (length(.GlobalEnv$.ShinyGR.hist) == 2 & !is.null(names(.GlobalEnv$.ShinyGR.hist[[1]]))) {
184
185
      isEqualSumQsim   <- sum(.GlobalEnv$.ShinyGR.hist[[1]]$Crit$Value) != sum(.GlobalEnv$.ShinyGR.hist[[2]]$Crit$Value)
      isEqualTypeModel <- .GlobalEnv$.ShinyGR.hist[[1]]$TypeModel == .GlobalEnv$.ShinyGR.hist[[2]]$TypeModel
186
      if (length(.GlobalEnv$.ShinyGR.hist[[1]]$Qsim) != length(.GlobalEnv$.ShinyGR.hist[[2]]$Qsim) |
187
          (isEqualSumQsim & isEqualTypeModel)) {
188
        OBSold <- getPrep()$PREP
189
        OBSold$TypeModel <- .GlobalEnv$.ShinyGR.hist[[1]]$TypeModel
190
        if (.TypeModelGR(OBSold)$CemaNeige & !.TypeModelGR(getPrep()$PREP)$CemaNeige | # present: No CemaNeige ; old: CemaNeige
191
            isEqualSumQsim & isEqualTypeModel) {
192
          if (input$Dataset == "Unnamed watershed") {
193
            ObsDF <- NULL
194
          } else {
195
196
            # ObsDF <- get(input$Dataset)
            ObsDF <- .ShinyGR.args$ObsDF[[input$Dataset]]
197
          }
198
          OBSold <- PrepGR(ObsDF = ObsDF,
199
                          Precip = .ShinyGR.args$Precip, PotEvap = .ShinyGR.args$PotEvap,
200
                          Qobs = .ShinyGR.args$Qobs, TempMean = .ShinyGR.args$TempMean, 
201
202
203
                          ZInputs = .ShinyGR.args$ZInputs[[input$Dataset]],
                          HypsoData = .ShinyGR.args$HypsoData[[input$Dataset]],
                          NLayers = .ShinyGR.args$NLayers[[input$Dataset]],
204
205
                          HydroModel = input$HydroModel,
                          CemaNeige = input$SnowModel == "CemaNeige")
206
        }
207
        SIMold <- SimGR(PrepGR = OBSold,
208
209
210
211
                        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)
212
213
        CRITold <- lapply(CRIT_opt$Crit, function(iCRIT) {
          SIM_transfo <- lapply(CRIT_opt$Transfo, function(iTRSF) {
214
215
            iTRSF <- gsub("NO", "", iTRSF)
            SIMold$OptionsCrit$transfo <- iTRSF
216
217
218
219
220
221
222
223
224
225
            iCRITold <- ErrorCrit(InputsCrit = SIMold$OptionsCrit, OutputsModel = SIMold$OutputsModel, FUN_CRIT = get(iCRIT), verbose = FALSE)
            iCRITold <- iCRITold[c("CritName", "CritValue")]
            return(iCRITold)
          })
        })
        CRITold <- as.data.frame(matrix(na.omit(unlist(CRITold)), ncol = 2, byrow = TRUE), stringsAsFactors = FALSE)
        colnames(CRITold) <- c("Criterion", "Value")
        rownames(CRITold) <- NULL    
        CRITold$Value     <- as.numeric(CRITold$Value)
        CRITold$Criterion <- gsub("\\[", " [", CRITold$Criterion)
226
        
227
        .GlobalEnv$.ShinyGR.hist[[1]]$Crit <- CRITold
228
        .GlobalEnv$.ShinyGR.hist[[1]]$Qsim <- SIMold$OutputsModel$Qsim
229
      }
230
231
    }
    
232
    return(list(PARAM = PARAM, SIM = SIM, SIMold = .GlobalEnv$.ShinyGR.hist, Crit = CRIT))
233
234
    
  })
235
  
236
  
237
238
239
240
  
  ## --------------- Plot
  
  ## Choice
241
242
243
244
245
246
247
  getPlotType <- reactive({
    switch(input$PlotType,
           "Model performance" = 1,
           "Flow time series"  = 2,
           "State variables"   = 3,
           "Model diagram"     = 4)
  })
248
  
249
250
251
  
  ## Models available considering the plot type
  observe({
252
    if (getPlotType() == 4) {
253
      updateSelectInput(session, inputId = "HydroModel", choice = c("GR4J", "GR5J", "GR6J"), selected = input$HydroModel)
254
      updateSelectInput(session, inputId = "SnowModel" , choice = c("None"))
255
    } else {
256
      updateSelectInput(session, inputId = "HydroModel", choice = c("GR4J", "GR5J", "GR6J"), selected = input$HydroModel)
257
      updateSelectInput(session, inputId = "SnowModel" , choice = c("None", "CemaNeige")   , selected = input$SnowModel)
258
    }
259
260
  })
  
unknown's avatar
unknown committed
261
  
262
  ## Plots available considering the model type
263
264
265
266
267
268
269
270
271
272
273
  # observe({
  #   if (input$HydroModel == "GR6J") {
  #     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)
  #   }
  # })
274
275
  
  
276
277
278
279
280
281
282
283
284
285
286
287
288
  # Formated simulation results
  getData <- reactive({
    OutputsModel <- getSim()$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(getSim()$SIM$Qobs[IndPlot]))
    
    if (length(OutputsModel2$DatesR) != 0) {
      data <- data.frame(DatesR  = OutputsModel2$DatesR,
                         precip. = OutputsModel2$Precip,
                         PET     = OutputsModel2$PotEvap,
                         prod.   = OutputsModel2$Prod,
                         rout.   = OutputsModel2$Rout,
289
290
291
                         # exp.    = rep(NA, length(OutputsModel2$DatesR)),
                         # 'exp. (+)'= rep(NA, length(OutputsModel2$DatesR)),
                         # 'exp. (-)'= rep(NA, length(OutputsModel2$DatesR)),
292
293
294
295
                         Qr      = OutputsModel2$QR,
                         Qd      = OutputsModel2$QD,
                         Qsim    = OutputsModel2$Qsim,
                         Qobs    = OutputsModel2$Qobs,
296
297
                         QsimOld = rep(NA, length(OutputsModel2$DatesR)))
                         # QrExp   = rep(NA, length(OutputsModel2$DatesR)))
298
299
300
301
302
      
      if (length(.GlobalEnv$.ShinyGR.hist) == 2 & input$ShowOldQsim == "Yes") {
        data$QsimOld <- .GlobalEnv$.ShinyGR.hist[[1]]$Qsim
      }
      if (input$HydroModel == "GR6J") {
303
304
305
        data$'exp.'    <- NULL
        data$'exp. (+)'<- ifelse(OutputsModel2$Exp >= 0, OutputsModel2$Exp, NA)
        data$'exp. (-)'<- ifelse(OutputsModel2$Exp <  0, OutputsModel2$Exp, NA)
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
        data$QrExp <- OutputsModel2$QRExp
      }
      
      return(list(OutputsModel = OutputsModel2, Tab = data))
    }
  })
  
  
  ## 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$dyPlotSVq_date_window) && getPlotType() == 3) {
      dateWindow <- as.POSIXct(strftime(input$dyPlotSVq_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")) {
      # if (dateWindow[1L] == dateWindow[2L]) {
      #   if (dateWindow[1L] == as.POSIXct(.ShinyGR.args$SimPer[2L], tz = "UTC")) {
      #     updateSliderInput(session, inputId = "Period",
      #                       value = dateWindow - c(1, 0) * .TypeModelGR(input$HydroModel)$TimeLag)
      #   } else {
      #     updateSliderInput(session, inputId = "Period",
      #                       value = dateWindow + c(0, 1) * .TypeModelGR(input$HydroModel)$TimeLag)
      #   }
      # } else {
      if (dateWindow[1L] != dateWindow[2L]) {
        updateSliderInput(session, inputId = "Period",
                          value = dateWindow + .TypeModelGR(input$HydroModel)$TimeLag)
      }
      # }
    }
  }, priority = +100)
  
  
  # observe({
  #   if (getPlotType() == 1) {
  #     if (input$Period[1L] == input$Period[2L]) {
  #       if (input$Period[1L] == as.POSIXct(.ShinyGR.args$SimPer[2L], tz = "UTC")) {
  #         updateSliderInput(session, inputId = "Period",
  #                           value = input$Period - c(1, 0) * .TypeModelGR(input$HydroModel)$TimeLag)
  #       } else {
  #         updateSliderInput(session, inputId = "Period",
  #                           value = input$Period + c(0, 1) * .TypeModelGR(input$HydroModel)$TimeLag)
  #       }
  #     }
  #   }
  # }, priority = +100)
  
  ## Disable all inputs if there is no data
  observe({
    if (input$Period[1L] == input$Period[2L]) {
      inputs <- gsub("Period", "CalButton", getInputs())
      lapply(inputs, shinyjs::disable)
    }
  }, priority = -100)
unknown's avatar
unknown committed
365
  
366
367
368
369
  
  ## Reset period slider responds to dygraphs to mouse clicks
  observeEvent({input$dyPlotTS_click}, {
    updateSliderInput(session, inputId = "Period",
370
                      value = as.POSIXct(.ShinyGR.args$SimPer[[input$Dataset]], tz = "UTC"))
371
372
373
  }, priority = +10)
  observeEvent({input$dyPlotSVs_click}, {
    updateSliderInput(session, inputId = "Period",
374
                      value = as.POSIXct(.ShinyGR.args$SimPer[[input$Dataset]], tz = "UTC"))
375
376
377
  }, priority = +10)
  observeEvent({input$dyPlotSVq_click}, {
    updateSliderInput(session, inputId = "Period",
378
                      value = as.POSIXct(.ShinyGR.args$SimPer[[input$Dataset]], tz = "UTC"))
379
380
381
  }, priority = +10)
  observeEvent({input$dyPlotMDp_click}, {
    updateSliderInput(session, inputId = "Period",
382
                      value = as.POSIXct(.ShinyGR.args$SimPer[[input$Dataset]], tz = "UTC"))
383
384
385
  }, priority = +10)
  observeEvent({input$dyPlotMDe_click}, {
    updateSliderInput(session, inputId = "Period",
386
                      value = as.POSIXct(.ShinyGR.args$SimPer[[input$Dataset]], tz = "UTC"))
387
388
389
  }, priority = +10)
  observeEvent({input$dyPlotMDq_click}, {
    updateSliderInput(session, inputId = "Period",
390
                      value = as.POSIXct(.ShinyGR.args$SimPer[[input$Dataset]], tz = "UTC"))
391
392
  }, priority = +10)
  
393
394

  ## Time window slider
395
396
397
398
399
  observeEvent({input$Dataset}, {
    updateSliderInput(session, inputId = "Period",
                      min = as.POSIXct(.ShinyGR.args$SimPer[[input$Dataset]][1L], tz = "UTC"),
                      max = as.POSIXct(.ShinyGR.args$SimPer[[input$Dataset]][2L], tz = "UTC"),
                      value = as.POSIXct(.ShinyGR.args$SimPer[[input$Dataset]], tz = "UTC"))
400
  })
401
402
403
  
  
  ## Target date slider
404
  eventReactive({input$Dataset}, {
405
406
407
408
409
    updateSliderInput(session, inputId = "Event", label = "Select the target date:",
                      min = as.POSIXct(.ShinyGR.args$SimPer[[input$Dataset]][1L], tz = "UTC") + .TypeModelGR(input$HydroModel)$TimeLag,
                      max = as.POSIXct(.ShinyGR.args$SimPer[[input$Dataset]][2L], tz = "UTC"),
                      value = as.POSIXct(.ShinyGR.args$SimPer[[input$Dataset]][1L], tz = "UTC"), + .TypeModelGR(input$HydroModel)$TimeLag)
  })
unknown's avatar
unknown committed
410
  observe({
411
    updateSliderInput(session, inputId = "Event", label = "Select the target date:",
412
                      min = input$Period[1L] + .TypeModelGR(input$HydroModel)$TimeLag,
unknown's avatar
unknown committed
413
                      max = input$Period[2L])
414
  })
unknown's avatar
unknown committed
415
416
  
  
417
418
419
420
421
422
423
424
  ## 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"
425
      col_fg <- "black"
426
427
428
429
430
431
432
433
      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)))
  })
434
  
435
436
437
  
  ## Plot model performance
  output$stPlotMP <- renderPlot({
438
439
440
    if (length(getSim()$SIM$OutputsModel$DatesR) < 2) {
      return(NULL)
    }
441
    OutputsModel <- getSim()$SIM$OutputsModel
442
    IndPlot <- which(OutputsModel$DatesR >= input$Period[1L] & OutputsModel$DatesR <= input$Period[2L])
443
444
    par(getPlotPar()$par)
    par(cex.axis = 1.2)
445
446
447
    if (input$SnowModel != "CemaNeige") {
      par(oma = c(20, 0, 0, 0))
    }
448
    plot(OutputsModel, Qobs = getSim()$SIM$Qobs, IndPeriod_Plot = IndPlot, cex.lab = 1.2, cex.axis = 1.4, cex.leg = 1.4)
449
  }, bg = "transparent")
450
  
unknown's avatar
unknown committed
451
  
452
  ## Plot flow time series
453
  output$dyPlotTS <- dygraphs::renderDygraph({
454
455
456
    if (length(getSim()$SIM$OutputsModel$DatesR) < 2) {
      return(NULL)
    }
457
458
    if (length(getSim()$SIMold) == 2 & input$ShowOldQsim == "Yes") {
      QsimOld <- getSim()$SIMold[[1]]$Qsim
459
460
461
    } else {
      QsimOld <- NULL
    }
462
    op <- getPlotPar()$par
463
    dg1 <- dyplot(getSim()$SIM, Qsup = QsimOld, Qsup.name = "Qold", RangeSelector = FALSE, LegendShow = "auto",
464
                  col.Q = c(op$fg, "orangered", "grey"), col.Precip = c("#428BCA", "lightblue"))
465
466
    dg1 <- dygraphs::dyOptions(dg1, axisLineColor = op$fg, axisLabelColor = op$fg, retainDateWindow = FALSE)
    dg1 <- dygraphs::dyLegend(dg1, show = "follow", width = 325)
467
468
  })
  
unknown's avatar
unknown committed
469
  
470
  ## Plot state variables stores
471
  output$dyPlotSVs <- dygraphs::renderDygraph({
472
473
474
    if (length(getSim()$SIM$OutputsModel$DatesR) < 2) {
      return(NULL)
    }
475
476
477
478
    # OutputsModel <- getSim()$SIM$OutputsModel
    # data <- data.frame(DatesR = OutputsModel$DatesR,
    #                    prod.  = OutputsModel$Prod,
    #                    rout.  = OutputsModel$Rout)
479

480
481
482
483
484
485
486
487
488
    data <- getData()$Tab[, c("DatesR", "prod.", "rout.", grep("^exp", colnames(getData()$Tab), value = TRUE))]
    data.xts <- xts::xts(data[, -1L], order.by = data$DatesR)

    if (input$HydroModel == "GR6J") {
      colors = c("#00008B", "#008B8B", "#10B510", "#FF0303")
    } else {
      colors = c("#00008B", "#008B8B")
    }
        
489
    op <- getPlotPar()$par
490
    dg2 <- dygraphs::dygraph(data.xts, group = "state_var", ylab = "store [mm]")
491
    dg2 <- dygraphs::dyOptions(dg2, colors = colors,
492
493
494
495
                               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")
496
497
  })
  
unknown's avatar
unknown committed
498
  
499
  ## Plot state variables Q
500
  output$dyPlotSVq <- dygraphs::renderDygraph({
501
502
503
    if (length(getSim()$SIM$OutputsModel$DatesR) < 2) {
      return(NULL)
    }
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
    # OutputsModel <- getSim()$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(getSim()$SIM$Qobs[IndPlot]))
    # 
    # data <- data.frame(DatesR = OutputsModel2$DatesR,
    #                    Qr     = OutputsModel2$QR,
    #                    Qd     = OutputsModel2$QD,
    #                    Qsim   = OutputsModel2$Qsim,
    #                    Qobs   = OutputsModel2$Qobs)
    # if (input$HydroModel == "GR6J") {
    #   data$QrExp <- OutputsModel2$QRExp
    # } else {
    #   data$QrExp <- NA
    # }
519
520
521
522
523
524
525
    
    colSelec <- c("DatesR", "Qr", "Qd", grep("^QrExp", colnames(getData()$Tab), value = TRUE), "Qsim", "Qobs")
    if (length(getSim()$SIMold) == 2 & input$ShowOldQsim == "Yes") {
      colSelec <- c(colSelec, "QsimOld")
    }
    
    data <- getData()$Tab[, colSelec]
526
    data.xts <- xts::xts(data[, -1L], order.by = data$DatesR)
527
    
528
529
530
531
532
533
534
    if (input$HydroModel == "GR6J") {
      names  <- c("Qd", "Qr", "QrExp")
      colors <- c("#FFD700", "#EE6300", "brown")
    } else {
      names  <- c("Qd", "Qr")
      colors <- c("#FFD700", "#EE6300")
    }
535

536
    op <- getPlotPar()$par
537
    dg3 <- dygraphs::dygraph(data.xts, group = "state_var", ylab = paste0("flow [mm/", getPrep()$TMGR$TimeUnit, "]"), main = " ")
538
539
540
    dg3 <- dygraphs::dyOptions(dg3, fillAlpha = 1.0,
                               axisLineColor = op$fg, axisLabelColor = op$fg,
                               titleHeight = 10, retainDateWindow = FALSE)
541
542
    dg3 <- dygraphs::dyStackedRibbonGroup(dg3, name = names,
                                color = colors, strokeBorderColor = "black")
543
544
    dg3 <- dygraphs::dySeries(dg3, name = "Qobs", fillGraph = FALSE, drawPoints = TRUE, color = op$fg)
    dg3 <- dygraphs::dySeries(dg3, name = "Qsim", fillGraph = FALSE, color = "orangered")
545
546
547
    if (length(getSim()$SIMold) == 2 & input$ShowOldQsim == "Yes") {
      dg3 <- dygraphs::dySeries(dg3, name = "QsimOld", label = "Qold", fillGraph = FALSE, color = "grey", strokePattern = "dashed")
    }
548
549
    dg3 <- dygraphs::dyCrosshair(dg3, direction = "vertical")
    dg3 <- dygraphs::dyLegend(dg3, show = "always", width = 325)
550
  })
551
  
unknown's avatar
unknown committed
552
  
553
  ## Plot model diagram precipitation
554
  output$dyPlotMDp <- dygraphs::renderDygraph({
555
556
557
558
559
560
    if (length(getSim()$SIM$OutputsModel$DatesR) < 2) {
      return(NULL)
    }
    data <- data.frame(DatesR  = getSim()$SIM$OutputsModel$DatesR,
                       precip. = getSim()$SIM$OutputsModel$Precip)
    # data <- getData()$Tab[, c("DatesR", "precip.")]
561
    data.xts <- xts::xts(data[, -1L, drop = FALSE], order.by = data$DatesR)
562
    
563
    dg4 <- dygraphs::dygraph(data.xts, group = "mod_diag", ylab = paste0("precip. [mm/", getPrep()$TMGR$TimeUnit, "]"))
564
    dg4 <- dygraphs::dyOptions(dg4, colors = "#428BCA", drawXAxis = FALSE, retainDateWindow = FALSE)
565
    dg4 <- dygraphs::dyBarSeries(dg4, name = "precip.")
566
567
568
569
    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")
570
  })
571
  
unknown's avatar
unknown committed
572
  
573
  ## Plot model diagram ETP
574
  output$dyPlotMDe <- dygraphs::renderDygraph({
575
576
577
    if (length(getSim()$SIM$OutputsModel$DatesR) < 2) {
      return(NULL)
    }
578
579
580
    # data <- data.frame(DatesR = getSim()$SIM$OutputsModel$DatesR,
    #                    PET    = getSim()$SIM$OutputsModel$PotEvap)
    data <- getData()$Tab[, c("DatesR", "PET")]
581
    data.xts <- xts::xts(data[, -1L, drop = FALSE], order.by = data$DatesR)
582
    
583
    op <- getPlotPar()$par
584
    dg5 <- dygraphs::dygraph(data.xts, group = "mod_diag", ylab = paste0("PET [mm/", getPrep()$TMGR$TimeUnit, "]"), main = " ")
585
586
587
588
589
590
591
    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")
592
  })
593
  
unknown's avatar
unknown committed
594
  
595
  ## Plot model diagram flow
596
  output$dyPlotMDq <- dygraphs::renderDygraph({
597
598
599
    if (length(getSim()$SIM$OutputsModel$DatesR) < 2) {
      return(NULL)
    }
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
    # if (length(getSim()$SIMold) == 2 & input$ShowOldQsim == "Yes") {
    #   QsimOld <- getSim()$SIMold[[1]]$Qsim
    # } else {
    #   QsimOld <- NA
    # }
    # OutputsModel <- getSim()$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(getSim()$SIM$Qobs[IndPlot]))
    # OutputsModel2$Qsim <- ifelse(format(OutputsModel2$DatesR, "%Y%m%d") > format(input$Event, "%Y%m%d"), NA, OutputsModel2$Qsim)
    # OutputsModel2$Qold <- ifelse(format(OutputsModel2$DatesR, "%Y%m%d") > format(input$Event, "%Y%m%d"), NA, QsimOld[IndPlot])
    # 
    # data <- data.frame(DatesR  = OutputsModel2$DatesR,
    #                    Qobs    = OutputsModel2$Qobs,
    #                    Qsim    = OutputsModel2$Qsim,
    #                    QsimOld = OutputsModel2$Qold)
    data <- getData()$Tab[, c("DatesR", "Qobs", "Qsim", "QsimOld")]
    data$Qsim    <- ifelse(format(data$DatesR, "%Y%m%d") > format(input$Event, "%Y%m%d"), NA, data$Qsim)
    data$QsimOld <- ifelse(format(data$DatesR, "%Y%m%d") > format(input$Event, "%Y%m%d"), NA, data$QsimOld)
619
    data.xts <- xts::xts(data[, -1L, drop = FALSE], order.by = data$DatesR)
620
621
    
    op <- getPlotPar()$par
622
    dg6 <- dygraphs::dygraph(data.xts, group = "mod_diag", ylab = paste0("flow [mm/", getPrep()$TMGR$TimeUnit, "]"), main = " ")
623
    dg6 <- dygraphs::dyOptions(dg6, colors = c(op$fg, "orangered", "grey"), drawPoints = TRUE,
624
625
626
627
                               axisLineColor = op$fg, axisLabelColor = op$fg,
                               titleHeight = 10, retainDateWindow = FALSE)
    dg6 <- dygraphs::dySeries(dg6, name = "Qsim"   , drawPoints = FALSE)
    dg6 <- dygraphs::dyEvent(dg6, input$Event, color = "orangered")
628
    dg6 <- dygraphs::dySeries(dg6, name = "QsimOld", label = "Qold", drawPoints = FALSE, strokePattern = "dashed")
629
630
    dg6 <- dygraphs::dyLegend(dg6, show = "onmouseover", width = 225)
    dg6 <- dygraphs::dyCrosshair(dg6, direction = "vertical")
631
  })
632
  
unknown's avatar
unknown committed
633
  
634
635
  ## Plot model diagram chart
  output$stPlotMD <- renderPlot({
636
637
638
    if (length(getSim()$SIM$OutputsModel$DatesR) < 2) {
      return(NULL)
    }
639
640
641
642
643
644
    # OutputsModel <- getSim()$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(getSim()$SIM$Qobs[IndPlot]))
    
    OutputsModel2 <- getData()$OutputsModel
645
    
646
    par(getPlotPar()$par)
647
648
649
    .DiagramGR(OutputsModel = OutputsModel2, Param = getSim()$PARAM,
               SimPer = input$Period, EventDate = input$Event,
               HydroModel = input$HydroModel)
650
  }, bg = "transparent")
651
652
  
  
653
  
654
655
  ## --------------- Criteria table
  
656
  output$Criteria <- renderTable({
657

658
    ## Table created in order to choose order the criteria in the table output
659
660
    tabCrit_gauge <- data.frame(Criterion = c("NSE [Q]", "NSE [sqrt(Q)]", "NSE [1/Q]",
                                              "KGE [Q]", "KGE [sqrt(Q)]", "KGE [1/Q]"),
661
                                ID        = 1:6, stringsAsFactors = FALSE)
662
    
663
664
665
    if (length(getSim()$SIMold) == 2 & input$ShowOldQsim == "Yes") {
      tabCrit_old <- getSim()$SIMold[[1]]$Crit$Value
      tabCrit_val <- cbind(getSim()$Crit, tabCrit_old)
666
      colnames(tabCrit_val) <- c(colnames(getSim()$Crit), "Qold")
667
668
      CellColHisto <- '<div style="color: #808080;"><span>9999</span></div>'
    } else {
669
      tabCrit_val <- getSim()$Crit
670
671
    }
    tabCrit_out <- merge(tabCrit_gauge, tabCrit_val, by = "Criterion", all.x = TRUE)
672
673
    tabCrit_out <- tabCrit_out[order(tabCrit_out$ID), ]
    tabCrit_out <- tabCrit_out[, !colnames(tabCrit_out) %in% "ID"]
674
    tabCrit_out[tabCrit_out <= -99.99] <- - 99.99
675
676
    tabCrit_out[, seq_len(ncol(tabCrit_out))[-1]] <- sapply(seq_len(ncol(tabCrit_out))[-1], function(x) sprintf("%7.2f", tabCrit_out[, x]))
    tabCrit_out <- as.data.frame(tabCrit_out)
677
    tabCrit_out[tabCrit_out == " -99.99"] <- "< - 99.99"
678
    colnames(tabCrit_out) <- gsub("Value", "Qsim", colnames(tabCrit_out))
679
680
681
    
    ## Color the cell of the crietaia uses during the calibration
    if (CAL_click$valueButton >= 0) {
682
      CellColCalib <- '<div style="color: #FFFFFF; background-color: #A4C400; border: 5px solid #A4C400; position:relative; top: 0px; left: 5px; padding: 0px; margin: -5px -0px -8px -10px;">
683
<span>9999</span></div>'
684
685
686
      CellColCalib_id  <- which(tabCrit_out$Criterion == input$TypeCrit)
      tabCrit_out[CellColCalib_id, 2] <- gsub("9999", tabCrit_out[CellColCalib_id, 2], CellColCalib)
    }
687
    if (input$ShowOldQsim == "Yes" & length(getSim()$SIMold) > 1) {
688
      tabCrit_out[, "Qold"] <- apply(tabCrit_out[, "Qold", drop = FALSE], 1, function(x) gsub("9999", x, CellColHisto))
689
690
691
    }
    
    return(tabCrit_out)
692
  }, align = c("r"), sanitize.text.function = function(x) x)
693
  
694
  
695
696
697
  
  ## --------------- Download buttons
  
698
  ## Download simulation table
699
700
701
  output$DownloadTab <- downloadHandler(
    filename = function() {
      filename <- "TabSim"
702
      filename <- sprintf("airGR_%s_%s.csv", filename, gsub("(.*)( )(\\d{2})(:)(\\d{2})(:)(\\d{2})", "\\1_\\3h\\5m\\7s", Sys.time()))
703
704
    },
    content = function(file) {
705
706
      PREP <- getPrep()$PREP
      SIM  <- getSim()$SIM
707
708
709
      if (input$SnowModel != "CemaNeige") {
        PrecipSim <- NA
        FracSolid <- NA
710
        TempMean  <- NA
711
      } else {
712
713
        PrecipSol <- rowMeans(as.data.frame(PREP$InputsModel$LayerPrecip) * as.data.frame(PREP$InputsModel$LayerFracSolidPrecip), na.rm = TRUE)
        PrecipSim <- rowMeans(as.data.frame(PREP$InputsModel$LayerPrecip), na.rm = TRUE)
714
715
716
717
        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]
718
        FracSolid <- round(FracSolid, digits = 3)
719
        TempMean  <- rowMeans(as.data.frame(PREP$InputsModel$LayerTempMean), na.rm = TRUE)
720
        TempMean  <- TempMean[SIM$OptionsSimul$IndPeriod_Run]
721
      }
722
723
724
725
726
727
728
729
      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)
730
      colnames(TabSim) <- sprintf("%s [%s]", colnames(TabSim), c("-", rep("mm", 3), "-", "°C", rep("mm", 2)))
731
      colnames(TabSim) <- ifelse(grepl("mm", colnames(TabSim)),
732
                                 gsub("mm", paste0("mm/", .TypeModelGR(PREP)$TimeUnit), colnames(TabSim)),
733
                                 colnames(TabSim))
734
735
736
737
      write.table(TabSim, file = file, row.names = FALSE, sep = ";")
    }
  )
  
738
  
739
  ## Download plots
740
741
742
  output$DownloadPlot <- downloadHandler(
    filename = function() {
      filename <- switch(input$PlotType,
743
744
745
746
                         "Model performance" = "PlotModelPerf",
                         "Flow time series"  = "PlotFlowTimeSeries",
                         "State variables"   = "PlotStateVar",
                         "Model diagram"     = "PlotModelDiag")
747
      filename <- sprintf("airGR_%s_%s.png", filename, gsub("(.*)( )(\\d{2})(:)(\\d{2})(:)(\\d{2})", "\\1_\\3h\\5m\\7s", Sys.time()))
748
749
750
    },
    content = function(file) {
      k <- 1.75
751
752
      ParamTitle <- c("X1", "X2"   , "X3", "X4", "X5", "X6")[seq_len(getPrep()$TMGR$NbParam)]
      ParamUnits <- c("mm", "mm/%s", "mm", "%s",   "", "mm")[seq_len(getPrep()$TMGR$NbParam)]
753
754
      if (input$SnowModel == "CemaNeige") {
        ParamTitle <- c(ParamTitle, "C1", "C2")
755
        ParamUnits <- c(ParamUnits,   "", "mm/°C/%s")
756
      }
757
      ParamTitle <- paste(ParamTitle, paste(getSim()$PARAM, sprintf(ParamUnits, getPrep()$TMGR$TimeUnit)), sep = " = ", collapse = ", ")
758
      ParamTitle <- gsub(" ,", ",", ParamTitle)
759
760
761
762
      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)
763
764
      if (getPlotType() == 1) {
        png(filename = file, width = 1000*k,  height = ifelse(input$SnowModel != "CemaNeige", 700*k, 1100*k), pointsize = 14, res = 150)
765
        par(oma = c(0, 0, 4, 0))
766
        plot(getSim()$SIM)
767
        mtext(text = PngTitle, side = 3, outer = TRUE, cex = 0.8, line = 1.2)
768
769
770
771
        dev.off()
      }
      if (getPlotType() == 2) {
        png(filename = file, width = 1000*k, height = 600*k, pointsize = 14, res = 150)