Utils.R 30.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
## =================================================================================
## function to generate a message to reinstall the dygraph package
## =================================================================================

.onAttach <- function(lib, pkg) {
  packageStartupMessage("\nThe airGRteaching needs the last version of the 'dygraphs' package.\n",
                        "Please, reinstall it with the following command:\n",
                        "\tdevtools::install_github(c(\"ramnathv/htmlwidgets\", \"rstudio/dygraphs\"), force = TRUE)")
}


12
13


14
## =================================================================================
15
## commands to avoid warnings during package checking when global variables are used
16
## =================================================================================
17
18

if (getRversion() >= "2.15.1") {
19
20
  utils::globalVariables(c(".ShinyGR.args"))
  utils::suppressForeignCheck(c(".ShinyGR.args"))
21
22
  utils::globalVariables(c(".ShinyGR.hist"))
  utils::suppressForeignCheck(c(".ShinyGR.hist"))
23
24
}

25
26


27

28
## =================================================================================
29
## function to manage the model units
30
31
## =================================================================================

32
.TypeModelGR <- function(x) {
33
  
34
35
  if (!is.list(x)) {
    x <- list(TypeModel = x)
36
  }
37
38
  if (any(class(x) %in% c("ObsGR", "CalGR", "SimGR")) || names(x) %in% "TypeModel") {
    x <- x$TypeModel
39
40
  }
  
41
  StrName    <- "(.*)(GR)(\\d{1})(\\D{1})"
42
  NameModel  <- gsub(StrName, "\\2\\3\\4", x)
43
  TimeUnitFR <- gsub(StrName, "\\4", x)
44
45
46
47
  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)
  NbParam    <- gsub(StrName, "\\3", x)
  isCN       <- grepl("CemaNeige"  , x)
48
  
49
50
51
  res <- list(TypeModel = x, NameModel = NameModel, CemaNeige = isCN, 
              NbParam = as.numeric(NbParam),
              TimeUnit = TimeUnit, TimeLag = TimeLag)
52
  return(res)
53
54
55
}


56
57


58
## =================================================================================
59
## function to plot the gr models diagrams (only GR4J and GR5J)
60
## =================================================================================
61

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

619
  
620
621
622
623
624
625
626
627
628
629
630
631
  # --------------------------------------------------------------------------------
  # RESERVOIR EXPONENTIEL
  # --------------------------------------------------------------------------------

  if (HydroModel == "GR6J") {
    
    # Triche pour la taille du reservoire exponentiel
    tmp_triche   <- 0#80
    
    # Reservoir exponentiel
    rect(xleft = xy_min_EXPO[1], xright = xy_min_EXPO[1]+base_res,
         ybottom = xy_min_EXPO[2], ytop = xy_min_EXPO[2]+OutputsModel$Exp[i_pdt]*fact_res+tmp_triche,
632
         col = ifelse(OutputsModel$Exp[i_pdt] > 0, "#10B510", "#FF0303"), border = NA)
633
634
635
636
637
638
639
640
641
642
643
644
645
    # rect(xleft = xy_min_EXPO[1], xright = xy_min_EXPO[1]+base_res,
    #      ybottom = xy_min_EXPO[2], ytop = xy_min_EXPO[2]-OutputsModel$Exp[i_pdt]*fact_res-tmp_triche,
    #      col = col_SR, border = NA)
    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],
             y0 = xy_min_EXPO[2], y1 = xy_min_EXPO[2]+max(abs(OutputsModel$Exp))*fact_res+tmp_triche)
    segments(x0 = xy_min_EXPO[1]+base_res, x1 = xy_min_EXPO[1]+base_res,
             y0 = xy_min_EXPO[2], y1 = xy_min_EXPO[2]+max(abs(OutputsModel$Exp))*fact_res+tmp_triche)
    segments(x0 = xy_min_EXPO[1], x1 = xy_min_EXPO[1],
             y0 = xy_min_EXPO[2], y1 = xy_min_EXPO[2]-max(abs(OutputsModel$Exp))*fact_res-tmp_triche)
    segments(x0 = xy_min_EXPO[1]+base_res, x1 = xy_min_EXPO[1]+base_res,
             y0 = xy_min_EXPO[2], y1 = xy_min_EXPO[2]-max(abs(OutputsModel$Exp))*fact_res-tmp_triche)
646
647
    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, # +
648
         cex = 2.0, col = "#10B510")
649
    points(x = 178, y = xy_min_EXPO[2]-20, pch = 95, # -
650
651
652
653
654
655
656
657
         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))
    }
658

659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
    # 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))
678
679
680
  }
  
  
681
682
683
684
685
686
687
688
  # --------------------------------------------------------------------------------
  # 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)
689
  pch = ifelse(OutputsModel$AExch1[i_pdt] < 0, tri_R, tri_L)
690
  points(x = xy_min_ROUT[1]+base_res+100, y = y_Ech_Q9,
691
692
         type = "p", pch = pch, col = col_P,
         cex = cex_tri(OutputsModel$AExch1[i_pdt], fact = fact_triangle, max = cex_max_poly))
693
  
694
695
696
697
  # 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)
698
  pch = ifelse(OutputsModel$AExch2[i_pdt] < 0, tri_R, tri_L)
699
  points(x = xy_Q1[1]+100, y = y_Ech_Q1,
700
701
         type = "p", pch = pch, col = col_P,
         cex = cex_tri(OutputsModel$AExch2[i_pdt], fact = fact_triangle, max = cex_max_poly))  
702
  
703
  if (HydroModel%in% c("GR4J", "GR6J")) {
704
705
706
707
    
    # --------------------------------------------------------------------------------
    # UH 1 & 2 PLOT
    # --------------------------------------------------------------------------------
708
709
    
    tmpUH1 <- PR_mat_UH1 %*% UH1[seq_len(PR_mat_UH1_lg)] * 0.75 # 0.75 au lieu de 0.90 pour reduire dif. visuelle
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
    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))
    
  }
  
  
727
  if (HydroModel == "GR5J") {
728
729
730
731
732
733
734
735
736
737
738
739
740
741
    
    # --------------------------------------------------------------------------------
    # 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))
    
  }
742
  
743
744
}

745