Utils.R 30.2 KB
Newer Older
1
## =================================================================================
2
## commands to avoid warnings during package checking when global variables are used
3
## =================================================================================
4
5

if (getRversion() >= "2.15.1") {
6
7
  utils::globalVariables(c(".ShinyGR.args"))
  utils::suppressForeignCheck(c(".ShinyGR.args"))
8
9
  utils::globalVariables(c(".ShinyGR.hist"))
  utils::suppressForeignCheck(c(".ShinyGR.hist"))
10
11
}

12
13


14

15
## =================================================================================
16
## function to manage the model units
17
18
## =================================================================================

19
.TypeModelGR <- function(x) {
20
  
21
22
  if (!is.list(x)) {
    x <- list(TypeModel = x)
23
  }
24
  if (any(class(x) %in% c("PrepGR", "CalGR", "SimGR")) || names(x) %in% "TypeModel") {
25
    x <- x$TypeModel
26
27
  }
  
28
  StrName    <- "(.*)(GR)(\\d{1})(\\D{1})"
29
  NameModel  <- gsub(StrName, "\\2\\3\\4", x)
30
  TimeUnitFR <- gsub(StrName, "\\4", x)
31
32
33
34
  # TimeUnit   <- switch(TimeUnitFR, H = "hour", J = "day", M = "month", A = "year")
  # TimeLag    <- switch(TimeUnit, "hour" = 3600, "day" = 3600*24, "month" = 3600*24*31, "year" = 366)
  TimeUnit   <- switch(TimeUnitFR, H = "h", J = "d", M = "month", A = "y")
  TimeLag    <- switch(TimeUnit, "h" = 3600, "d" = 3600*24, "month" = 3600*24*31, "y" = 366)
35
36
  NbParam    <- gsub(StrName, "\\3", x)
  isCN       <- grepl("CemaNeige"  , x)
37
  
38
39
40
  res <- list(TypeModel = x, NameModel = NameModel, CemaNeige = isCN, 
              NbParam = as.numeric(NbParam),
              TimeUnit = TimeUnit, TimeLag = TimeLag)
41
  return(res)
42
43
44
}


45
46


47
## =================================================================================
48
## function to plot the gr models diagrams (only GR4J and GR5J)
49
## =================================================================================
50

51
.DiagramGR <- function(OutputsModel, Param, SimPer, EventDate, HydroModel, CemaNeige) {
52
53
  
  
54
  # --------------------------------------------------------------------------------
55
  # PARAMETRES
56
  # --------------------------------------------------------------------------------
57
58
  
  # Parametres
59
  mgp             <- c(0, 0.75, 0)
60
  col_P           <- rgb(066, 139, 202, maxColorValue = 255) #"royalblue"
61
  col_E           <- rgb(164, 196, 000, maxColorValue = 255) #"forestgreen"
62
  col_Q           <- "orangered"
63
64
  col_SP          <- adjustcolor("cyan4"   , alpha.f = 0.60)
  col_SR          <- adjustcolor("darkblue", alpha.f = 0.60)
65
66
67
  col_R           <- rgb(066, 139, 202, maxColorValue = 255) #rgb(037, 155, 210, maxColorValue = 255)
  col_mod_bg      <- rgb(245, 245, 245, maxColorValue = 255)
  col_mod_bd      <- rgb(231, 231, 231, maxColorValue = 255)
68
  xy_E            <- c(250, 980)
69
70
  xy_PE           <- c(250, 940)
  xy_AE           <- c(250, 860)
71
  xy_P            <- c(600, 980)
72
  xy_Precip       <- c(600, 950)
73
  xy_Q            <- c(700,  30)
74
75
76
77
78
  x_Ps            <- 440
  x_PnPs          <- 700
  y_interception  <- 900
  y_rendement     <- 815
  y_percolation   <- 575
79
80
  xy_Q9           <- c(400, 310)
  xy_Q1           <- c(800, 310)
81
  y_routage       <- 100
82
  fact_res        <- 200/800/3
83
  fact_resExp     <- 1
84
85
86
87
  base_res        <- 300
  NH              <- 10
  xy_min_PROD     <- c(200, 610)
  xy_min_ROUT     <- c(250, 150)
88
  xy_min_EXPO     <- c(200, 250)
89
  y_entreeUH      <- 500
90
91
92
93
  xy_UH1          <- c(500, 420)
  xy_UH2          <- c(900, 420)
  y_Ech_Q1        <- 170 #200
  y_Ech_Q9        <- 150 #180
94
95
96
  D               <- 5/2
  xpad            <- 1.5
  ypad            <- 1.5
97
98
  max_triangle    <- max(unlist(OutputsModel[c("Perc", "PR", "Q9", "Q1", "QR", "QD", "Pn", "Ps", 
                                               "AE", "Precip", "PotEvap", "AExch1", "AExch2")]))
99
  fact_var        <- 40
100
  fact_triangle   <- 100#25
101
102
  cex_max_poly    <- 2 # 0.005
  cex_tri         <- function(cex, fact = 25, max) suppressWarnings(log(abs(cex) * fact + 1) / max)
103
104
  radius1         <- 0
  radius2         <- 60
105
106
107
  tri_R           <- -0x25BA
  tri_B           <- -0x25BC
  tri_L           <- -0x25C4
108
  tri_T           <- -0x25B2
109
  
110
  par(col.axis = par("fg"), cex.axis = 1.3, cex.lab = 1.3, cex = 0.7, mgp = mgp)
111
  
112
  if (.GlobalEnv$.ShinyGR.args$theme == "Cyborg") {
113
114
    col_mod_bg    <- rgb(255-245, 255-245, 255-245, maxColorValue = 255)
    col_mod_bd    <- rgb(255-231, 255-231, 255-231, maxColorValue = 255)
115
  }
116
  if (.GlobalEnv$.ShinyGR.args$theme == "Flatly") {
117
118
    col_mod_bg    <- "#ECF0F1"
    col_mod_bd    <- "#ECF0F1"
119
  }
120
  
121
  # Pas de temps
122
123
124
  dates_deb       <- EventDate
  n_pdt           <- length(which(OutputsModel$DatesR >= EventDate & OutputsModel$DatesR <= SimPer[2L]))
  i_pdt           <- which(format(OutputsModel$DatesR, "%Y%m%d") == format(EventDate, "%Y%m%d"))
125
126
127
  
  
  
128
  # --------------------------------------------------------------------------------
129
  # UH 1 & 2
130
  # --------------------------------------------------------------------------------
131
  
132
  if (HydroModel %in% c("GR4J", "GR6J")) {
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
    # Calcul des ordonnees SH1 de l' "hydrogramme unitaire cumule" UH1
    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" UH1
    UH1     <- array(NA, NH)
    for (j in 1:NH) {
      if (j == 1) {
        UH1[j] <- SH1[j]
      } else {
        UH1[j] <- SH1[j] - SH1[j-1]
      }  
    }
    
    # Parametres
    max_UH1 <- log(sqrt(max(max(UH1)*OutputsModel$PR*0.9))+1)
    
154
  }
155
156
  
  
157
  # Calcul des ordonnees SH2 de l' "hydrogramme unitaire cumule" UH2
158
159
  SH2     <- array(NA, 2*NH)
  for (i in 1:(2*NH)) {
160
161
    if (i <= 0)                           SH2[i] <- 0
    if (i > 0 & i < Param[4])             SH2[i] <- 0.5*(i/Param[4])^(D)
162
    if (i >= Param[4] & i < 2*Param[4])   SH2[i] <- 1 - (0.5*(2-i/Param[4])^(D))
163
    if (i >= 2*Param[4])                  SH2[i] <- 1
164
165
  }
  
166
  # Calcul des ordonnees UH2 de l' "hydrogramme unitaire discret" UH2
167
168
169
  UH2     <- array(NA, 2*NH)
  for (j in 1:(2*NH)) {
    if (j == 1) {
170
      UH2[j] <- SH2[j]
171
    } else {
172
      UH2[j] <- SH2[j] - SH2[j-1]
173
174
175
176
    }  
  }
  
  # Parametres
177
  max_UH2 <- log(sqrt(max(max(UH2)*OutputsModel$PR*0.1))+1) 
178
179
  
  
180
181
182
  # --------------------------------------------------------------------------------
  # PARTITIONNEMENT FENETRE GRAPHIQUE
  # --------------------------------------------------------------------------------
183
  
184
  # layout(matrix(c(1:4, 4, 4), nrow = 3, ncol = 2, byrow = FALSE), widths = c(1.0, 0.6))
185
  
186
  # --------------------------------------------------------------------------------
187
  # PLUIE ET ETP
188
  # --------------------------------------------------------------------------------
189
190
  
  # P
191
192
193
194
195
  # par(mar = c(2, 4, 1, 1), mgp = mgp)
  # plot(OutputsModel$Dates, OutputsModel$Precip, type = "h", col = col_P, ylim = rev(range(OutputsModel$Precip)), xaxt = "n", ylab = "precip. [mm/d]")
  # rect(xleft = EventDate, ybottom = par("usr")[3], xright =  par("usr")[2], ytop =  par("usr")[4], col = adjustcolor(par("bg"), alpha.f = 0.75), border = NA)
  # abline(v = EventDate, col = "grey", lwd = 2, lty = 2)
  # box()
196
197
  
  # ETP
198
199
200
201
202
  # par(mar = c(2, 4, 1, 1), mgp = mgp)
  # plot(OutputsModel$Dates, OutputsModel$PotEvap, pch = 19, col = col_E, xaxt = "n", ylab = "evapo. [mm/d]")
  # rect(xleft = EventDate, ybottom = par("usr")[3], xright =  par("usr")[2], ytop =  par("usr")[4], col = adjustcolor(par("bg"), alpha.f = 0.75), border = NA)
  # abline(v = EventDate, col = "grey", lwd = 2, lty = 2)
  # box()  
203
204
  
  
205
206
207
  # --------------------------------------------------------------------------------
  # DEBIT
  # --------------------------------------------------------------------------------
208
  
209
  # Q
210
211
212
213
214
215
  # par(mar = c(2, 4, 1, 1), mgp = mgp)
  # plot(OutputsModel$Dates, OutputsModel$Qobs, type = "l", ylab = "flow [mm/d]")
  # lines(OutputsModel$Dates[1:i_pdt], OutputsModel$Qsim[1:i_pdt], type = "l", col = "orangered")
  # rect(xleft = EventDate, ybottom = par("usr")[3], xright =  par("usr")[2], ytop =  par("usr")[4], col = adjustcolor(par("bg"), alpha.f = 0.75), border = NA)
  # abline(v = EventDate, col = "grey", lwd = 2, lty = 2)
  # box()  
216
  
217
  
218
  # --------------------------------------------------------------------------------
219
  # SCHEMAS MODELES
220
  # --------------------------------------------------------------------------------
221
222
  
  # Cadre
223
  par(mar = rep(0.2, 4))
224
  par(fg = par("fg"))
225
  plot(x = 0, type = "n", xlab = "", ylab = "", axes = FALSE, ylim = c(0, 1000), xlim = c(0, 1000))
226
  
227
  # Le modele
228
  rect(xleft = 0, xright = 1000, ybottom = 50, ytop = 970, col = col_mod_bg, border = col_mod_bd)
229
230
  
  
231
  # --------------------------------------------------------------------------------
232
  # ENTREES / SORTIES
233
  # --------------------------------------------------------------------------------
234
  
235
  # Entrees P et ETP
236
237
238
  if (CemaNeige) {
    text(x = xy_P[1]*1.65, y = xy_P[2]*0.98, labels = "+ CemaNeige", adj = c(1, 1), font = 2, col = "grey40", cex = 1.6)
  }
239
240
  text(x = xy_P[1], y = xy_P[2], labels = "P", pos = 3, font = 2, col = col_P, cex = 1.8)
  text(x = xy_E[1], y = xy_E[2], labels = "E", pos = 3, font = 2, col = col_E, cex = 1.8)
241
242
  
  # Sorties Q
243
  text(x = xy_Q[1], y = xy_Q[2], labels = "Q", pos = 1, font = 2, col = col_Q, cex = 1.8)
244
  
245
246
  # Parametres
  tmp_decal   <- 20
247
  
248
  # --------------------------------------------------------------------------------
249
  # NEUTRALISATION DE P
250
  # --------------------------------------------------------------------------------
251
252
  
  # Interception
253
254
255
  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)
256
257
  
  # P vers Pn
258
  segments(x0 = xy_P[1], x1 = xy_P[1], y0 = xy_P[2], y1 = y_interception+tmp_decal)
259
260
  
  # Pn vers Ps
261
262
263
264
265
266
267
268
269
270
271
272
273
  segments(x0 = xy_P[1], x1 = xy_P[1],
           y0 = y_interception, y1 = y_rendement+2*tmp_decal)
  segments(x0 = x_Ps, x1 = xy_P[1],
           y0 = y_rendement+2*tmp_decal, y1 = y_rendement+2*tmp_decal)
  segments(x0 = x_Ps, x1 = x_Ps,
           y0 = y_rendement+2*tmp_decal, y1 = y_rendement)
  
  # Pn vers Pn - Ps
  segments(x0 = xy_P[1], x1 = x_PnPs,
           y0 = y_rendement+2*tmp_decal, y1 = y_rendement+2*tmp_decal)
  segments(x0 = x_PnPs , x1 = x_PnPs,
           y0 = y_rendement+2*tmp_decal, y1 = y_rendement)
  
274
  # Pn - Ps vers Pr
275
276
  segments(x0 = x_PnPs, x1 = x_PnPs,
           y0 = y_rendement, y1 = y_percolation)
277
278
  
  # E vers En puis Es
279
280
281
282
  segments(x0 = xy_E[1], x1 = xy_E[1],
           y0 = xy_E[2], y1 = y_interception+tmp_decal)
  segments(x0 = xy_E[1], x1 = xy_E[1],
           y0 = y_interception, y1 = y_rendement)
283
284
  
  # Ecriture
285
286
287
288
  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)
289
  
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
  # ETP
  if (OutputsModel$PotEvap[i_pdt] != 0) {
    points(x = xy_PE[1], y =  xy_PE[2],
           type = "p", pch = tri_T, col = col_E,
           cex = cex_tri(OutputsModel$PotEvap[i_pdt], fact = fact_triangle, max = cex_max_poly))
  }
  
  # Precipitation
  if (OutputsModel$Precip[i_pdt] != 0) {
    points(x = xy_Precip[1], y =  xy_Precip[2],
           type = "p", pch = tri_B, col = col_P,
           cex = cex_tri(OutputsModel$Precip[i_pdt], fact = fact_triangle, max = cex_max_poly))
  }
  
  # 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))
  
  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))
  
314
  
315
  # --------------------------------------------------------------------------------
316
  # FONCTION DE RENDEMENT
317
  # --------------------------------------------------------------------------------
318
319
  
  # Es
320
321
  plotrix::boxed.labels(x = xy_E[1], y = y_rendement, labels = "Es",
                        bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
322
  
323
324
325
326
327
328
329
  # Evaporation reelle
  if (OutputsModel$AE[i_pdt] != 0) {
    points(x = xy_AE[1], y =  xy_AE[2],
           type = "p", pch = tri_T, col = col_P,
           cex = cex_tri(OutputsModel$AE[i_pdt], fact = fact_triangle, max = cex_max_poly))
  }
  
330
  # Ps et Pn - Ps
331
332
333
334
  plotrix::boxed.labels(x = x_Ps  , y = y_rendement, labels = "Ps"   ,
                        bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
  plotrix::boxed.labels(x = x_PnPs, y = y_rendement, labels = "Pn - Ps",
                        bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
335
  
336
  # Reservoir de production
337
338
  rect(xleft   = xy_min_PROD[1], xright = xy_min_PROD[1]+base_res,
       ybottom = xy_min_PROD[2], ytop   = xy_min_PROD[2]+OutputsModel$Prod[i_pdt]*fact_res,
339
       col = col_SP, border = NA)
340
341
342
343
344
345
  segments(x0 = xy_min_PROD[1], x1 = xy_min_PROD[1]+base_res,
           y0 = xy_min_PROD[2], y1 = xy_min_PROD[2])
  segments(x0 = xy_min_PROD[1], x1 = xy_min_PROD[1],
           y0 = xy_min_PROD[2], y1 = xy_min_PROD[2]+Param[1]*fact_res)
  segments(x0 = xy_min_PROD[1]+base_res, x1 = xy_min_PROD[1]+base_res,
           y0 = xy_min_PROD[2], y1 = xy_min_PROD[2]+Param[1]*fact_res)
346
  text(x = 30, y = xy_min_PROD[2]+15, labels = "Prod.\nstore", cex = 1.4, pos = 4)
347
348
  
  
349
  # --------------------------------------------------------------------------------
350
  # PERCOLATION
351
  # --------------------------------------------------------------------------------
352
  
353
  # Reservoir de production vers Pr
354
355
356
357
  segments(x0 = xy_min_PROD[1]+base_res/2, x1 = xy_min_PROD[1]+base_res/2,
           y0 = xy_min_PROD[2], y1 = y_percolation)
  segments(x0 = xy_min_PROD[1]+base_res/2, x1 = x_PnPs,
           y0 = y_percolation, y1 = y_percolation)
358
359
  
  # Perc
360
361
  plotrix::boxed.labels(x = xy_min_PROD[1]+base_res/2, y = y_percolation, labels = "Perc.",
                        bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
362
363
  
  # Valeur de Perc
364
  if (OutputsModel$Perc[i_pdt] != 0) {
365
366
367
    points(x = xy_min_PROD[1]+base_res+75, y = y_percolation,
           type = "p", pch = tri_R, col = col_P,
           cex = cex_tri(OutputsModel$Perc[i_pdt], fact = fact_triangle, max = cex_max_poly))
368
369
  }
  
370
  # parametres
371
  tmp_decal   <- (y_percolation - y_entreeUH) / 2 
372
  
373
  # Pr vers UH
374
  segments(x0 = x_PnPs, x1 = x_PnPs,
375
           y0 = y_percolation, y1 = y_entreeUH+tmp_decal/2)
376
  
377
  if (HydroModel %in% c("GR4J", "GR6J")) {
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
    
    # --------------------------------------------------------------------------------
    # SEPARATION DE PR
    # --------------------------------------------------------------------------------
    
    # Pr vers UH1
    segments(x0 = xy_Q9[1], x1 = x_PnPs,
             y0 = y_entreeUH+tmp_decal/2, y1 = y_entreeUH+tmp_decal/2)
    segments(x0 = xy_Q9[1], x1 = xy_Q9[1],
             y0 = y_entreeUH+tmp_decal/2, y1 = xy_Q9[2])
    
    # Pr vers UH2
    segments(x0 = x_PnPs, x1 = xy_Q1[1],
             y0 = y_entreeUH+tmp_decal/2, y1 = y_entreeUH+tmp_decal/2)
    segments(x0 = xy_Q1[1], x1 = xy_Q1[1],
             y0 = y_entreeUH+tmp_decal/2, y1 = y_routage)
    
    # Pr
396
397
    plotrix::boxed.labels(x = x_PnPs, y = y_percolation, labels = "Pr",
                          bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
    
    # Pr
    if (OutputsModel$PR[i_pdt] != 0) {
      points(x = x_PnPs[1], y = y_entreeUH+tmp_decal,
             type = "p", pch = tri_B, col = col_P,
             cex = cex_tri(OutputsModel$PR[i_pdt], fact = fact_triangle, max = cex_max_poly))
    }
    
    # --------------------------------------------------------------------------------
    # HYDROGRAMME UNITAIRE 1
    # --------------------------------------------------------------------------------
    
    # Entree de UH1
    if (OutputsModel$PR[i_pdt] != 0) {
      points(x = xy_Q9[1], y =y_entreeUH,
             type = "p", pch = tri_B, col = col_P,
             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
    
    # --------------------------------------------------------------------------------
    # HYDROGRAMME UNITAIRE 2
    # --------------------------------------------------------------------------------
    
    # Entree de UH2
    if (OutputsModel$PR[i_pdt] != 0) {
      points(x = xy_Q1[1], y = y_entreeUH,
             type = "p", pch = tri_B, col = col_P,
             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
    
  } 
446
  
447
  if (HydroModel == "GR5J") {
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
    
    # --------------------------------------------------------------------------------
    # SEPARATION DE PR
    # --------------------------------------------------------------------------------
    
    # sortie UH
    segments(x0 = x_PnPs, x1 = x_PnPs,
             y0 = y_entreeUH-2*tmp_decal, y1 = y_entreeUH-3*tmp_decal)
    
    # sortie UH vers branche 1
    segments(x0 = xy_Q9[1], x1 = x_PnPs,
             y0 = y_entreeUH-3*tmp_decal, y1 = y_entreeUH-3*tmp_decal)
    segments(x0 = xy_Q9[1], x1 = xy_Q9[1],
             y0 = y_entreeUH-3*tmp_decal, y1 = xy_Q9[2])
    
    # sortie UH vers branche 2
    segments(x0 = x_PnPs, x1 = xy_Q1[1],
             y0 = y_entreeUH-3*tmp_decal, y1 = y_entreeUH-3*tmp_decal)
    segments(x0 = xy_Q1[1], x1 = xy_Q1[1],
             y0 = y_entreeUH-3*tmp_decal, y1 = y_routage)
    
    # Pr
470
471
    plotrix::boxed.labels(x = x_PnPs, y = y_percolation, labels = "Pr",
                          bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
    
    
    # --------------------------------------------------------------------------------
    # HYDROGRAMME UNITAIRE
    # --------------------------------------------------------------------------------
    
    # Entree de UH (PR)
    if (OutputsModel$PR[i_pdt] != 0) {
      points(x = x_PnPs[1], y = y_entreeUH+tmp_decal,
             type = "p", pch = tri_B, col = col_P,
             cex = cex_tri(OutputsModel$PR[i_pdt], 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
    
    
    # Sorties de UH
    if (PR_mat_UH2[1] != 0) {
      points(x = x_PnPs[1], y =  y_entreeUH-5*tmp_decal/2,
             type = "p", pch = tri_B, col = col_P,
             cex = cex_tri(PR_mat_UH2[1], fact = fact_triangle, max = cex_max_poly))
    }
    
  }   
503
  
504
505
506
507
508
509
  # sortie de UH1 vers reservoirs exponentiel et de routage
  if (HydroModel == "GR6J") {
    segments(x0 = xy_Q9[1], x1 = xy_Q9[1],
             y0 = y_entreeUH-3*tmp_decal, y1 = xy_Q9[2])
    segments(x0 = xy_Q9[1]*0.80, x1 = xy_Q9[1]*1.30,
             y0 = xy_Q9[2], y1 = xy_Q9[2])
510
511
    segments(x0 = xy_Q9[1]*1.30, x1 = xy_Q9[1]*1.30,
             y0 = xy_Q9[2], y1 = xy_Q9[2]*0.65)    
512
513
    segments(x0 = xy_Q9[1]*0.80, x1 = xy_Q9[1]*0.80,
             y0 = xy_Q9[2], y1 = xy_Q9[2]*0.90)
514
515
516
517
    segments(x0 = xy_Q9[1]*0.55, x1 = xy_Q9[1]*0.55,
             y0 = xy_Q9[2]*0.70, y1 = y_routage)
    segments(x0 = xy_Q9[1]*0.55, x1 = xy_min_ROUT[1]+base_res/2,
             y0 = y_routage, y1 = y_routage)
518
519
  }
  
520
  # Q9
521
522
523
524
  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))
525
526
527
528
529
530
531
532
533
534
535
536
537
    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)
    }
  }
538
  
539
540
  plotrix::boxed.labels(x = xy_Q9[1], y = xy_Q9[2], labels = "Q9",
                        bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
541
  
542
  # Q1
543
  if (OutputsModel$Q1[i_pdt] != 0) {
544
545
546
547
    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)
548
  }
549
  
550
  plotrix::boxed.labels(x = xy_Q1[1], y = xy_Q1[2], labels = "Q1", bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
551
  
552
  # Valeur de Qd
553
  if (OutputsModel$QD[i_pdt] != 0) {
554
555
556
    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))
557
  }
558
  
559
  # Qd
560
  plotrix::boxed.labels(x = xy_Q1[1], y = y_routage, labels = "Qd", bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
561
  
562
  
563
  # --------------------------------------------------------------------------------
564
  # RESERVOIR DE ROUTAGE
565
  # --------------------------------------------------------------------------------
566
  
567
  # Triche pour la taille du reservoire de routage
568
  tmp_triche   <- 0#80
569
  
570
  # Reservoir de routage
571
572
  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,
573
       col = col_SR, border = NA)
574
575
576
577
578
579
  segments(x0 = xy_min_ROUT[1], x1 = xy_min_ROUT[1]+base_res, 
           y0 = xy_min_ROUT[2], y1 = xy_min_ROUT[2])
  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)
  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)
580
  text(x = 30, y = xy_min_ROUT[2]+15, labels = "Routing\nstore", cex = 1.4, pos = 4)
581
  
582
  # Sorties du reservoir
583
584
585
586
  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)
  segments(x0 = xy_min_ROUT[1]+base_res/2, x1 = xy_Q[1],
           y0 = y_routage, y1 = y_routage)
587
  
588
  # Qr
589
590
591
592
593
594
595
596
  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)
  }
597
  
598
  # Valeur de Qr
599
  if (OutputsModel$QR[i_pdt] != 0) {
600
    if (HydroModel != "GR6J") {
601
602
603
    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))
604
605
606
607
608
    } 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))
    }
609
  }
610
  
611

612
  
613
614
615
616
617
618
619
620
621
  # --------------------------------------------------------------------------------
  # RESERVOIR EXPONENTIEL
  # --------------------------------------------------------------------------------

  if (HydroModel == "GR6J") {
    
    # Triche pour la taille du reservoire exponentiel
    tmp_triche   <- 0#80
    
622
623
624
625
626
    # Exp en log
    signExp <- ifelse(OutputsModel$Exp[i_pdt] > 0, +1, -1)
    logExpIpdt <- log(abs(OutputsModel$Exp[i_pdt])+1e-6) * signExp
    logExpMax  <- log(max(abs(OutputsModel$Exp  ))+1e-6)

627
628
    # Reservoir exponentiel
    rect(xleft = xy_min_EXPO[1], xright = xy_min_EXPO[1]+base_res,
629
         ybottom = xy_min_EXPO[2], ytop = xy_min_EXPO[2]+logExpIpdt*fact_resExp+tmp_triche,
630
         col = ifelse(OutputsModel$Exp[i_pdt] > 0, "#10B510", "#FF0303"), border = NA)
631
632
633
    segments(x0 = xy_min_EXPO[1], x1 = xy_min_EXPO[1]+base_res, 
             y0 = xy_min_EXPO[2], y1 = xy_min_EXPO[2])
    segments(x0 = xy_min_EXPO[1], x1 = xy_min_EXPO[1],
634
             y0 = xy_min_EXPO[2], y1 = xy_min_EXPO[2]+logExpMax*fact_resExp+tmp_triche)
635
    segments(x0 = xy_min_EXPO[1]+base_res, x1 = xy_min_EXPO[1]+base_res,
636
             y0 = xy_min_EXPO[2], y1 = xy_min_EXPO[2]+logExpMax*fact_resExp+tmp_triche)
637
    segments(x0 = xy_min_EXPO[1], x1 = xy_min_EXPO[1],
638
             y0 = xy_min_EXPO[2], y1 = xy_min_EXPO[2]-logExpMax*fact_resExp-tmp_triche)
639
    segments(x0 = xy_min_EXPO[1]+base_res, x1 = xy_min_EXPO[1]+base_res,
640
             y0 = xy_min_EXPO[2], y1 = xy_min_EXPO[2]-logExpMax*fact_resExp-tmp_triche)
641
642
    text(x = 30, y = xy_min_EXPO[2]+00, labels = "Exp.\nstore", cex = 1.4, pos = 4)
    points(x = 180, y = xy_min_EXPO[2]+20, pch = 43, # +
643
         cex = 2.0, col = "#10B510")
644
    points(x = 178, y = xy_min_EXPO[2]-20, pch = 95, # -
645
646
647
648
649
650
651
652
         cex = 1.6, col = "#FF0303")

    # Valeur de QrExp
    if (OutputsModel$QR[i_pdt] != 0) {
      points(x = xy_Q[1]-350, y = y_routage,
             type = "p", pch = tri_R, col = col_P,
             cex = cex_tri(OutputsModel$QRExp[i_pdt], fact = fact_triangle, max = cex_max_poly))
    }
653

654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
    # Valeur de Qr + QrExp
    if (OutputsModel$QR[i_pdt] != 0) {    
      points(x = xy_Q[1]-100, y = y_routage,
             type = "p", pch = tri_R, col = col_P,
             cex = cex_tri(OutputsModel$QR[i_pdt]+OutputsModel$QRExp[i_pdt], fact = fact_triangle, max = cex_max_poly))
    }
  }
  
  
  # --------------------------------------------------------------------------------
  # Q FINAL
  # -------------------------------------------------------------------------------
  
  # Q final
  segments(x0 = xy_Q[1], x1 = xy_Q[1], y0 = y_routage, y1 = xy_Q[2]+10)
  if (OutputsModel$Qsim[i_pdt] != 0) {
    points(x = xy_Q[1], y =  mean(c(y_routage, xy_Q[2]+10)),
           type = "p", pch = tri_B, col = col_Q,
           cex = cex_tri(OutputsModel$Qsim[i_pdt], fact = fact_triangle, max = cex_max_poly))
673
674
675
  }
  
  
676
677
678
679
680
681
682
683
  # --------------------------------------------------------------------------------
  # EXCHANGE
  # --------------------------------------------------------------------------------
  
  # Actual exchange Q9
  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)
684
685
  pch <- ifelse(OutputsModel$AExch1[i_pdt] < 0, tri_R, tri_L)
  points(x = xy_min_ROUT[1]+base_res+130, y = y_Ech_Q9,
686
687
         type = "p", pch = pch, col = col_P,
         cex = cex_tri(OutputsModel$AExch1[i_pdt], fact = fact_triangle, max = cex_max_poly))
688
  
689
690
691
692
  # Actual exchange Q1
  arrows(x0 = xy_Q1[1], x1 = 1025,
         y0 = y_Ech_Q1, y1 = y_Ech_Q1,
         length = 0.075, angle = 20, code = 3)
693
  pch <- ifelse(OutputsModel$AExch2[i_pdt] < 0, tri_R, tri_L)
694
  points(x = xy_Q1[1]+100, y = y_Ech_Q1,
695
         type = "p", pch = pch, col = col_P,
696
697
698
699
700
701
702
703
704
705
706
707
         cex = cex_tri(OutputsModel$AExch2[i_pdt], fact = fact_triangle, max = cex_max_poly)) 
  
  # Actual exchange Q9 exp.
  if (HydroModel%in% c("GR6J")) {
    arrows(x0 = xy_min_EXPO[1]+base_res[1], x1 = 1025,
           y0 = xy_min_EXPO[2], y1 = xy_min_EXPO[2],
           length = 0.075, angle = 20, code = 3)
    pch <- ifelse(OutputsModel$Exch[i_pdt] < 0, tri_R, tri_L)
    points(x = xy_Q1[1]+100, y = xy_min_EXPO[2],
           type = "p", pch = pch, col = col_P,
           cex = cex_tri(OutputsModel$Exch[i_pdt], fact = fact_triangle, max = cex_max_poly)) 
  }
708
  
709
  if (HydroModel%in% c("GR4J", "GR6J")) {
710
711
712
713
    
    # --------------------------------------------------------------------------------
    # UH 1 & 2 PLOT
    # --------------------------------------------------------------------------------
714
715
    
    tmpUH1 <- PR_mat_UH1 %*% UH1[seq_len(PR_mat_UH1_lg)] * 0.75 # 0.75 au lieu de 0.90 pour reduire dif. visuelle
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
    tmpUH2 <- PR_mat_UH2 %*% UH2[seq_len(PR_mat_UH2_lg)] * 0.25 # 0.25 au lieu de 0.10 pour reduire dif. visuelle
    palUH0 <- colorRampPalette(c("#428BCA", "#B5D2EA"))
    par(mar = rep(1, 4))
    par(new = TRUE, plt = c(0.35, 0.46,  0.40, 0.45), yaxt = "n", bg = "red")
    rect(xleft = 325, xright = 475,  ybottom = 390, ytop = 445, col = col_mod_bg, border = NA)
    par(new = TRUE)
    barplot(tmpUH1, beside = TRUE, space = 0, col = palUH0(length(tmpUH1)),
            ylim = range(c(0, 5/length(PR_mat_UH2_lg)), na.rm = TRUE))
    par(new = TRUE, plt = c(0.67, 0.89,  0.40, 0.45), yaxt = "n")
    rect(xleft = 0, xright = 900,  ybottom = 0, ytop = 470, col = col_mod_bg, border = NA)
    par(new = TRUE)
    barplot(tmpUH2, beside = TRUE, space = 0, col = palUH0(length(tmpUH2)),
            ylim = range(c(0, 5/length(PR_mat_UH2_lg)), na.rm = TRUE))
    
  }
  
  
733
  if (HydroModel == "GR5J") {
734
735
736
737
738
739
740
741
742
743
744
745
746
747
    
    # --------------------------------------------------------------------------------
    # UH PLOT
    # --------------------------------------------------------------------------------
    
    tmpUH2 <- PR_mat_UH2 %*% UH2[seq_len(PR_mat_UH2_lg)] * 1 / 2 # * 1/2 car seul HU et pour que ca rentre
    palUH0 <- colorRampPalette(c("#428BCA", "#B5D2EA"))
    par(mar = rep(1, 4))
    par(new = TRUE, plt = c(0.57, 0.79,  0.45, 0.5), yaxt = "n")
    par(new = TRUE)
    barplot(tmpUH2, beside = TRUE, space = 0, col = palUH0(length(tmpUH2)),
            ylim = range(c(0, 5/length(PR_mat_UH2_lg)), na.rm = TRUE))
    
  }
748
  
749
750
}

751