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)
    
    return(list(TMGR = TMGR, OBS = OBS))
    
  })
  
  
  
  ## --------------- 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 = .ShinyGR.args$WupPer, 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 = .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)
        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)
    
    return(list(PARAM = PARAM, SIM = SIM$ErrorCrit_KGE$NO$SIM, 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)
  })
  
  
  ## 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, # 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({
    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)
    data <- data.frame(DatesR = OutputsModel2$DatesR,
                       Qobs   = OutputsModel2$Qobs,
                       Qsim   = OutputsModel2$Qsim)
    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)
    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)
  })
  
  
  
  ## --------------- 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)
})