server.R

# server.R

shinyServer(function(input, output, session) {
  
  
  ## --------------- Data preparation
  
  getPrep <- reactive({
    
    TMGR  <- .TypeModelGR(input$HydroModel)
    PARAM <- c(input$X1, input$X2, input$X3, input$X4, input$X5, input$X6)[seq_len(TMGR$NbParam)]
    
    if (input$SnowModel == "CemaNeige") {
      PARAM <- c(PARAM, input$C1, input$C2)
    }
    
    OBS <- ObsGR(ObsBV = get(input$Dataset), HydroModel = input$HydroModel,
                 CemaNeige = input$SnowModel == "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)
    
    WUPPER <- c(OBS$InputsModel$DatesR[1L], input$Period[1]-.TypeModelGR(OBS)$TimeLag)
    if (WUPPER[2] < WUPPER[1]) {
      WUPPER[2] <- WUPPER[1]
    }
    
    return(list(TMGR = TMGR, OBS = OBS, WUPPER = WUPPER))
    
  })
  
  
  
  ## --------------- Calibration
  
  ## If the user calibrate the model
  CAL_click <- reactiveValues(valueButton = 0)
  
  
  ## Automatic calibration
  observeEvent(input$CalButton, {

    CAL_opt <- list(Crit    = gsub(" .*", "", input$TypeCrit),
                    Transfo = gsub("(\\D{3} \\[)(\\w{0,4})(\\W*Q\\W*\\])", "\\2", input$TypeCrit))
    
    CAL     <- CalGR(ObsGR = getPrep()$OBS, CalCrit = CAL_opt$Crit, transfo = CAL_opt$Transfo,
                     WupPer = substr(getPrep()$WUPPER, 1, 10),
                     CalPer = substr(c(input$Period[1], input$Period[2]), 1, 10), verbose = FALSE)
    PARAM   <- CAL$OutputsCalib$ParamFinalR

    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])
    if (getPrep()$TMGR$NbParam >= 5) {
      updateSliderInput(session, inputId = "X5", value = PARAM[5L])
    }
    if (getPrep()$TMGR$NbParam >= 6) {
      updateSliderInput(session, inputId = "X6", value = PARAM[6L])
    }
    if (input$SnowModel == "CemaNeige") {
      updateSliderInput(session, inputId = "C1", value = PARAM[length(PARAM)-1])
      updateSliderInput(session, inputId = "C2", value = PARAM[length(PARAM)])
    }
    updateActionButton(session, inputId = "CalButton", label = "Model calibrated", icon = icon("check"))
    CAL_click$valueButton <- 1
  })
  
  
  ## Manual calibration
  observeEvent({input$Dataset ; input$HydroModel ; input$SnowModel ;
    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"))
      }
      
    })
  
  
  
  ## --------------- Simulation
  
  getRES <- reactive({

    PARAM <- c(input$X1, input$X2, input$X3, input$X4, input$X5, input$X6)[seq_len(getPrep()$TMGR$NbParam)]
    if (input$SnowModel == "CemaNeige") {
      PARAM <- c(PARAM, input$C1, input$C2)
    }
    
    ## 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) {
        iTRSF <- gsub("NO", "", iTRSF)
        iSIM  <- SimGR(ObsGR = getPrep()$OBS, Param = PARAM,
                       WupPer = substr(getPrep()$WUPPER, 1, 10),
                       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)
        iCRIT <- iCRIT[c("CritName", "CritValue")]
        return(list(SIM = iSIM, CRIT = iCRIT))
      })
      names(SIM_transfo) <- SIM_opt$Transfo
      return(SIM_transfo)
    })
    names(SIM) <- SIM_opt$Crit

    ## Criteria computation
    CRIT <- lapply(SIM, function(iCRIT) {
      lapply(SIM_opt$Transfo, function(iTRSF) {
        iCRIT[[iTRSF]][["CRIT"]]
      })
    })
    CRIT <- as.data.frame(matrix(na.omit(unlist(CRIT)), ncol = 2, byrow = TRUE), stringsAsFactors = FALSE)
    colnames(CRIT) <- c("Criterion", "Value")
    rownames(CRIT) <- NULL    
    CRIT$Value     <- as.numeric(CRIT$Value)
    CRIT$Criterion <- gsub("\\[", " [", CRIT$Criterion)
    
    .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)

    .GlobalEnv$.ShinyGR.hist <- .GlobalEnv$.ShinyGR.hist[!(duplicated(sapply(.GlobalEnv$.ShinyGR.hist, function(x) sum(x$Param)), fromLast = TRUE) & 
                                                           duplicated(sapply(.GlobalEnv$.ShinyGR.hist, function(x) x$TypeModel ), fromLast = TRUE))]
    .GlobalEnv$.ShinyGR.hist <- tail(.GlobalEnv$.ShinyGR.hist, n = 2)
    if (length(.GlobalEnv$.ShinyGR.hist) == 2 &  is.null(names(.GlobalEnv$.ShinyGR.hist[[1]]))) {
      .GlobalEnv$.ShinyGR.hist[[1]] <- NULL
    }
    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)) {
        OBSold <- getPrep()$OBS
        OBSold$TypeModel <- .GlobalEnv$.ShinyGR.hist[[1]]$TypeModel
        if (.TypeModelGR(OBSold)$CemaNeige & !.TypeModelGR(getPrep()$OBS)$CemaNeige) { # present: No CemaNeige ; old: CemaNeige
          OBSold <- ObsGR(ObsBV = get(input$Dataset), HydroModel = .TypeModelGR(OBSold)$NameModel,
                       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)
        }

        SIMold <- SimGR(ObsGR = OBSold,
                        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
      }
    }
    
    return(list(PARAM = PARAM, SIM = SIM$ErrorCrit_KGE$NO$SIM, SIMold = .GlobalEnv$.ShinyGR.hist, Crit = CRIT))
    
  })
  
  
  
  ## --------------- Plot
  
  ## Choice
  getPlotType <- reactive({
    switch(input$PlotType,
           "Model performance" = 1,
           "Flow time series"  = 2,
           "State variables"   = 3,
           "Model diagram"     = 4)
  })
  
  
  ## Models available considering the plot type
  observe({
  if (getPlotType() == 4) {
      updateSelectInput(session, inputId = "HydroModel", choice = c("GR4J", "GR5J"), selected = input$HydroModel)
      updateSelectInput(session, inputId = "SnowModel" , choice = c("None"))
  } else {
      updateSelectInput(session, inputId = "HydroModel", choice = c("GR4J", "GR5J", "GR6J"), selected = input$HydroModel)
      updateSelectInput(session, inputId = "SnowModel" , choice = c("None", "CemaNeige")   , selected = input$SnowModel)
    }
  })
  
  
  ## Plots available considering the model type
  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)
    }
  })
  
  
  ## 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",
                        value = dateWindow + .TypeModelGR(input$HydroModel)$TimeLag)
    }
  })
  
  
  ## 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)
  
  
  ## Target date slider
  observe({
    updateSliderInput(session, inputId = "Event",
                      min = input$Period[1L] + .TypeModelGR(input$HydroModel)$TimeLag,
                      max = input$Period[2L])
  })  
  
  
  ## 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)))
  })

  
  ## Plot model performance
  output$stPlotMP <- renderPlot({
    OutputsModel <- getRES()$SIM$OutputsModel
    IndPlot <- which(OutputsModel$DatesR >= input$Period[1L] & OutputsModel$DatesR <= input$Period[2L])
    par(getPlotPar()$par)
    par(cex.axis = 1.2)
    if (input$SnowModel != "CemaNeige") {
      par(oma = c(20, 0, 0, 0))
    }
    plot(OutputsModel, Qobs = getRES()$SIM$Qobs, IndPeriod_Plot = IndPlot, cex.lab = 1.2, cex.axis = 1.4, cex.leg = 1.4)
  }, bg = "transparent")
  
  
  ## 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)
  })
  
  
  ## 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")
  })
  
  
  ## 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]))

    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
    }
    data.xts <- xts(data[, -1L], order.by = data$DatesR)
    
    op <- getPlotPar()$par
    dg3 <- dygraph(data.xts, group = "state_var", ylab = "flow [mm/d]", main = " ")
    dg3 <- dyOptions(dg3, fillAlpha = 1.0,
                     axisLineColor = op$fg, axisLabelColor = op$fg,
                     titleHeight = 10, retainDateWindow = FALSE)
    dg3 <- dyStackedRibbonGroup(dg3, name = c("Qd", "Qr", "QrExp"),
                                color = c("#FFD700", "#EE6300", "brown"), strokeBorderColor = "black")
    dg3 <- dySeries(dg3, name = "Qobs", fillGraph = FALSE, drawPoints = TRUE, color = op$fg)
    dg3 <- dySeries(dg3, name = "Qsim", fillGraph = FALSE, color = "orangered")
    dg3 <- dyCrosshair(dg3, direction = "vertical")
    dg3 <- dyLegend(dg3, show = "always", width = 325)
  })
  
  
  ## 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)
    dg4 <- dyAxis(dg4, name = "y", valueRange = c(max(data.xts[, "precip."], na.rm = TRUE), -1e-3))
    dg4 <- dyEvent(dg4, input$Event, color = "orangered")
    dg4 <- dyLegend(dg4, show = "onmouseover", width = 225)
    dg4 <- dyCrosshair(dg4, direction = "vertical")
  })
  
  
  ## Plot model diagram ETP
  output$dyPlotMDe <- renderDygraph({
    op <- getPlotPar()$par
    data <- data.frame(DatesR = getRES()$SIM$OutputsModel$DatesR,
                       PET    = getRES()$SIM$OutputsModel$PotEvap)
    data.xts <- xts(data[, -1L, drop = FALSE], order.by = data$DatesR)
    
    dg5 <- dygraph(data.xts, group = "mod_diag", ylab = "PET [mm/d]", main = " ")
    dg5 <- dyOptions(dg5, colors = "#A4C400", drawPoints = TRUE,
                     strokeWidth = 0, pointSize = 2, drawXAxis = FALSE,
                     axisLineColor = op$fg, axisLabelColor = op$fg,
                     titleHeight = 10, retainDateWindow = FALSE)
    dg5 <- dyEvent(dg5, input$Event, color = "orangered")
    dg5 <- dyLegend(dg5, show = "onmouseover", width = 225)
    dg5 <- dyCrosshair(dg5, direction = "vertical")
  })
  
  
  ## Plot model diagram flow
  output$dyPlotMDq <- renderDygraph({
    if (length(.GlobalEnv$.ShinyGR.hist) == 2 & input$ShowOldQsim == "Yes") {
      QsimOld <- getRES()$SIMold[[1]]$Qsim
    } else {
      QsimOld <- NA
    }
    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]))
    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.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, "grey", "orangered"), drawPoints = TRUE,
                     axisLineColor = op$fg, axisLabelColor = op$fg,
                     titleHeight = 10, retainDateWindow = FALSE)
    dg6 <- dySeries(dg6, name = "QsimOld", drawPoints = FALSE, strokePattern = "dashed")
    dg6 <- dySeries(dg6, name = "Qsim"   , drawPoints = FALSE)
    dg6 <- dyEvent(dg6, input$Event, color = "orangered")
    dg6 <- dyLegend(dg6, show = "onmouseover", width = 225)
    dg6 <- dyCrosshair(dg6, direction = "vertical")
  })
  
  
  ## 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,
                              HydroModel = input$HydroModel)
  }, bg = "transparent")
  
  
  
  ## --------------- Criteria table
  
  output$Criteria <- renderTable({
    
    ## Table created in order to choose order the criteria in the table output
    tabCrit_gauge <- data.frame(Criterion = c("NSE [Q]", "NSE [sqrt(Q)]", "NSE [log(Q)]",
                                              "KGE [Q]", "KGE [sqrt(Q)]", "KGE [log(Q)]"),
                                ID        = 1:6, stringsAsFactors = FALSE)
    
    tabCrit_out <- merge(tabCrit_gauge, getRES()$Crit, by = "Criterion", all.x = TRUE)
    tabCrit_out <- tabCrit_out[order(tabCrit_out$ID), ]
    tabCrit_out <- tabCrit_out[, !colnames(tabCrit_out) %in% "ID"]
    
    ## Color the cell of the crietaia uses during the calibration
    if (CAL_click$valueButton >= 0) {
      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;">
<span>9999</span></div>'
      CellCol_id  <- which(tabCrit_out$Criterion == input$TypeCrit)
      tabCrit_out[CellCol_id, 1] <- gsub("9999", tabCrit_out[CellCol_id, 1], CellCol)
    }
    
    return(tabCrit_out)
  }, sanitize.text.function = function(x) x)
  
  
  
  ## --------------- Download buttons
  
  ## simulation table
  output$DownloadTab <- downloadHandler(
    filename = function() {
      filename <- "TabSim"
      filename <- sprintf("airGR_%s_%s.csv", filename, gsub("(.*)( )(\\d{2})(:)(\\d{2})(:)(\\d{2})", "\\1_\\3h\\5m\\7s", Sys.time()))
    },
    content = function(file) {
      OBS <- getPrep()$OBS
      SIM <- getRES()$SIM
      if (input$SnowModel != "CemaNeige") {
        PrecipSim <- NA
        FracSolid <- NA
        TempMean  <- NA
      } 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]
        FracSolid <- round(FracSolid, digits = 3)
        TempMean  <- rowMeans(as.data.frame(OBS$InputsModel$LayerTempMean), na.rm = TRUE)
        TempMean  <- TempMean[SIM$OptionsSimul$IndPeriod_Run]
      }
      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))
      write.table(TabSim, file = file, row.names = FALSE, sep = ";")
    }
  )
  
  ## plots
  output$DownloadPlot <- downloadHandler(
    filename = function() {
      filename <- switch(input$PlotType,
                            "Model performance" = "PlotModelPerf",
                            "Flow time series"  = "PlotFlowTimeSeries",
                            "State variables"   = "PlotStateVar",
                            "Model diagram"     = "PlotModelDiag")
      filename <- sprintf("airGR_%s_%s.png", filename, gsub("(.*)( )(\\d{2})(:)(\\d{2})(:)(\\d{2})", "\\1_\\3h\\5m\\7s", Sys.time()))
    },
    content = function(file) {
      k <- 1.75
      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)
      if (getPlotType() == 1) {
        png(filename = file, width = 1000*k,  height = ifelse(input$SnowModel != "CemaNeige", 700*k, 1100*k), pointsize = 14, res = 150)
        par(oma = c(0, 0, 4, 0))
        plot(getRES()$SIM)
        mtext(text = PngTitle, side = 3, outer = TRUE, cex = 0.8, line = 1.2)
        dev.off()
      }
      if (getPlotType() == 2) {
        png(filename = file, width = 1000*k, height = 600*k, pointsize = 14, res = 150)
        par(oma = c(0, 0, 4, 0))
        plot(getRES()$SIM, which = c( "Precip", "Flows"))
        mtext(text = PngTitle, side = 3, outer = TRUE, cex = 0.8, line = 1.2)
        dev.off()
      }
    }
  )
  
  
})