Commit bbad13a9 authored by Delaigue Olivier's avatar Delaigue Olivier
Browse files

v0.2.10.13 NEW: first test of GR2M model diagram in the GUI #14

Needs to be greatly improved!
Showing with 132 additions and 104 deletions
+132 -104
Package: airGRteaching Package: airGRteaching
Type: Package Type: Package
Title: Teaching Hydrological Modelling with the GR Rainfall-Runoff Models ('Shiny' Interface Included) Title: Teaching Hydrological Modelling with the GR Rainfall-Runoff Models ('Shiny' Interface Included)
Version: 0.2.10.12 Version: 0.2.10.13
Date: 2020-04-14 Date: 2020-04-14
Authors@R: c( Authors@R: c(
person("Olivier", "Delaigue", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-7668-8468"), email = "airGR@inrae.fr"), person("Olivier", "Delaigue", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-7668-8468"), email = "airGR@inrae.fr"),
......
...@@ -4,7 +4,7 @@ ...@@ -4,7 +4,7 @@
### 0.2.10.12 Release Notes (2020-04-14) ### 0.2.10.13 Release Notes (2020-04-14)
#### New features #### New features
......
...@@ -211,30 +211,30 @@ if (getRversion() >= "2.15.1") { ...@@ -211,30 +211,30 @@ if (getRversion() >= "2.15.1") {
} }
if (HydroModel != "GR2M") {
# Calcul des ordonnees SH2 de l' "hydrogramme unitaire cumule" UH2 # Calcul des ordonnees SH2 de l' "hydrogramme unitaire cumule" UH2
SH2 <- array(NA, 2*NH) SH2 <- array(NA, 2*NH)
for (i in 1:(2*NH)) { for (i in 1:(2*NH)) {
if (i <= 0) SH2[i] <- 0 if (i <= 0) SH2[i] <- 0
if (i > 0 & i < Param[4]) SH2[i] <- 0.5*(i/Param[4])^(D) 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 >= Param[4] & i < 2*Param[4]) SH2[i] <- 1 - (0.5*(2-i/Param[4])^(D))
if (i >= 2*Param[4]) SH2[i] <- 1 if (i >= 2*Param[4]) SH2[i] <- 1
} }
# Calcul des ordonnees UH2 de l' "hydrogramme unitaire discret" UH2 # Calcul des ordonnees UH2 de l' "hydrogramme unitaire discret" UH2
UH2 <- array(NA, 2*NH) UH2 <- array(NA, 2*NH)
for (j in 1:(2*NH)) { for (j in 1:(2*NH)) {
if (j == 1) { if (j == 1) {
UH2[j] <- SH2[j] UH2[j] <- SH2[j]
} else { } else {
UH2[j] <- SH2[j] - SH2[j-1] UH2[j] <- SH2[j] - SH2[j-1]
} }
}
# Parametres
max_UH2 <- log(sqrt(max(max(UH2)*OutputsModel$PR*0.1))+1)
} }
# Parametres
max_UH2 <- log(sqrt(max(max(UH2)*OutputsModel$PR*0.1))+1)
# -------------------------------------------------------------------------------- # --------------------------------------------------------------------------------
# PARTITIONNEMENT FENETRE GRAPHIQUE # PARTITIONNEMENT FENETRE GRAPHIQUE
# -------------------------------------------------------------------------------- # --------------------------------------------------------------------------------
...@@ -308,10 +308,12 @@ if (getRversion() >= "2.15.1") { ...@@ -308,10 +308,12 @@ if (getRversion() >= "2.15.1") {
# NEUTRALISATION DE P # NEUTRALISATION DE P
# -------------------------------------------------------------------------------- # --------------------------------------------------------------------------------
# Interception if (HydroModel != "GR2M") {
segments(x0 = xy_E[1]-50, x1 = xy_P[1]+50, # Interception
y0 = y_interception+tmp_decal, y1 = y_interception+tmp_decal) segments(x0 = xy_E[1]-50, x1 = xy_P[1]+50,
text(x = xy_P[1]+50, y = y_interception+20, labels = "Interception", pos = 4, font = 1, cex = 1.4) y0 = y_interception+tmp_decal, y1 = y_interception+tmp_decal)
text(x = xy_P[1]+50, y = y_interception+20, labels = "Interception", pos = 4, font = 1, cex = 1.4)
}
# P vers Pn # P vers Pn
segments(x0 = xy_P[1], x1 = xy_P[1], y0 = xy_P[2], y1 = y_interception+tmp_decal) segments(x0 = xy_P[1], x1 = xy_P[1], y0 = xy_P[2], y1 = y_interception+tmp_decal)
...@@ -340,11 +342,13 @@ if (getRversion() >= "2.15.1") { ...@@ -340,11 +342,13 @@ if (getRversion() >= "2.15.1") {
segments(x0 = xy_E[1], x1 = xy_E[1], segments(x0 = xy_E[1], x1 = xy_E[1],
y0 = y_interception, y1 = y_rendement) y0 = y_interception, y1 = y_rendement)
if (HydroModel != "GR2M") {
# Ecriture # Ecriture
plotrix::boxed.labels(x = xy_P[1], y = y_interception, labels = "Pn", plotrix::boxed.labels(x = xy_P[1], y = y_interception, labels = "Pn",
bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad) bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
plotrix::boxed.labels(x = xy_E[1], y = y_interception, labels = "En", plotrix::boxed.labels(x = xy_E[1], y = y_interception, labels = "En",
bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad) bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
}
# ETP # ETP
if (OutputsModel$PotEvap[i_pdt] != 0) { if (OutputsModel$PotEvap[i_pdt] != 0) {
...@@ -361,11 +365,11 @@ if (getRversion() >= "2.15.1") { ...@@ -361,11 +365,11 @@ if (getRversion() >= "2.15.1") {
} }
# Pn et Ps # Pn et Ps
if (HydroModel != "GR2M") {
points(x = x_Ps, y = y_rendement+1.2*tmp_decal, points(x = x_Ps, y = y_rendement+1.2*tmp_decal,
type = "p", pch = tri_B, col = col_P, type = "p", pch = tri_B, col = col_P,
cex = cex_tri(OutputsModel$Ps[i_pdt], fact = fact_triangle, max = cex_max_poly)) cex = cex_tri(OutputsModel$Ps[i_pdt], fact = fact_triangle, max = cex_max_poly))
}
points(x = x_PnPs, y = y_rendement+1.2*tmp_decal, points(x = x_PnPs, y = y_rendement+1.2*tmp_decal,
type = "p", pch = tri_B, col = col_P, type = "p", pch = tri_B, col = col_P,
cex = cex_tri(OutputsModel$Pn[i_pdt] - OutputsModel$Ps[i_pdt], fact = fact_triangle, max = cex_max_poly)) cex = cex_tri(OutputsModel$Pn[i_pdt] - OutputsModel$Ps[i_pdt], fact = fact_triangle, max = cex_max_poly))
...@@ -474,13 +478,15 @@ if (getRversion() >= "2.15.1") { ...@@ -474,13 +478,15 @@ if (getRversion() >= "2.15.1") {
cex = cex_tri(OutputsModel$PR[i_pdt]*0.9, fact = fact_triangle, max = cex_max_poly)) cex = cex_tri(OutputsModel$PR[i_pdt]*0.9, fact = fact_triangle, max = cex_max_poly))
} }
# Remplissage de UH1 if (HydroModel != "GR2M") {
PR_mat_UH1_lg <- ceiling(Param[4]) # Remplissage de UH1
PR_mat_UH1_id <- max(i_pdt-PR_mat_UH1_lg+1, 1):i_pdt PR_mat_UH1_lg <- ceiling(Param[4])
PR_mat_UH1 <- matrix(rep(c(rep(0, times = PR_mat_UH1_lg-length(PR_mat_UH1_id)+1), PR_mat_UH1_id <- max(i_pdt-PR_mat_UH1_lg+1, 1):i_pdt
OutputsModel$PR[PR_mat_UH1_id]), times = PR_mat_UH1_lg), PR_mat_UH1 <- matrix(rep(c(rep(0, times = PR_mat_UH1_lg-length(PR_mat_UH1_id)+1),
ncol = PR_mat_UH1_lg+1)[, -1L] OutputsModel$PR[PR_mat_UH1_id]), times = PR_mat_UH1_lg),
PR_mat_UH1[lower.tri(PR_mat_UH1)] <- 0 ncol = PR_mat_UH1_lg+1)[, -1L]
PR_mat_UH1[lower.tri(PR_mat_UH1)] <- 0
}
# -------------------------------------------------------------------------------- # --------------------------------------------------------------------------------
...@@ -494,14 +500,15 @@ if (getRversion() >= "2.15.1") { ...@@ -494,14 +500,15 @@ if (getRversion() >= "2.15.1") {
cex = cex_tri(OutputsModel$PR[i_pdt]*0.1, fact = fact_triangle, max = cex_max_poly)) cex = cex_tri(OutputsModel$PR[i_pdt]*0.1, fact = fact_triangle, max = cex_max_poly))
} }
if (HydroModel != "GR2M") {
# Remplissage de UH2 # Remplissage de UH2
PR_mat_UH2_lg <- ceiling(Param[4]*2) PR_mat_UH2_lg <- ceiling(Param[4]*2)
PR_mat_UH2_id <- max(i_pdt-PR_mat_UH2_lg+1, 1):i_pdt PR_mat_UH2_id <- max(i_pdt-PR_mat_UH2_lg+1, 1):i_pdt
PR_mat_UH2 <- matrix(rep(c(rep(0, times = PR_mat_UH2_lg-length(PR_mat_UH2_id)+1), PR_mat_UH2 <- matrix(rep(c(rep(0, times = PR_mat_UH2_lg-length(PR_mat_UH2_id)+1),
OutputsModel$PR[PR_mat_UH2_id]), times = PR_mat_UH2_lg), OutputsModel$PR[PR_mat_UH2_id]), times = PR_mat_UH2_lg),
ncol = PR_mat_UH2_lg+1)[, -1L] ncol = PR_mat_UH2_lg+1)[, -1L]
PR_mat_UH2[lower.tri(PR_mat_UH2)] <- 0 PR_mat_UH2[lower.tri(PR_mat_UH2)] <- 0
}
} }
...@@ -578,49 +585,51 @@ if (getRversion() >= "2.15.1") { ...@@ -578,49 +585,51 @@ if (getRversion() >= "2.15.1") {
y0 = y_routage, y1 = y_routage) y0 = y_routage, y1 = y_routage)
} }
# Q9 if (HydroModel != "GR2M") {
if (OutputsModel$Q9[i_pdt] != 0) { # Q9
points(x = xy_Q9[1], y = xy_Q9[2]+tmp_decal, if (OutputsModel$Q9[i_pdt] != 0) {
type = "p", pch = tri_B, col = col_P, points(x = xy_Q9[1], y = xy_Q9[2]+tmp_decal,
cex = cex_tri(OutputsModel$Q9[i_pdt], fact = fact_triangle, max = cex_max_poly))
if (HydroModel == "GR6J") {
# Q9 exp
points(x = xy_Q9[1]*0.80, y = xy_Q9[1]*0.73,
type = "p", pch = tri_B, col = col_P, type = "p", pch = tri_B, col = col_P,
cex = cex_tri(OutputsModel$Q9[i_pdt]*0.4, fact = fact_triangle, max = cex_max_poly)) cex = cex_tri(OutputsModel$Q9[i_pdt], fact = fact_triangle, max = cex_max_poly))
# Q9 rout if (HydroModel == "GR6J") {
points(x = xy_Q9[1]*1.30, y = xy_Q9[1]*0.73, # Q9 exp
points(x = xy_Q9[1]*0.80, y = xy_Q9[1]*0.73,
type = "p", pch = tri_B, col = col_P,
cex = cex_tri(OutputsModel$Q9[i_pdt]*0.4, fact = fact_triangle, max = cex_max_poly))
# Q9 rout
points(x = xy_Q9[1]*1.30, y = xy_Q9[1]*0.73,
type = "p", pch = tri_B, col = col_P,
cex = cex_tri(OutputsModel$Q9[i_pdt]*0.6, fact = fact_triangle, max = cex_max_poly))
# QrExp
plotrix::boxed.labels(x = xy_Q9[1]*0.55, y = y_routage, labels = "QrExp", bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
}
}
plotrix::boxed.labels(x = xy_Q9[1], y = xy_Q9[2], labels = "Q9",
bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
# Q1
if (OutputsModel$Q1[i_pdt] != 0) {
points(x = xy_Q1[1], y = xy_Q1[2]+tmp_decal,
type = "p", pch = tri_B, col = col_P, type = "p", pch = tri_B, col = col_P,
cex = cex_tri(OutputsModel$Q9[i_pdt]*0.6, fact = fact_triangle, max = cex_max_poly)) cex = cex_tri(OutputsModel$Q1[i_pdt], fact = fact_triangle, max = cex_max_poly))
# QrExp segments(x0 = xy_Q[1], x1 = xy_Q1[1], y0 = y_routage, y1 = y_routage)
plotrix::boxed.labels(x = xy_Q9[1]*0.55, y = y_routage, labels = "QrExp", bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
} }
plotrix::boxed.labels(x = xy_Q1[1], y = xy_Q1[2], labels = "Q1", bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
# Valeur de Qd
if (OutputsModel$QD[i_pdt] != 0) {
points(x = xy_Q[1]+30, y = y_routage,
type = "p", pch = tri_L, col = col_P,
cex = cex_tri(OutputsModel$QD[i_pdt], fact = fact_triangle, max = cex_max_poly))
}
# Qd
plotrix::boxed.labels(x = xy_Q1[1], y = y_routage, labels = "Qd", bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
} }
plotrix::boxed.labels(x = xy_Q9[1], y = xy_Q9[2], labels = "Q9",
bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
# Q1
if (OutputsModel$Q1[i_pdt] != 0) {
points(x = xy_Q1[1], y = xy_Q1[2]+tmp_decal,
type = "p", pch = tri_B, col = col_P,
cex = cex_tri(OutputsModel$Q1[i_pdt], fact = fact_triangle, max = cex_max_poly))
segments(x0 = xy_Q[1], x1 = xy_Q1[1], y0 = y_routage, y1 = y_routage)
}
plotrix::boxed.labels(x = xy_Q1[1], y = xy_Q1[2], labels = "Q1", bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
# Valeur de Qd
if (OutputsModel$QD[i_pdt] != 0) {
points(x = xy_Q[1]+30, y = y_routage,
type = "p", pch = tri_L, col = col_P,
cex = cex_tri(OutputsModel$QD[i_pdt], fact = fact_triangle, max = cex_max_poly))
}
# Qd
plotrix::boxed.labels(x = xy_Q1[1], y = y_routage, labels = "Qd", bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
# -------------------------------------------------------------------------------- # --------------------------------------------------------------------------------
# RESERVOIR DE ROUTAGE # RESERVOIR DE ROUTAGE
# -------------------------------------------------------------------------------- # --------------------------------------------------------------------------------
...@@ -629,6 +638,9 @@ if (getRversion() >= "2.15.1") { ...@@ -629,6 +638,9 @@ if (getRversion() >= "2.15.1") {
tmp_triche <- 0#80 tmp_triche <- 0#80
# Reservoir de routage # Reservoir de routage
if (HydroModel == "GR2M") {
Param[3] <- 600
}
rect(xleft = xy_min_ROUT[1], xright = xy_min_ROUT[1]+base_res, 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, ybottom = xy_min_ROUT[2], ytop = xy_min_ROUT[2]+OutputsModel$Rout[i_pdt]*fact_res+tmp_triche,
col = col_SR, border = NA) col = col_SR, border = NA)
...@@ -646,26 +658,29 @@ if (getRversion() >= "2.15.1") { ...@@ -646,26 +658,29 @@ if (getRversion() >= "2.15.1") {
segments(x0 = xy_min_ROUT[1]+base_res/2, x1 = xy_Q[1], segments(x0 = xy_min_ROUT[1]+base_res/2, x1 = xy_Q[1],
y0 = y_routage, y1 = y_routage) y0 = y_routage, y1 = y_routage)
# Qr
if (HydroModel != "GR6J") {
plotrix::boxed.labels(x = xy_min_ROUT[1]+base_res/2, y = y_routage, labels = "Qr",
bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
}
if (HydroModel == "GR6J") {
plotrix::boxed.labels(x = xy_min_ROUT[1]+base_res/1.5, y = (xy_min_ROUT[2]+y_routage)/2, labels = "Qr",
bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
}
# Valeur de Qr if (HydroModel != "GR2M") {
if (OutputsModel$QR[i_pdt] != 0) { # Qr
if (HydroModel != "GR6J") { if (HydroModel != "GR6J") {
points(x = xy_Q[1]-100, y = y_routage, plotrix::boxed.labels(x = xy_min_ROUT[1]+base_res/2, y = y_routage, labels = "Qr",
type = "p", pch = tri_R, col = col_P, bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
cex = cex_tri(OutputsModel$QR[i_pdt], fact = fact_triangle, max = cex_max_poly)) }
} else { if (HydroModel == "GR6J") {
points(x = xy_min_ROUT[1]+base_res/2, y = (xy_min_ROUT[2]+y_routage)/2, plotrix::boxed.labels(x = xy_min_ROUT[1]+base_res/1.5, y = (xy_min_ROUT[2]+y_routage)/2, labels = "Qr",
type = "p", pch = tri_B, col = col_P, bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
cex = cex_tri(OutputsModel$QR[i_pdt], fact = fact_triangle, max = cex_max_poly)) }
# Valeur de Qr
if (OutputsModel$QR[i_pdt] != 0) {
if (HydroModel != "GR6J") {
points(x = xy_Q[1]-100, y = y_routage,
type = "p", pch = tri_R, col = col_P,
cex = cex_tri(OutputsModel$QR[i_pdt], fact = fact_triangle, max = cex_max_poly))
} else {
points(x = xy_min_ROUT[1]+base_res/2, y = (xy_min_ROUT[2]+y_routage)/2,
type = "p", pch = tri_B, col = col_P,
cex = cex_tri(OutputsModel$QR[i_pdt], fact = fact_triangle, max = cex_max_poly))
}
} }
} }
...@@ -737,6 +752,7 @@ if (getRversion() >= "2.15.1") { ...@@ -737,6 +752,7 @@ if (getRversion() >= "2.15.1") {
# EXCHANGE # EXCHANGE
# -------------------------------------------------------------------------------- # --------------------------------------------------------------------------------
if (HydroModel != "GR2M") {
# Actual exchange Q9 # Actual exchange Q9
arrows(x0 = xy_min_ROUT[1]+base_res, x1 = 1025, arrows(x0 = xy_min_ROUT[1]+base_res, x1 = 1025,
y0 = y_Ech_Q9 , y1 = y_Ech_Q9, y0 = y_Ech_Q9 , y1 = y_Ech_Q9,
...@@ -754,6 +770,18 @@ if (getRversion() >= "2.15.1") { ...@@ -754,6 +770,18 @@ if (getRversion() >= "2.15.1") {
points(x = xy_Q1[1]+100, y = y_Ech_Q1, points(x = xy_Q1[1]+100, y = y_Ech_Q1,
type = "p", pch = pch, col = col_P, type = "p", pch = pch, col = col_P,
cex = cex_tri(OutputsModel$AExch2[i_pdt], fact = fact_triangle, max = cex_max_poly)) cex = cex_tri(OutputsModel$AExch2[i_pdt], fact = fact_triangle, max = cex_max_poly))
}
if (HydroModel == "GR2M") {
# Actual exchange
arrows(x0 = xy_min_ROUT[1]+base_res, x1 = 1025,
y0 = y_Ech_Q9 , y1 = y_Ech_Q9,
length = 0.075, angle = 20, code = 3)
pch <- ifelse(OutputsModel$Exch[i_pdt] < 0, tri_R, tri_L)
points(x = xy_min_ROUT[1]+base_res+130, y = y_Ech_Q9,
type = "p", pch = pch, col = col_P,
cex = cex_tri(OutputsModel$Exch[i_pdt], fact = fact_triangle, max = cex_max_poly))
}
# Actual exchange Q9 exp. # Actual exchange Q9 exp.
if (HydroModel%in% c("GR6J")) { if (HydroModel%in% c("GR6J")) {
......
Supports Markdown
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