Commit 094c10f2 authored by unknown's avatar unknown
Browse files

v0.1.1.10 shiny.SimGR(à updated (the package now needs airGR 1.0.5.22)

parent e85382a9
Package: airGRteaching
Type: Package
Title: Tools to Simplify the Use of the airGR Hydrological Package for Education (including a Shiny Application)
Version: 0.1.1.9
Title: Tools to Simplify the Use of the airGR Hydrological Package for Education (Including a Shiny Application)
Version: 0.1.1.10
Date: 2017-02-16
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")))
Depends: airGR (>= 1.0.5.4)
Depends: airGR (>= 1.0.5.22)
Imports: xts, dygraphs, shiny, plotrix
Description: Providing 3 functions which allows to run very simply the airGR hydrological modelling steps. Also providing plotting function to help student to explore observed data and to interprete the results of calibration and simulation of the GR hydrological models.
Description: Providing 3 functions which allows to run very simply the airGR hydrological modelling steps. Also providing plotting function to help student to explore observed data and to interpret the results of calibration and simulation of the GR hydrological models.
License: GPL-2
NeedsCompilation: no
URL: https://webgr.irstea.fr/airGR/?lang=en
......
......@@ -6,6 +6,15 @@
##################################################################################
##################################################################################
if (getRversion() >= "2.15.1") {
utils::globalVariables(c(".SimGR.args"))
utils::suppressForeignCheck(c(".SimGR.args"))
utils::globalVariables(c(".SimGR4J.args"))
utils::suppressForeignCheck(c(".SimGR4J.args"))
}
.TypeModelGR <- function(x) {
if (!is.list(x)) {
x <- list(TypeModel = x)
......@@ -35,146 +44,146 @@
DiagramGR4J <- function(OutputsModel, Param, SimPer, EventDate) {
##################################################################################
# PARAMETRES
##################################################################################
# Parametres
mgp <- c(3,0.5,0)
col_P <- "black"
col_E <- "forestgreen"
col_Q <- "black"
col_modele <- "gray95"
nom_modele <- "GR4J"
xy_E <- c(250, 970)
xy_P <- c(600, 970)
xy_Q <- c(750, 30)
x_Ps <- 440
x_PnPs <- 700
y_interception <- 900
y_rendement <- 815
y_percolation <- 575
xy_Q9 <- c(400,310)
xy_Q1 <- c(800,310)
y_routage <- 100
fact_res <- 200/800
base_res <- 300
NH <- 10
xy_min_PROD <- c(200, 610)
xy_min_ROUT <- c(250, 150)
y_entreeHU <- 500
xy_HU1 <- c(500, 420)
xy_HU2 <- c(900, 420)
y_Ech_Q1 <- 200
y_Ech_Q9 <- 180
y_out_aubes <- 400
D <- 5/2
xpad <- 1.5
ypad <- 1.5
max_triangle <- max(unlist(OutputsModel[c("Perc", "PR", "Q9", "Q1", "QR", "QD")]))
fact_var <- 40
fact_triangle <- 25
radius1 <- 0
radius2 <- 60
cex_max_poly <- 0.005
# Pas de temps
dates_deb <- EventDate#as.Date(OutputsModel$DatesR)[1]
n_pdt <- length(which(OutputsModel$DatesR >= EventDate & OutputsModel$DatesR <= SimPer[2L]))#365
dates_epi <- OutputsModel$DatesR[which(OutputsModel$DatesR >= EventDate & OutputsModel$DatesR <= SimPer[2L])]#as.Date(OutputsModel$DatesR[1:n_pdt])
dates_epi <- as.Date(dates_epi)
i_pdt <- which(format(OutputsModel$DatesR, "%Y%m%d") == format(EventDate, "%Y%m%d")) #150#
version <- "V5"
create_polygon <- function(x_center = x_center, y_center = y_center, cex = cex, dir = "D") {
xy <- matrix(NA, nrow = 3, ncol = 2)
cex <- log(cex * fact_triangle + 1)
if(dir == "D") {
xy[,1] <- c(x_center-cex/2, x_center-cex/2, x_center+cex/2)
xy[,2] <- c(y_center-cex/2, y_center+cex/2, y_center)
##################################################################################
# PARAMETRES
##################################################################################
# Parametres
mgp <- c(3,0.5,0)
col_P <- "black"
col_E <- "forestgreen"
col_Q <- "black"
col_modele <- "gray95"
nom_modele <- "GR4J"
xy_E <- c(250, 970)
xy_P <- c(600, 970)
xy_Q <- c(750, 30)
x_Ps <- 440
x_PnPs <- 700
y_interception <- 900
y_rendement <- 815
y_percolation <- 575
xy_Q9 <- c(400,310)
xy_Q1 <- c(800,310)
y_routage <- 100
fact_res <- 200/800
base_res <- 300
NH <- 10
xy_min_PROD <- c(200, 610)
xy_min_ROUT <- c(250, 150)
y_entreeHU <- 500
xy_HU1 <- c(500, 420)
xy_HU2 <- c(900, 420)
y_Ech_Q1 <- 200
y_Ech_Q9 <- 180
y_out_aubes <- 400
D <- 5/2
xpad <- 1.5
ypad <- 1.5
max_triangle <- max(unlist(OutputsModel[c("Perc", "PR", "Q9", "Q1", "QR", "QD")]))
fact_var <- 40
fact_triangle <- 25
radius1 <- 0
radius2 <- 60
cex_max_poly <- 0.005
# Pas de temps
dates_deb <- EventDate#as.Date(OutputsModel$DatesR)[1]
n_pdt <- length(which(OutputsModel$DatesR >= EventDate & OutputsModel$DatesR <= SimPer[2L]))#365
dates_epi <- OutputsModel$DatesR[which(OutputsModel$DatesR >= EventDate & OutputsModel$DatesR <= SimPer[2L])]#as.Date(OutputsModel$DatesR[1:n_pdt])
dates_epi <- as.Date(dates_epi)
i_pdt <- which(format(OutputsModel$DatesR, "%Y%m%d") == format(EventDate, "%Y%m%d")) #150#
version <- "V5"
create_polygon <- function(x_center = x_center, y_center = y_center, cex = cex, dir = "D") {
xy <- matrix(NA, nrow = 3, ncol = 2)
cex <- log(cex * fact_triangle + 1)
if(dir == "D") {
xy[,1] <- c(x_center-cex/2, x_center-cex/2, x_center+cex/2)
xy[,2] <- c(y_center-cex/2, y_center+cex/2, y_center)
}
if(dir == "G") {
xy[,1] <- c(x_center+cex/2, x_center+cex/2, x_center-cex/2)
xy[,2] <- c(y_center-cex/2, y_center+cex/2, y_center)
}
if(dir == "H") {
xy[,1] <- c(x_center-cex/2, x_center+cex/2, x_center)
xy[,2] <- c(y_center-cex/2, y_center-cex/2, y_center+cex/2)
}
if(dir == "B") {
xy[,1] <- c(x_center-cex/2, x_center+cex/2, x_center)
xy[,2] <- c(y_center+cex/2, y_center+cex/2, y_center-cex/2)
}
return(xy)
}
if(dir == "G") {
xy[,1] <- c(x_center+cex/2, x_center+cex/2, x_center-cex/2)
xy[,2] <- c(y_center-cex/2, y_center+cex/2, y_center)
# plot(x=0, type="n", xlim=c(0,max_triangle), ylim=c(0,max_triangle), xlab="", ylab="")
# for(i in seq(1,as.integer(max_triangle), 2)) {
# tmp_triangle <- create_polygon(x_center=i, y_center=i,
# cex=log((i/2)+1), dir="H")
# polygon(x=tmp_triangle[,1], y=tmp_triangle[,2], border=NA, col="royalblue")
# }
##################################################################################
# HUs
##################################################################################
# Calcul des ordonnees SH1 de l' "hydrogramme unitaire cumule" HU1
SH1 <- array(NA, NH)
for (i in 1:NH) {
if (i <= 0) SH1[i] <- 0
if (i > 0 & i < Param[4]) SH1[i] <- (i/Param[4])^(D)
if (i >= Param[4]) SH1[i] <- 1
}
if(dir == "H") {
xy[,1] <- c(x_center-cex/2, x_center+cex/2, x_center)
xy[,2] <- c(y_center-cex/2, y_center-cex/2, y_center+cex/2)
# Calcul des ordonnees UH1 de l' "hydrogramme unitaire discret" HU1
UH1 <- array(NA, NH)
for (j in 1:NH) {
if (j == 1) {
UH1[j] <- SH1[j]
} else {
UH1[j] <- SH1[j] - SH1[j-1]
}
}
if(dir == "B") {
xy[,1] <- c(x_center-cex/2, x_center+cex/2, x_center)
xy[,2] <- c(y_center+cex/2, y_center+cex/2, y_center-cex/2)
# Calcul des ordonnees SH2 de l' "hydrogramme unitaire cumule" HU2
SH2 <- array(NA, 2*NH)
for (i in 1:(2*NH)) {
if(i <= 0) SH2[i] <- 0
if(i > 0 & i < Param[4]) SH2[i] <- 0.5*(i/Param[4])^(D)
if(i > Param[4] & i < 2*Param[4]) SH2[i] <- 1 - (0.5*(2-i/Param[4])^(D))
if(i >= 2*Param[4]) SH2[i] <- 1
}
return(xy)
}
# plot(x=0, type="n", xlim=c(0,max_triangle), ylim=c(0,max_triangle), xlab="", ylab="")
# for(i in seq(1,as.integer(max_triangle), 2)) {
# tmp_triangle <- create_polygon(x_center=i, y_center=i,
# cex=log((i/2)+1), dir="H")
# polygon(x=tmp_triangle[,1], y=tmp_triangle[,2], border=NA, col="royalblue")
# }
##################################################################################
# HUs
##################################################################################
# Calcul des ordonnees SH1 de l' "hydrogramme unitaire cumule" HU1
SH1 <- array(NA, NH)
for (i in 1:NH) {
if (i <= 0) SH1[i] <- 0
if (i > 0 & i < Param[4]) SH1[i] <- (i/Param[4])^(D)
if (i >= Param[4]) SH1[i] <- 1
}
# Calcul des ordonnees UH1 de l' "hydrogramme unitaire discret" HU1
UH1 <- array(NA, NH)
for (j in 1:NH) {
if (j == 1) {
UH1[j] <- SH1[j]
} else {
UH1[j] <- SH1[j] - SH1[j-1]
}
}
# Calcul des ordonnees SH2 de l' "hydrogramme unitaire cumule" HU2
SH2 <- array(NA, 2*NH)
for (i in 1:(2*NH)) {
if(i <= 0) SH2[i] <- 0
if(i > 0 & i < Param[4]) SH2[i] <- 0.5*(i/Param[4])^(D)
if(i > Param[4] & i < 2*Param[4]) SH2[i] <- 1 - (0.5*(2-i/Param[4])^(D))
if(i >= 2*Param[4]) SH2[i] <- 1
}
# Calcul des ordonnees UH2 de l' "hydrogramme unitaire discret" HU2
UH2 <- array(NA, 2*NH)
for (j in 1:(2*NH)) {
if (j == 1) {
UH2[j] <- SH2[j]
} else {
UH2[j] <- SH2[j] - SH2[j-1]
}
}
# Parametres
max_UH1 <- log(sqrt(max(max(UH1)*OutputsModel$PR*0.9))+1)
max_UH2 <- log(sqrt(max(max(UH2)*OutputsModel$PR*0.1))+1)
##################################################################################
# BOUCLE SUR PAS DE TEMPS
##################################################################################
# for(i_pdt in 1:n_pdt) {
# Calcul des ordonnees UH2 de l' "hydrogramme unitaire discret" HU2
UH2 <- array(NA, 2*NH)
for (j in 1:(2*NH)) {
if (j == 1) {
UH2[j] <- SH2[j]
} else {
UH2[j] <- SH2[j] - SH2[j-1]
}
}
# Parametres
max_UH1 <- log(sqrt(max(max(UH1)*OutputsModel$PR*0.9))+1)
max_UH2 <- log(sqrt(max(max(UH2)*OutputsModel$PR*0.1))+1)
##################################################################################
# BOUCLE SUR PAS DE TEMPS
##################################################################################
# for(i_pdt in 1:n_pdt) {
# layout(matrix(1:3, nrow=3, ncol=1, byrow=TRUE), height=c(5,12,5))
layout(matrix(c(1:4, 4, 4), nrow = 3, ncol = 2, byrow = FALSE), widths = c(1.0, 0.6))
# layout(matrix(1:3, nrow=3, ncol=1, byrow=TRUE), height=c(5,12,5))
layout(matrix(c(1:4, 4, 4), nrow = 3, ncol = 2, byrow = FALSE), widths = c(1.0, 0.6))
##################################################################################
......@@ -433,35 +442,35 @@ layout(matrix(c(1:4, 4, 4), nrow = 3, ncol = 2, byrow = FALSE), widths = c(1.0,
segments(x0=xy_HU1[1], y0=xy_HU1[2], x1=(cos(angle1)*radius2) + xy_HU1[1], y1=(sin(angle1)*radius2) + xy_HU1[2])
segments(x0=xy_HU1[1], y0=xy_HU1[2], x1=(cos(angle2)*radius2) + xy_HU1[1], y1=(sin(angle2)*radius2) + xy_HU1[2])
}
# Fleche sortant de HU1
arrows(x0=xy_Q9[1]+30, y0=y_out_aubes+10, x1=xy_Q9[1], y1=y_out_aubes, length=0.075, angle=20)
# Sorties de HU1
if(OutputsModel$Q9[i_pdt] != 0) {
tmp_triangle <- create_polygon(x_center=xy_Q9[1], y_center=xy_Q9[2]+tmp_decal,
cex=OutputsModel$Q9[i_pdt]*max_triangle / cex_max_poly, dir="B")
polygon(x=tmp_triangle[,1], y=tmp_triangle[,2], border=NA, col="royalblue")
}
# Q9
boxed.labels(x=xy_Q9[1], y=xy_Q9[2], labels="Q9", col="black", bg=col_modele, border=NA, xpad=xpad, ypad=ypad)
##################################################################################
# HYDROGRAMME UNITAIRE 2
##################################################################################
# Entree de HU2
if(OutputsModel$PR[i_pdt] != 0) {
tmp_triangle <- create_polygon(x_center=xy_Q1[1], y_center=y_entreeHU+tmp_decal/2,
cex=OutputsModel$PR[i_pdt]*0.1*max_triangle / cex_max_poly, dir="B")
polygon(x=tmp_triangle[,1], y=tmp_triangle[,2], border=NA, col="royalblue")
}
# Fleche vers HU2
arrows(x0=xy_Q1[1], y0=y_entreeHU, x1=xy_Q1[1]+30, y1=y_entreeHU-10, length=0.075, angle=20)
# Remplissage de HU2
tmp_UH2 <- OutputsModel$PR[i_pdt]*0.1*fact_var*UH2
tmp_UH2[1] <- 0
......@@ -478,42 +487,42 @@ layout(matrix(c(1:4, 4, 4), nrow = 3, ncol = 2, byrow = FALSE), widths = c(1.0,
angle1 <- angle2
angle2 <- temp
}
# if (radius1 > radius2) {
# temp <- radius1
# radius1 <- radius2
# radius2 <- temp
# }
# if (radius1 > radius2) {
# temp <- radius1
# radius1 <- radius2
# radius2 <- temp
# }
angles <- seq(angle1, angle2, by=0.01)
angles[length(angles)] <- angle2
xpos <- c(cos(angles) * radius1, cos(rev(angles)) * radius2) + xy_HU2[1]
ypos <- c(sin(angles) * radius1, sin(rev(angles)) * radius2) + xy_HU2[2]
polygon(xpos, ypos, col="royalblue", border=NA)
}
#
# # Roue vide
# for (i_dir in 1:(2*NH)) {
# angle1 <- sector_edges[i_dir]
# angle2 <- sector_edges[i_dir+1]
# radius1 <- 0
# radius2 <- 60
# if (angle1 > angle2) {
# temp <- angle1
# angle1 <- angle2
# angle2 <- temp
# }
# if (radius1 > radius2) {
# temp <- radius1
# radius1 <- radius2
# radius2 <- temp
# }
#
# segments(x0=xy_HU2[1], y0=xy_HU2[2], x1=(cos(angle1)*radius2) + xy_HU2[1], y1=(sin(angle1)*radius2) + xy_HU2[2])
# segments(x0=xy_HU2[1], y0=xy_HU2[2], x1=(cos(angle2)*radius2) + xy_HU2[1], y1=(sin(angle2)*radius2) + xy_HU2[2])
# }
#
#
# # Roue vide
# for (i_dir in 1:(2*NH)) {
# angle1 <- sector_edges[i_dir]
# angle2 <- sector_edges[i_dir+1]
# radius1 <- 0
# radius2 <- 60
# if (angle1 > angle2) {
# temp <- angle1
# angle1 <- angle2
# angle2 <- temp
# }
# if (radius1 > radius2) {
# temp <- radius1
# radius1 <- radius2
# radius2 <- temp
# }
#
# segments(x0=xy_HU2[1], y0=xy_HU2[2], x1=(cos(angle1)*radius2) + xy_HU2[1], y1=(sin(angle1)*radius2) + xy_HU2[2])
# segments(x0=xy_HU2[1], y0=xy_HU2[2], x1=(cos(angle2)*radius2) + xy_HU2[1], y1=(sin(angle2)*radius2) + xy_HU2[2])
# }
#
# Fleche sortant de HU2
arrows(x0=xy_Q1[1]+30, y0=y_out_aubes+10, x1=xy_Q1[1], y1=y_out_aubes, length=0.075, angle=20)
# Sorties de HU2
if(OutputsModel$Q1[i_pdt] != 0) {
tmp_triangle <- create_polygon(x_center=xy_Q1[1], y_center=xy_Q1[2]+tmp_decal,
......@@ -521,32 +530,32 @@ layout(matrix(c(1:4, 4, 4), nrow = 3, ncol = 2, byrow = FALSE), widths = c(1.0,
polygon(x=tmp_triangle[,1], y=tmp_triangle[,2], border=NA, col="royalblue")
}
segments(x0=xy_Q[1], x1=xy_Q1[1], y0=y_routage, y1=y_routage)
# Q1
boxed.labels(x=xy_Q1[1], y=xy_Q1[2], labels="Q1", col="black", bg=col_modele, border=NA, xpad=xpad, ypad=ypad)
# Echange
arrows(x0=xy_Q1[1], y0=y_Ech_Q1, x1=1025,y1=y_Ech_Q1, length=0.075, angle=20)
arrows(x0=1025, y0=y_Ech_Q1, x1=xy_Q1[1],y1=y_Ech_Q1, length=0.075, angle=20)
# Valeur de Qd
if(OutputsModel$QD[i_pdt] != 0) {
tmp_triangle <- create_polygon(x_center=xy_Q[1]+10, y_center=y_routage,
cex=OutputsModel$QD[i_pdt]*max_triangle / cex_max_poly, dir="G")
polygon(x=tmp_triangle[,1], y=tmp_triangle[,2], border=NA, col="royalblue")
}
# Qd
boxed.labels(x=xy_Q1[1], y=y_routage, labels="Qd", col="black", bg=col_modele, border=NA, xpad=xpad, ypad=ypad)
##################################################################################
# RESERVOIR DE ROUTAGE
##################################################################################
# Triche pour la taille du reservoire de routage
tmp_triche <- 80
# Reservoir de routage
rect(xleft=xy_min_ROUT[1], xright=xy_min_ROUT[1]+base_res, ybottom=xy_min_ROUT[2],
ytop=xy_min_ROUT[2]+OutputsModel$Rout[i_pdt]*fact_res+tmp_triche, col="royalblue", border=NA)
......@@ -554,46 +563,46 @@ layout(matrix(c(1:4, 4, 4), nrow = 3, ncol = 2, byrow = FALSE), widths = c(1.0,
segments(x0=xy_min_ROUT[1], x1=xy_min_ROUT[1], y0=xy_min_ROUT[2], y1=xy_min_ROUT[2]+Param[3]*fact_res+tmp_triche, col="black")
segments(x0=xy_min_ROUT[1]+base_res, x1=xy_min_ROUT[1]+base_res, y0=xy_min_ROUT[2], y1=xy_min_ROUT[2]+Param[3]*fact_res+tmp_triche, col="black")
text(x=50, y=xy_min_ROUT[2]+Param[1]*fact_res/3, "routing store", cex=1.4, col="black", pos=4)
# Sorties du reservoir
segments(x0=xy_min_ROUT[1]+base_res/2, x1=xy_min_ROUT[1]+base_res/2, y0=xy_min_ROUT[2], y1=y_routage, col="black")
segments(x0=xy_min_ROUT[1]+base_res/2, x1=xy_Q[1], y0=y_routage, y1=y_routage)
# Echange
arrows(x0=xy_min_ROUT[1]+base_res, y0=y_Ech_Q9, x1=1025,y1=y_Ech_Q9, length=0.075, angle=20)
arrows(x0=1025, y0=y_Ech_Q9, x1=xy_min_ROUT[1]+base_res,y1=y_Ech_Q9, length=0.075, angle=20)
# Qr
boxed.labels(x=xy_min_ROUT[1]+base_res/2, y=y_routage, labels="Qr", col="black", bg=col_modele, border=NA, xpad=xpad, ypad=ypad)
# Valeur de Qr
if(OutputsModel$QR[i_pdt] != 0) {
tmp_triangle <- create_polygon(x_center=xy_Q[1]-50, y_center=y_routage,
cex=OutputsModel$QR[i_pdt]*max_triangle / cex_max_poly, dir="D")
polygon(x=tmp_triangle[,1], y=tmp_triangle[,2], border=NA, col="royalblue")
}
# Q final
segments(x0=xy_Q[1], x1=xy_Q[1], y0=y_routage, y1=xy_Q[2]+10)
# Echange
tmp_Exch <- Param[2]*(OutputsModel$Rout[i_pdt]/Param[3])^(7/2)
if(tmp_Exch < 0) {
tmp_sens <- "G"
} else {
tmp_sens <- "D"
}
segments(x0=xy_min_ROUT[1]+base_res/2, x1=-10, y0=y_routage+40, y1=y_routage+40, col="black")
tmp_triangle <- create_polygon(x_center=xy_min_ROUT[1]-base_res/2, y_center=y_routage+40, cex=tmp_Exch*50, dir=tmp_sens)
polygon(x=tmp_triangle[,1], y=tmp_triangle[,2], border=NA, col="royalblue")
# Echange
tmp_Exch <- Param[2]*(OutputsModel$Rout[i_pdt]/Param[3])^(7/2)
if(tmp_Exch < 0) {
tmp_sens <- "G"
} else {
tmp_sens <- "D"
}
segments(x0=xy_min_ROUT[1]+base_res/2, x1=-10, y0=y_routage+40, y1=y_routage+40, col="black")
tmp_triangle <- create_polygon(x_center=xy_min_ROUT[1]-base_res/2, y_center=y_routage+40, cex=tmp_Exch*50, dir=tmp_sens)
polygon(x=tmp_triangle[,1], y=tmp_triangle[,2], border=NA, col="royalblue")
}
......@@ -3,7 +3,7 @@ shiny.SimGR <- function(ObsBV = NULL, DatesR = NULL, Precip = NULL, PotEvap = NU
Param = c(200, 0, 100, 2), WupPer = NULL, SimPer = NULL) {
.GlobalEnv$.SimGR.args <- list(ObsBV = ObsBV,
.GlobalEnv$.SimGR.args <- list(ObsBV = as.list(ObsBV),
DatesR = DatesR, Precip = Precip, PotEvap = PotEvap, Qobs = Qobs, TempMean = TempMean,
ZInputs = ZInputs, HypsoData = HypsoData, NLayers = NLayers,
Param = Param, WupPer = WupPer, SimPer = SimPer)
......
......@@ -2,133 +2,217 @@
shinyServer(function(input, output, session) {
# vars <- reactiveValues(counter = 0)
#
#
# observe({
# if (!is.null(input$CalButton)) {
# input$CalButton
# isolate({
# vars$counter <- vars$counter + 1
# print(vars$counter)
# })
# }
# })
#
# label <- reactive({
# if (!is.null(input$CalButton)) {
# if (vars$counter >= 2) {
# label <- "new label"
# updateActionButton(session, "CalButton", label = label)
# print(label)
# } else {
# label <- "old label"
# updateActionButton(session, "CalButton", label = label)
# print(label)
# }
# }
# })
#
output$myPlot <- renderPlot({
getRES <- reactive({