Commit 46afbd2a authored by unknown's avatar unknown
Browse files

v0.1.3.0 its is now possible to choose the objective function to calibrate the model #4670

parent b116f036
Package: airGRteaching
Type: Package
Title: Tools to Simplify the Use of the airGR Hydrological Package for Education (Including a Shiny Application)
Version: 0.1.2.52
Date: 2017-06-02
Version: 0.1.3.0
Date: 2017-06-13
Authors@R: c(person("Olivier", "Delaigue", role = c("aut", "cre"), email = "airGR@irstea.fr"), person("Laurent", "Coron", role = c("aut")), person("Pierre", "Brigode", role = c("aut")), person("Guillaume", "Thirel", role = c("ctb")))
Depends: airGR (>= 1.0.5.22)
Imports: xts, dygraphs, shiny, plotrix
......
......@@ -3,9 +3,13 @@
shinyServer(function(input, output, session) {
## --------------- Data preparation
getPrep <- reactive({
TMGR <- .TypeModelGR(input$TypeModel)
PARAM <- c(input$X1, input$X2, input$X3, input$X4, input$X5, input$X6)[seq_len(TMGR$NbParam)]
if (input$CemaNeige == "CemaNeige") {
PARAM <- c(PARAM, input$C1, input$C2)
}
......@@ -15,15 +19,30 @@ shinyServer(function(input, output, session) {
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
##
CAL_click <- reactiveValues(valueButton = 0)
## Automatic calibration
observeEvent(input$CalButton, {
TMGR <- getPrep()$TMGR
OBS <- getPrep()$OBS
CAL <- CalGR(ObsGR = OBS, WupPer = .ShinyGR.args$WupPer, CalPer = substr(c(input$Period[1], input$Period[2]), 1, 10), verbose = FALSE)
PARAM <- CAL$OutputsCalib$ParamFinalR
TMGR <- getPrep()$TMGR
OBS <- getPrep()$OBS
CAL_opt <- list(Crit = gsub(" .*", "", input$TypeCrit),
Transfo = gsub("(\\D{3} \\[)(\\w{0,4})(\\W*Q\\W*\\])", "\\2", input$TypeCrit))
CAL <- CalGR(ObsGR = 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])
......@@ -38,16 +57,30 @@ shinyServer(function(input, output, session) {
updateSliderInput(session, inputId = "C1", value = PARAM[length(PARAM)-1])
updateSliderInput(session, inputId = "C2", value = PARAM[length(PARAM)])
}
updateActionButton(session, inputId = "CalButton", label = "Model calibrated")
updateActionButton(session, inputId = "CalButton", label = "Model calibrated", icon = icon("check"))
CAL_click$valueButton <- 1
})
observeEvent({input$X1 ; input$X2 ; input$X3 ; input$X4 ; input$X5 ; input$X6}, {
updateActionButton(session, inputId = "CalButton", label = "Run")
## Manual calibration
observeEvent({input$Dataset ; input$TypeModel ; input$CemaNeige ;
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({
TMGR <- getPrep()$TMGR
OBS <- getPrep()$OBS
PARAM <- c(input$X1, input$X2, input$X3, input$X4, input$X5, input$X6)[seq_len(TMGR$NbParam)]
......@@ -55,41 +88,43 @@ shinyServer(function(input, output, session) {
PARAM <- c(PARAM, input$C1, input$C2)
}
list_CRIT <- c("ErrorCrit_NSE", "ErrorCrit_KGE")
list_TRSF <- c("NO", "log", "sqrt")
SIM <- lapply(list_CRIT, function(iCRIT) {
SIM_transfo <- lapply(list_TRSF, function(iTRSF) {
## 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 = 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) <- list_TRSF
names(SIM_transfo) <- SIM_opt$Transfo
return(SIM_transfo)
})
names(SIM) <- list_CRIT
names(SIM) <- SIM_opt$Crit
## Criteria computation
CRIT <- lapply(SIM, function(iCRIT) {
lapply(list_TRSF, function(iTRSF) {
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("Criteria", "Value")
rownames(CRIT) <- NULL
CRIT$Value <- as.numeric(CRIT$Value)
CRIT$Criteria <- gsub("\\[", " [", CRIT$Criteria)
# suppressWarnings(rm(list = ls()[!ls() %in% c("PARAM", "SIM", "CRIT")]))
# gc()
return(list(PARAM = PARAM, SIM = SIM$ErrorCrit_KGE$NO$SIM, Crit = CRIT))
})
## --------------- Plot
## Choice
getPlotType <- reactive({
switch(input$PlotType,
"Model performance" = 1,
......@@ -98,7 +133,7 @@ shinyServer(function(input, output, session) {
"Model diagram" = 4)
})
## Plotting
output$myPlot <- renderPlot({
OutputsModel <- getRES()$SIM$OutputsModel
......@@ -111,18 +146,19 @@ shinyServer(function(input, output, session) {
par(fg = "#2C3E50", col.axis = "#2C3E50", col.lab = "#2C3E50")
}
## Target date slider
observe({
updateSliderInput(session, inputId = "Event",
min = input$Period[1L]+.TypeModelGR(input$TypeModel)$TimeLag,
max = input$Period[2L])
})
## Add Qobs to outputs model
IndPlot <- which(OutputsModel$DatesR >= input$Period[1L] & OutputsModel$DatesR <= input$Period[2L])
aaa <- sapply(OutputsModel[ seq_len(which(names(OutputsModel) == "Qsim"))], function(x) x[IndPlot])
aaa <- c(aaa, Qobs = list(getRES()$SIM$Qobs[IndPlot]))
OutputsModel2 <- sapply(OutputsModel[seq_len(which(names(OutputsModel) == "Qsim"))], function(x) x[IndPlot])
OutputsModel2 <- c(OutputsModel2, Qobs = list(getRES()$SIM$Qobs[IndPlot]))
## Models available considering the plot type
if (getPlotType() == 4) {
observe({
updateSelectInput(session, inputId = "TypeModel", choice = c("GR4J"))
......@@ -135,46 +171,67 @@ shinyServer(function(input, output, session) {
})
}
## Flow time series
if (getPlotType() == 2) {
par(mfrow = c(1, 1), par(oma = c(7, 2, 3, 0)))
plot(OutputsModel, Qobs = getRES()$SIM$Qobs, IndPeriod_Plot = IndPlot, which = "Flows", cex.lab = 1.4, cex.axis = 1.5, cex.leg = 1.5)
}
## Model performance
if (getPlotType() == 1) {
par(cex.axis = 1.2)
plot(OutputsModel, Qobs = getRES()$SIM$Qobs, IndPeriod_Plot = IndPlot, cex.lab = 1.2, cex.axis = 1.4, cex.leg = 1.4)
}
## Model diagram
if (getPlotType() == 4) {
airGRteaching:::DiagramGR4J(OutputsModel = aaa, Param = getRES()$PARAM, SimPer = input$Period, EventDate = input$Event)
airGRteaching:::DiagramGR4J(OutputsModel = OutputsModel2, Param = getRES()$PARAM, SimPer = input$Period, EventDate = input$Event)
}
## State variables
if (getPlotType() == 3) {
par(mfrow = c(2, 1), cex.axis = 1.2, cex.lab = 1.1, cex.lab = 1.1)
par(mar = c(2.5, 4.0, 2.5, 2.0), xaxt = "n")
plot(range(aaa$Dates), range(aaa$Prod, aaa$Rout), type = "n", xlab = "", ylab = "store [mm]")
polygon(c(aaa$Dates, rev(range(aaa$Dates))), c(aaa$Prod, rep(min(aaa$Prod, aaa$Rout), 2)), border = "darkblue", col = adjustcolor("darkblue", alpha.f = 0.30))
polygon(c(aaa$Dates, rev(range(aaa$Dates))), c(aaa$Rout, rep(min(aaa$Prod, aaa$Rout), 2)), border = "cyan4" , col = adjustcolor("cyan4" , alpha.f = 0.30))
plot(range(OutputsModel2$Dates), range(OutputsModel2$Prod, OutputsModel2$Rout), type = "n", xlab = "", ylab = "store [mm]")
polygon(c(OutputsModel2$Dates, rev(range(OutputsModel2$Dates))), c(OutputsModel2$Prod, rep(min(OutputsModel2$Prod, OutputsModel2$Rout), 2)), border = "darkblue", col = adjustcolor("darkblue", alpha.f = 0.30))
polygon(c(OutputsModel2$Dates, rev(range(OutputsModel2$Dates))), c(OutputsModel2$Rout, rep(min(OutputsModel2$Prod, OutputsModel2$Rout), 2)), border = "cyan4" , col = adjustcolor("cyan4" , alpha.f = 0.30))
legend("topright", bty = "n", legend = c("prod.", "rout."), fill = adjustcolor(c("darkblue", "cyan4"), alpha.f = 0.30))
par(mar = c(5.0, 4.0, 0.0, 2.0), xaxt = "s")
plot(aaa$DatesR, aaa$Qobs, type = "n", xlab = "", ylab = "flow [mm/d]")
polygon(c(aaa$Dates, rev(range(aaa$Dates))), c(aaa$QR+aaa$QD, rep(0, 2)), col = adjustcolor("orange2", alpha.f = 0.75), border = NA)
polygon(c(aaa$Dates, rev(range(aaa$Dates))), c(aaa$QR , rep(0, 2)), col = adjustcolor("gold" , alpha.f = 0.70), border = NA)
lines(aaa$DatesR, aaa$Qsim, lwd = 1, col = "orangered")
lines(aaa$DatesR, aaa$Qobs, type = "l", lwd = 1, col = par("fg"))
plot(OutputsModel2$DatesR, OutputsModel2$Qobs, type = "n", xlab = "", ylab = "flow [mm/d]")
polygon(c(OutputsModel2$Dates, rev(range(OutputsModel2$Dates))), c(OutputsModel2$QR+OutputsModel2$QD, rep(0, 2)), col = adjustcolor("orange2", alpha.f = 0.75), border = NA)
polygon(c(OutputsModel2$Dates, rev(range(OutputsModel2$Dates))), c(OutputsModel2$QR , rep(0, 2)), col = adjustcolor("gold" , alpha.f = 0.70), border = NA)
lines(OutputsModel2$DatesR, OutputsModel2$Qsim, lwd = 1, col = "orangered")
lines(OutputsModel2$DatesR, OutputsModel2$Qobs, type = "l", lwd = 1, col = par("fg"))
legend("topright", bty = "n", legend = c("obs.", "sim.", "Qd", "Qr"), col = c(par("fg"), "orangered", "orange2", "gold"), lwd = c(2, 2, NA, NA), pch = c(NA, NA, 15, 15))
}
})
## --------------- Criteria table
output$Criteria <- renderTable({
tabCrit_gauge <- data.frame(Criteria = c("NSE [Q]", "NSE [log(Q)]", "NSE [sqrt(Q)]",
"KGE [Q]", "KGE [log(Q)]", "KGE [sqrt(Q)]"),
ID = 1:6)
## Table created in order to choose order the criteria in the table output
tabCrit_gauge <- data.frame(Criteria = 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 = "Criteria", 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="width: 087%; height: 99%; z-index: 0; color: #FFFFFF; background-color: #A4C400; position:absolute; top: 0px; left: 7px; padding:5px;">
<span>9999</span></div>'
CellCol_id <- which(tabCrit_out$Criteria == input$TypeCrit)
tabCrit_out[CellCol_id, 1] <- gsub("9999", tabCrit_out[CellCol_id, 1], CellCol)
}
return(tabCrit_out)
}, sanitize.text.function = function(x) x)
# gc()
})
......@@ -17,7 +17,7 @@ navbarPage(title = div("airGRteaching",
Flatly = "bootstrap.min_Flatly.css",
United = "bootstrap.min_United.css",
Yeti = "bootstrap.min_Yeti.css"
),
),
tabPanel(title = "Interface",
......@@ -30,13 +30,11 @@ navbarPage(title = div("airGRteaching",
sidebarPanel(width = 3,
h4("Choose a dataset:"),
fluidRow(
column(12, selectInput("Dataset", label = NULL, choices = .ShinyGR.args$ObsBV))
),
h4("Choose a model:"),
fluidRow(
column(6, selectInput("TypeModel", label = "Hydrological model", choices = c("GR4J", "GR5J", "GR6J"))),
column(6, selectInput("CemaNeige", label = "Snow module" , choices = c("None", "CemaNeige")))
......@@ -44,7 +42,6 @@ navbarPage(title = div("airGRteaching",
h4("Parameters values:"),
conditionalPanel(condition = "input.TypeModel == 'GR4J' || input.TypeModel =='GR5J' || input.TypeModel =='GR6J'",
sliderInput("X1", label = "X1 (production store capacity)",
post = " [mm]",
......@@ -70,7 +67,6 @@ navbarPage(title = div("airGRteaching",
max = 10,
step = 0.1,
value = .ShinyGR.args$Param[4L])),
conditionalPanel(condition = "input.TypeModel == 'GR5J' || input.TypeModel =='GR6J'",
sliderInput("X5", label = "X5 (intercatchment exchange threshold)",
post = " [-]",
......@@ -78,7 +74,6 @@ navbarPage(title = div("airGRteaching",
max = 3,
step = 0.05,
value = .ShinyGR.args$Param[5L])),
conditionalPanel(condition = "input.TypeModel == 'GR6J'",
sliderInput("X6", label = "X6 (coeff. for emptying exponential store)",
post = " [mm]",
......@@ -86,7 +81,6 @@ navbarPage(title = div("airGRteaching",
max = 20,
step = 0.5,
value = .ShinyGR.args$Param[6L])),
conditionalPanel(condition = "input.CemaNeige == 'CemaNeige'",
sliderInput("C1", label = "C1 (weighting coeff. for snow pack thermal state)",
post = " [-]",
......@@ -94,7 +88,6 @@ navbarPage(title = div("airGRteaching",
max = 1,
step = 0.01,
value = .ShinyGR.args$Param[length(.ShinyGR.args$Param)-1]),
sliderInput("C2", label = "C2 (degree-day melt coefficient)",
post = " [mm/°C/d]",
min = 0,
......@@ -102,14 +95,17 @@ navbarPage(title = div("airGRteaching",
step = 0.5,
value = .ShinyGR.args$Param[length(.ShinyGR.args$Param)])),
h4("Automatic calibration (with KGE [Q]):"),
actionButton("CalButton", label = "Run automatic calibration", width = "100%",
icon = icon("refresh"),
style = ifelse(.GlobalEnv$.ShinyGR.args$theme != "Cerulean",
"color: #fff; background-color: #A4C400; border-color: #A4C400",
""))
h4("Automatic calibration:"),
fluidRow(
column(6, selectInput("TypeCrit", label = "Objective function",
choices = c("NSE [Q]", "NSE [sqrt(Q)]", "NSE [log(Q)]",
"KGE [Q]", "KGE [sqrt(Q)]", "KGE [log(Q)]"))),
column(6, actionButton("CalButton", label = "Run", width = "100%",
icon = icon("refresh"),
style = ifelse(.GlobalEnv$.ShinyGR.args$theme != "Cerulean",
"color: #fff; background-color: #A4C400; border-color: #A4C400; margin-top: 25px; padding:6px;",
"")))
)
),
......@@ -120,50 +116,16 @@ navbarPage(title = div("airGRteaching",
selectInput("PlotType", label = "Choose a plot:",
choices = c("Model performance", "Flow time series", "State variables", "Model diagram"))#, "user-defined"))
),
# conditionalPanel(condition = "input.PlotType != 'Model diagram'",
column(4, offset = 1,
sliderInput("Period", label = "Select the time window:",
min = as.POSIXct(.ShinyGR.args$SimPer[1L], tz = "UTC"),
max = as.POSIXct(.ShinyGR.args$SimPer[2L], tz = "UTC"),
value = as.POSIXct(.ShinyGR.args$SimPer, tz = "UTC"),
timeFormat = "%F",
timezone = "+0000",
animate = FALSE)
),
# ),
# conditionalPanel(condition = "input.PlotType == 'Model diagram'",
# column(4, offset = 1,
# sliderInput("Period", label = "Select the time window:",
# min = as.POSIXct(.ShinyGR.args$SimPer[1L], tz = "UTC"),
# max = as.POSIXct(.ShinyGR.args$SimPer[2L], tz = "UTC"),
# value = as.POSIXct(.ShinyGR.args$SimPer, tz = "UTC"),
# timeFormat = "%F",
# timezone = "+0000",
# animate = FALSE)
# )
# ),
# conditionalPanel(condition = "input.PlotType == 'Model diagram' & input.TypeModel =='GR4H'",
# column(4, offset = 0,
# sliderInput("Event", label = "Select the target date:",
# min = as.POSIXct(.ShinyGR.args$SimPer[1L], tz = "UTC"),
# max = as.POSIXct(.ShinyGR.args$SimPer[2L], tz = "UTC"),
# value = as.POSIXct(.ShinyGR.args$SimPer[1L], tz = "UTC"),
# timeFormat = "%F",
# timezone = "+0000",
# animate = animationOptions(interval = 500),
# step = 3600)
# )
# ),
column(4, offset = 1,
sliderInput("Period", label = "Select the time window:",
min = as.POSIXct(.ShinyGR.args$SimPer[1L], tz = "UTC"),
max = as.POSIXct(.ShinyGR.args$SimPer[2L], tz = "UTC"),
value = as.POSIXct(.ShinyGR.args$SimPer, tz = "UTC"),
timeFormat = "%F",
timezone = "+0000",
animate = FALSE)
),
conditionalPanel(condition = "input.PlotType == 'Model diagram' & (input.TypeModel =='GR4J' || input.TypeModel =='GR5J' || input.TypeModel =='GR6J')",
# column(4, offset = 1,
# sliderInput("Period", label = "Select the time window:",
# min = as.POSIXct(.ShinyGR.args$SimPer[1L], tz = "UTC"),
# max = as.POSIXct(.ShinyGR.args$SimPer[2L], tz = "UTC"),
# value = as.POSIXct(.ShinyGR.args$SimPer, tz = "UTC"),
# timeFormat = "%F",
# timezone = "+0000",
# animate = FALSE)
# ),
column(4, offset = 0,
sliderInput("Event", label = "Select the target date:",
min = as.POSIXct(.ShinyGR.args$SimPer[1L], tz = "UTC"),
......@@ -174,46 +136,22 @@ navbarPage(title = div("airGRteaching",
animate = animationOptions(interval = 500),
step = 3600 * 24)
)
# ),
# conditionalPanel(condition = "input.PlotType == 'Model diagram' & input.TypeModel =='GR2M'",
# column(4, offset = 0,
# sliderInput("Event", label = "Select the target date:",
# min = as.POSIXct(.ShinyGR.args$SimPer[1L], tz = "UTC"),
# max = as.POSIXct(.ShinyGR.args$SimPer[2L], tz = "UTC"),
# value = as.POSIXct(.ShinyGR.args$SimPer[1L], tz = "UTC"),
# timeFormat = "%F",
# timezone = "+0000",
# animate = animationOptions(interval = 500),
# step = 3600 * 24 * 31)
# )
# ),
# conditionalPanel(condition = "input.PlotType == 'Model diagram' & input.TypeModel =='GR1A'",
# column(4, offset = 0,
# sliderInput("Event", label = "Select the target date:",
# min = as.POSIXct(.ShinyGR.args$SimPer[1L], tz = "UTC"),
# max = as.POSIXct(.ShinyGR.args$SimPer[2L], tz = "UTC"),
# value = as.POSIXct(.ShinyGR.args$SimPer[1L], tz = "UTC"),
# timeFormat = "%F",
# timezone = "+0000",
# animate = animationOptions(interval = 500),
# step = 3600 * 24 * 31 * 12)
# )
)
),
fluidRow(column(10, plotOutput("myPlot", width = "100%", height = "665px") ), #800
),
fluidRow(column(10, plotOutput("myPlot", width = "100%", height = "665px") ),
column(02, tableOutput("Criteria"))
)
)
)
),
# tabPanel(title = "How to use interface", fluidRow(column(6, includeMarkdown("wwwtab_ui.md"))),
# icon = icon("question-circle-o")),
tabPanel(title = "Functionalities" , fluidRow(column(6, includeMarkdown("www/tab_fun.md"))),
icon = icon("cog")),
tabPanel(title = "About" , fluidRow(column(6, includeMarkdown("www/tab_about.md")),
column(5, includeMarkdown("www/tab_authors.md"))),
icon = icon("navicon"))
icon = icon("navicon")),
tags$head(tags$style("#Criteria td{ position:relative; };"))
)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment