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

3
shinyServer(function(input, output, session) {
4
5
  
  
6
7
  ## --------------- Data preparation
  
8
  getPrep <- reactive({
9
    
10
    TMGR  <- .TypeModelGR(input$HydroModel)
11
    PARAM <- c(input$X1, input$X2, input$X3, input$X4, input$X5, input$X6)[seq_len(TMGR$NbParam)]
12
    
13
    if (input$SnowModel == "CemaNeige") {
14
15
      PARAM <- c(PARAM, input$C1, input$C2)
    }
16
    
17
18
    OBS <- ObsGR(ObsBV = get(input$Dataset), HydroModel = input$HydroModel,
                 CemaNeige = input$SnowModel == "CemaNeige",
19
20
                 Precip = .ShinyGR.args$Precip, PotEvap = .ShinyGR.args$PotEvap,
                 Qobs = get(input$Dataset), TempMean = .ShinyGR.args$TempMean, 
21
22
                 ZInputs = .ShinyGR.args$ZInputs, HypsoData = .ShinyGR.args$HypsoData,
                 NLayers = .ShinyGR.args$NLayers)
23
    
24
    return(list(TMGR = TMGR, OBS = OBS))
25
    
26
  })
27
  
28
  
29
30
31
  
  ## --------------- Calibration
  
32
  ## If the user calibrate the model
33
  CAL_click <- reactiveValues(valueButton = 0)
34
  
unknown's avatar
unknown committed
35
  
36
  ## Automatic calibration
37
  observeEvent(input$CalButton, {
38
39
40
    CAL_opt <- list(Crit    = gsub(" .*", "", input$TypeCrit),
                    Transfo = gsub("(\\D{3} \\[)(\\w{0,4})(\\W*Q\\W*\\])", "\\2", input$TypeCrit))
    
41
    CAL     <- CalGR(ObsGR = getPrep()$OBS, CalCrit = CAL_opt$Crit, transfo = CAL_opt$Transfo,
42
                     WupPer = .ShinyGR.args$WupPer, CalPer = substr(c(input$Period[1], input$Period[2]), 1, 10), verbose = FALSE)
43
44
    PARAM   <- CAL$OutputsCalib$ParamFinalR
    
45
46
47
48
    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])
49
    if (getPrep()$TMGR$NbParam >= 5) {
50
51
      updateSliderInput(session, inputId = "X5", value = PARAM[5L])
    }
52
    if (getPrep()$TMGR$NbParam >= 6) {
53
54
      updateSliderInput(session, inputId = "X6", value = PARAM[6L])
    }
55
    if (input$SnowModel == "CemaNeige") {
56
57
58
      updateSliderInput(session, inputId = "C1", value = PARAM[length(PARAM)-1])
      updateSliderInput(session, inputId = "C2", value = PARAM[length(PARAM)])
    }
59
60
    updateActionButton(session, inputId = "CalButton", label = "Model calibrated", icon = icon("check"))
    CAL_click$valueButton <- 1
61
  })
62
  
unknown's avatar
unknown committed
63
  
64
  ## Manual calibration
65
  observeEvent({input$Dataset ; input$HydroModel ; input$SnowModel ;
66
67
68
69
70
71
72
73
74
75
76
77
    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"))
      }
      
    })
  
  
78
  
79
  ## --------------- Simulation
80
81
  
  getRES <- reactive({
82
    
83
    PARAM <- c(input$X1, input$X2, input$X3, input$X4, input$X5, input$X6)[seq_len(getPrep()$TMGR$NbParam)]
84
    if (input$SnowModel == "CemaNeige") {
85
      PARAM <- c(PARAM, input$C1, input$C2)
86
    }
87
    
88
89
90
91
92
    ## 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) {
93
        iTRSF <- gsub("NO", "", iTRSF)
94
95
96
97
        iSIM  <- SimGR(ObsGR = getPrep()$OBS, Param = PARAM, WupPer = .ShinyGR.args$WupPer,
                       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)
98
99
100
        iCRIT <- iCRIT[c("CritName", "CritValue")]
        return(list(SIM = iSIM, CRIT = iCRIT))
      })
101
      names(SIM_transfo) <- SIM_opt$Transfo
102
      return(SIM_transfo)
103
    })
104
    names(SIM) <- SIM_opt$Crit
105
    
106
    ## Criteria computation
107
    CRIT <- lapply(SIM, function(iCRIT) {
108
      lapply(SIM_opt$Transfo, function(iTRSF) {
109
110
        iCRIT[[iTRSF]][["CRIT"]]
      })
111
    })
112
    CRIT <- as.data.frame(matrix(na.omit(unlist(CRIT)), ncol = 2, byrow = TRUE), stringsAsFactors = FALSE)
113
    colnames(CRIT) <- c("Criterion", "Value")
114
    rownames(CRIT) <- NULL    
115
116
    CRIT$Value     <- as.numeric(CRIT$Value)
    CRIT$Criterion <- gsub("\\[", " [", CRIT$Criterion)
117
    
118
    return(list(PARAM = PARAM, SIM = SIM$ErrorCrit_KGE$NO$SIM, Crit = CRIT))
119
120
    
  })
121
  
122
  
123
124
125
126
  
  ## --------------- Plot
  
  ## Choice
127
128
129
130
131
132
133
  getPlotType <- reactive({
    switch(input$PlotType,
           "Model performance" = 1,
           "Flow time series"  = 2,
           "State variables"   = 3,
           "Model diagram"     = 4)
  })
134
  
135
136
137
138
  
  ## Models available considering the plot type
  observe({
  if (getPlotType() == 4) {
139
140
      updateSelectInput(session, inputId = "HydroModel", choice = c("GR4J", "GR5J"), selected = input$HydroModel)
      updateSelectInput(session, inputId = "SnowModel", choice = c("None"))
141
  } else {
142
143
      updateSelectInput(session, inputId = "HydroModel", choice = c("GR4J", "GR5J", "GR6J"), selected = input$HydroModel)
      updateSelectInput(session, inputId = "SnowModel", choice = c("None", "CemaNeige")   , selected = input$SnowModel)
144
    }
145
146
  })
  
unknown's avatar
unknown committed
147
  
148
149
  ## Plots available considering the model type
  observe({
150
    if (input$HydroModel == "GR6J") {
151
152
153
154
155
156
157
      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)
158
    }
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
  })
  
  
  ## 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",
175
                        value = dateWindow + .TypeModelGR(input$HydroModel)$TimeLag)
176
    }
unknown's avatar
unknown committed
177
178
  })
  
179
180
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
  
  ## 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
207
208
209
  ## Target date slider
  observe({
    updateSliderInput(session, inputId = "Event",
210
                      min = input$Period[1L] + .TypeModelGR(input$HydroModel)$TimeLag,
unknown's avatar
unknown committed
211
212
213
214
                      max = input$Period[2L])
  })  
  
  
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
  ## 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)))
  })
unknown's avatar
unknown committed
232

233
234
235
236
  
  ## Plot model performance
  output$stPlotMP <- renderPlot({
    OutputsModel <- getRES()$SIM$OutputsModel
237
    IndPlot <- which(OutputsModel$DatesR >= input$Period[1L] & OutputsModel$DatesR <= input$Period[2L])
238
239
    par(getPlotPar()$par)
    par(cex.axis = 1.2)
240
241
242
    if (input$SnowModel != "CemaNeige") {
      par(oma = c(20, 0, 0, 0))
    }
243
    plot(OutputsModel, Qobs = getRES()$SIM$Qobs, IndPeriod_Plot = IndPlot, cex.lab = 1.2, cex.axis = 1.4, cex.leg = 1.4)
244
245
  })
  
unknown's avatar
unknown committed
246
  
247
248
249
250
251
252
253
  ## Plot flow time series
  output$dyPlotTS <- renderDygraph({
    op <- getPlotPar()$par
    dg1 <- dyplot(getRES()$SIM, RangeSelector = FALSE, LegendShow = "auto",
                  col.Q = c(op$fg, "orangered"), col.Precip = "#428BCA")
    dg1 <- dyOptions(dg1, axisLineColor = op$fg, axisLabelColor = op$fg, retainDateWindow = FALSE)
    dg1 <- dyLegend(dg1, show = "follow", width = 325)
254
255
  })
  
unknown's avatar
unknown committed
256
  
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
  ## Plot state variables stores
  output$dyPlotSVs <- renderDygraph({
    OutputsModel <- getRES()$SIM$OutputsModel
    data <- data.frame(DatesR = OutputsModel$DatesR,
                       prod.  = OutputsModel$Prod,
                       rout.  = OutputsModel$Rout)
    data.xts <- xts(data[, -1L], order.by = data$DatesR)
    
    op <- getPlotPar()$par
    dg2 <- dygraph(data.xts, group = "state_var", ylab = "store [mm]")
    dg2 <- dyOptions(dg2, colors = c("#00008B", "#008B8B"),
                     fillGraph = TRUE, fillAlpha = 0.3,
                     drawXAxis = FALSE, axisLineColor = op$fg, axisLabelColor = op$fg, retainDateWindow = FALSE)
    dg2 <- dyLegend(dg2, show = "always", width = 325)
    dg2 <- dyCrosshair(dg2, direction = "vertical")
272
273
  })
  
unknown's avatar
unknown committed
274
  
275
276
277
278
279
280
  ## Plot state variables Q
  output$dyPlotSVq <- renderDygraph({
    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]))
281
282
283
284
285
286

    data <- data.frame(DatesR = OutputsModel2$DatesR,
                       Qr     = OutputsModel2$QR,
                       Qd     = OutputsModel2$QD, # Qd
                       Qsim   = OutputsModel2$Qsim,
                       Qobs   = OutputsModel2$Qobs)
287
    if (input$HydroModel == "GR6J") {
288
      data$QrExp <- OutputsModel2$QRExp
289
    } else {
290
      data$QrExp <- NA
291
    }
292
293
294
295
    data.xts <- xts(data[, -1L], order.by = data$DatesR)
    
    op <- getPlotPar()$par
    dg3 <- dygraph(data.xts, group = "state_var", ylab = "flow [mm/d]", main = " ")
296
    dg3 <- dyOptions(dg3, fillAlpha = 1.0,
297
298
                     axisLineColor = op$fg, axisLabelColor = op$fg,
                     titleHeight = 10, retainDateWindow = FALSE)
299
    dg3 <- dyStackedRibbonGroup(dg3, name = c("Qd", "Qr", "QrExp"),
300
                                color = c("#FFD700", "#EE6300", "brown"), strokeBorderColor = "black")
301
302
    dg3 <- dySeries(dg3, name = "Qobs", fillGraph = FALSE, drawPoints = TRUE, color = op$fg)
    dg3 <- dySeries(dg3, name = "Qsim", fillGraph = FALSE, color = "orangered")
303
    dg3 <- dyCrosshair(dg3, direction = "vertical")
304
    dg3 <- dyLegend(dg3, show = "always", width = 325)
305
  })
306
  
unknown's avatar
unknown committed
307
  
308
309
310
311
312
313
314
315
316
317
  ## Plot model diagram precipitation
  output$dyPlotMDp <- renderDygraph({
    barChartPrecip <- scan(file = system.file("plugins/barChartPrecip.js", package = "airGRteaching"),
                           what = "character", quiet = TRUE)
    data <- data.frame(DatesR  = getRES()$SIM$OutputsModel$DatesR,
                       precip. = getRES()$SIM$OutputsModel$Precip)
    data.xts <- xts(data[, -1L, drop = FALSE], order.by = data$DatesR)
    
    dg4 <- dygraph(data.xts, group = "mod_diag", ylab = "precip. [mm/d]")
    dg4 <- dyOptions(dg4, colors = "#428BCA", drawXAxis = FALSE, plotter = barChartPrecip, retainDateWindow = FALSE)
318
    dg4 <- dyAxis(dg4, name = "y", valueRange = c(max(data.xts[, "precip."], na.rm = TRUE), -1e-3))
unknown's avatar
unknown committed
319
    dg4 <- dyEvent(dg4, input$Event, color = "orangered")
320
321
    dg4 <- dyLegend(dg4, show = "onmouseover", width = 225)
    dg4 <- dyCrosshair(dg4, direction = "vertical")
322
  })
323
  
unknown's avatar
unknown committed
324
  
325
326
327
328
  ## Plot model diagram ETP
  output$dyPlotMDe <- renderDygraph({
    op <- getPlotPar()$par
    data <- data.frame(DatesR = getRES()$SIM$OutputsModel$DatesR,
329
                       PET    = getRES()$SIM$OutputsModel$PotEvap)
330
331
    data.xts <- xts(data[, -1L, drop = FALSE], order.by = data$DatesR)
    
332
    dg5 <- dygraph(data.xts, group = "mod_diag", ylab = "PET [mm/d]", main = " ")
333
334
335
336
    dg5 <- dyOptions(dg5, colors = "#A4C400", drawPoints = TRUE,
                     strokeWidth = 0, pointSize = 2, drawXAxis = FALSE,
                     axisLineColor = op$fg, axisLabelColor = op$fg,
                     titleHeight = 10, retainDateWindow = FALSE)
unknown's avatar
unknown committed
337
    dg5 <- dyEvent(dg5, input$Event, color = "orangered")
338
339
    dg5 <- dyLegend(dg5, show = "onmouseover", width = 225)
    dg5 <- dyCrosshair(dg5, direction = "vertical")
340
  })
341
  
unknown's avatar
unknown committed
342
  
343
344
345
346
347
348
  ## Plot model diagram flow
  output$dyPlotMDq <- renderDygraph({
    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
349
    OutputsModel2$Qsim <- ifelse(format(OutputsModel2$DatesR, "%Y%m%d") > format(input$Event, "%Y%m%d"), NA, OutputsModel2$Qsim)
350
    data <- data.frame(DatesR = OutputsModel2$DatesR,
351
352
                       Qobs   = OutputsModel2$Qobs,
                       Qsim   = OutputsModel2$Qsim)
353
354
355
356
357
358
359
    data.xts <- xts(data[, -1L, drop = FALSE], order.by = data$DatesR)
    
    op <- getPlotPar()$par
    dg6 <- dygraph(data.xts, group = "mod_diag", ylab = "flow [mm/d]", main = " ")
    dg6 <- dyOptions(dg6, colors = c(op$fg, "orangered"), drawPoints = TRUE,
                     axisLineColor = op$fg, axisLabelColor = op$fg,
                     titleHeight = 10, retainDateWindow = FALSE)
360
    dg6 <- dySeries(dg6, name = "Qsim", drawPoints = FALSE)
unknown's avatar
unknown committed
361
    dg6 <- dyEvent(dg6, input$Event, color = "orangered")
362
363
    dg6 <- dyLegend(dg6, show = "onmouseover", width = 225)
    dg6 <- dyCrosshair(dg6, direction = "vertical")
364
  })
365
  
unknown's avatar
unknown committed
366
  
367
368
369
370
371
372
373
374
375
376
  ## 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]))
    
    par(getPlotPar()$par)
    airGRteaching:::DiagramGR(OutputsModel = OutputsModel2, Param = getRES()$PARAM,
                              SimPer = input$Period, EventDate = input$Event,
377
                              HydroModel = input$HydroModel)
378
  })
379
380
  
  
381
  
382
383
  ## --------------- Criteria table
  
384
  output$Criteria <- renderTable({
385
386
    
    ## Table created in order to choose order the criteria in the table output
387
    tabCrit_gauge <- data.frame(Criterion = c("NSE [Q]", "NSE [sqrt(Q)]", "NSE [log(Q)]",
388
                                             "KGE [Q]", "KGE [sqrt(Q)]", "KGE [log(Q)]"),
389
                                ID        = 1:6, stringsAsFactors = FALSE)
390
    
391
    tabCrit_out <- merge(tabCrit_gauge, getRES()$Crit, by = "Criterion", all.x = TRUE)
392
393
    tabCrit_out <- tabCrit_out[order(tabCrit_out$ID), ]
    tabCrit_out <- tabCrit_out[, !colnames(tabCrit_out) %in% "ID"]
394
395
396
    
    ## Color the cell of the crietaia uses during the calibration
    if (CAL_click$valueButton >= 0) {
397
      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;">
398
<span>9999</span></div>'
399
      CellCol_id  <- which(tabCrit_out$Criterion == input$TypeCrit)
400
401
402
403
404
      tabCrit_out[CellCol_id, 1] <- gsub("9999", tabCrit_out[CellCol_id, 1], CellCol)
    }
    
    return(tabCrit_out)
  }, sanitize.text.function = function(x) x)
405
406
})