From bbad13a93ee5a70a647cc52f77667fb07f42d517 Mon Sep 17 00:00:00 2001
From: Delaigue Olivier <olivier.delaigue@irstea.fr>
Date: Tue, 14 Apr 2020 17:02:04 +0200
Subject: [PATCH] v0.2.10.13 NEW: first test of GR2M model diagram in the GUI
#14 Needs to be greatly improved!
---
DESCRIPTION | 2 +-
NEWS.md | 2 +-
R/Utils.R | 232 +++++++++++++++++++++++++++++-----------------------
3 files changed, 132 insertions(+), 104 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 3f365d0..f166493 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
Package: airGRteaching
Type: Package
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
Authors@R: c(
person("Olivier", "Delaigue", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-7668-8468"), email = "airGR@inrae.fr"),
diff --git a/NEWS.md b/NEWS.md
index a3bed4b..2a5a960 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,7 +4,7 @@
-### 0.2.10.12 Release Notes (2020-04-14)
+### 0.2.10.13 Release Notes (2020-04-14)
#### New features
diff --git a/R/Utils.R b/R/Utils.R
index 0f1df34..19b41f4 100644
--- a/R/Utils.R
+++ b/R/Utils.R
@@ -211,30 +211,30 @@ if (getRversion() >= "2.15.1") {
}
-
- # Calcul des ordonnees SH2 de l' "hydrogramme unitaire cumule" UH2
- 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" UH2
- 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]
- }
+ if (HydroModel != "GR2M") {
+ # Calcul des ordonnees SH2 de l' "hydrogramme unitaire cumule" UH2
+ 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" UH2
+ 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_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
# --------------------------------------------------------------------------------
@@ -308,10 +308,12 @@ if (getRversion() >= "2.15.1") {
# NEUTRALISATION DE P
# --------------------------------------------------------------------------------
- # Interception
- segments(x0 = xy_E[1]-50, x1 = xy_P[1]+50,
- 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)
+ if (HydroModel != "GR2M") {
+ # Interception
+ segments(x0 = xy_E[1]-50, x1 = xy_P[1]+50,
+ 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
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") {
segments(x0 = xy_E[1], x1 = xy_E[1],
y0 = y_interception, y1 = y_rendement)
+ if (HydroModel != "GR2M") {
# Ecriture
plotrix::boxed.labels(x = xy_P[1], y = y_interception, labels = "Pn",
bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
plotrix::boxed.labels(x = xy_E[1], y = y_interception, labels = "En",
bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
+ }
# ETP
if (OutputsModel$PotEvap[i_pdt] != 0) {
@@ -361,11 +365,11 @@ if (getRversion() >= "2.15.1") {
}
# Pn et Ps
-
- points(x = x_Ps, y = y_rendement+1.2*tmp_decal,
- type = "p", pch = tri_B, col = col_P,
- cex = cex_tri(OutputsModel$Ps[i_pdt], fact = fact_triangle, max = cex_max_poly))
-
+ if (HydroModel != "GR2M") {
+ points(x = x_Ps, y = y_rendement+1.2*tmp_decal,
+ type = "p", pch = tri_B, col = col_P,
+ 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,
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))
@@ -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))
}
- # Remplissage de UH1
- PR_mat_UH1_lg <- ceiling(Param[4])
- PR_mat_UH1_id <- max(i_pdt-PR_mat_UH1_lg+1, 1):i_pdt
- PR_mat_UH1 <- matrix(rep(c(rep(0, times = PR_mat_UH1_lg-length(PR_mat_UH1_id)+1),
- OutputsModel$PR[PR_mat_UH1_id]), times = PR_mat_UH1_lg),
- ncol = PR_mat_UH1_lg+1)[, -1L]
- PR_mat_UH1[lower.tri(PR_mat_UH1)] <- 0
+ if (HydroModel != "GR2M") {
+ # Remplissage de UH1
+ PR_mat_UH1_lg <- ceiling(Param[4])
+ PR_mat_UH1_id <- max(i_pdt-PR_mat_UH1_lg+1, 1):i_pdt
+ PR_mat_UH1 <- matrix(rep(c(rep(0, times = PR_mat_UH1_lg-length(PR_mat_UH1_id)+1),
+ OutputsModel$PR[PR_mat_UH1_id]), times = PR_mat_UH1_lg),
+ 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") {
cex = cex_tri(OutputsModel$PR[i_pdt]*0.1, fact = fact_triangle, max = cex_max_poly))
}
-
- # Remplissage de UH2
- 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 <- 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),
- ncol = PR_mat_UH2_lg+1)[, -1L]
- PR_mat_UH2[lower.tri(PR_mat_UH2)] <- 0
+ if (HydroModel != "GR2M") {
+ # Remplissage de UH2
+ 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 <- 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),
+ ncol = PR_mat_UH2_lg+1)[, -1L]
+ PR_mat_UH2[lower.tri(PR_mat_UH2)] <- 0
+ }
}
@@ -578,49 +585,51 @@ if (getRversion() >= "2.15.1") {
y0 = y_routage, y1 = y_routage)
}
- # Q9
- if (OutputsModel$Q9[i_pdt] != 0) {
- points(x = xy_Q9[1], y = xy_Q9[2]+tmp_decal,
- type = "p", pch = tri_B, col = col_P,
- 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,
+ if (HydroModel != "GR2M") {
+ # Q9
+ if (OutputsModel$Q9[i_pdt] != 0) {
+ points(x = xy_Q9[1], y = xy_Q9[2]+tmp_decal,
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,
+ 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,
+ 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,
- 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)
+ 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)
+
}
- 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
# --------------------------------------------------------------------------------
@@ -629,6 +638,9 @@ if (getRversion() >= "2.15.1") {
tmp_triche <- 0#80
# Reservoir de routage
+ if (HydroModel == "GR2M") {
+ Param[3] <- 600
+ }
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 = col_SR, border = NA)
@@ -646,26 +658,29 @@ if (getRversion() >= "2.15.1") {
segments(x0 = xy_min_ROUT[1]+base_res/2, x1 = xy_Q[1],
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 (OutputsModel$QR[i_pdt] != 0) {
+ if (HydroModel != "GR2M") {
+ # Qr
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))
+ 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 (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") {
# EXCHANGE
# --------------------------------------------------------------------------------
+ if (HydroModel != "GR2M") {
# Actual exchange Q9
arrows(x0 = xy_min_ROUT[1]+base_res, x1 = 1025,
y0 = y_Ech_Q9 , y1 = y_Ech_Q9,
@@ -754,6 +770,18 @@ if (getRversion() >= "2.15.1") {
points(x = xy_Q1[1]+100, y = y_Ech_Q1,
type = "p", pch = pch, col = col_P,
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.
if (HydroModel%in% c("GR6J")) {
--
GitLab