Utils.R 30.2 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]+Param[1]*fact_res/3, labels = "Prod.\nstore",
354
       cex = 1.4, pos = 4)
355
356
  
  
357
  # --------------------------------------------------------------------------------
358
  # PERCOLATION
359
  # --------------------------------------------------------------------------------
360
  
361
  # Reservoir de production vers Pr
362
363
364
365
  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)
366
367
  
  # Perc
368
369
  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)
370
371
  
  # Valeur de Perc
372
  if (OutputsModel$Perc[i_pdt] != 0) {
373
374
375
    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))
376
377
  }
  
378
  # parametres
379
  tmp_decal   <- (y_percolation - y_entreeUH) / 2 
380
  
381
  # Pr vers UH
382
  segments(x0 = x_PnPs, x1 = x_PnPs,
383
           y0 = y_percolation, y1 = y_entreeUH+tmp_decal/2)
384
  
385
  if (HydroModel %in% c("GR4J", "GR6J")) {
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
    
    # --------------------------------------------------------------------------------
    # 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
404
405
    plotrix::boxed.labels(x = x_PnPs, y = y_percolation, labels = "Pr",
                          bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
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
453
    
    # 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
    
  } 
454
  
455
  if (HydroModel == "GR5J") {
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
    
    # --------------------------------------------------------------------------------
    # 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
478
479
    plotrix::boxed.labels(x = x_PnPs, y = y_percolation, labels = "Pr",
                          bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
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
510
    
    
    # --------------------------------------------------------------------------------
    # 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))
    }
    
  }   
511
  
512
513
514
515
516
517
  # 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])
518
519
    segments(x0 = xy_Q9[1]*1.30, x1 = xy_Q9[1]*1.30,
             y0 = xy_Q9[2], y1 = xy_Q9[2]*0.65)    
520
521
    segments(x0 = xy_Q9[1]*0.80, x1 = xy_Q9[1]*0.80,
             y0 = xy_Q9[2], y1 = xy_Q9[2]*0.90)
522
523
524
525
    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)
526
527
  }
  
528
  # Q9
529
530
531
532
  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))
533
534
535
536
537
538
539
540
541
542
543
544
545
    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)
    }
  }
546
  
547
548
  plotrix::boxed.labels(x = xy_Q9[1], y = xy_Q9[2], labels = "Q9",
                        bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
549
  
550
  # Q1
551
  if (OutputsModel$Q1[i_pdt] != 0) {
552
553
554
555
    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)
556
  }
557
  
558
  plotrix::boxed.labels(x = xy_Q1[1], y = xy_Q1[2], labels = "Q1", bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
559
  
560
  # Valeur de Qd
561
  if (OutputsModel$QD[i_pdt] != 0) {
562
563
564
    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))
565
  }
566
  
567
  # Qd
568
  plotrix::boxed.labels(x = xy_Q1[1], y = y_routage, labels = "Qd", bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
569
  
570
  
571
  # --------------------------------------------------------------------------------
572
  # RESERVOIR DE ROUTAGE
573
  # --------------------------------------------------------------------------------
574
  
575
  # Triche pour la taille du reservoire de routage
576
  tmp_triche   <- 0#80
577
  
578
  # Reservoir de routage
579
580
  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,
581
       col = col_SR, border = NA)
582
583
584
585
586
587
  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)
588
  text(x = 50, y = xy_min_ROUT[2]+Param[3]*fact_res/3, labels = "Routing\nstore",
589
       cex = 1.4, pos = 4)
590
  
591
  # Sorties du reservoir
592
593
594
595
  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)
596
  
597
  # Qr
598
599
600
601
602
603
604
605
  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)
  }
606
  
607
  # Valeur de Qr
608
  if (OutputsModel$QR[i_pdt] != 0) {
609
    if (HydroModel != "GR6J") {
610
611
612
    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))
613
614
615
616
617
    } 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))
    }
618
  }
619
  
620

621
  
622
623
624
625
626
627
628
629
630
631
632
633
  # --------------------------------------------------------------------------------
  # 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,
634
         col = ifelse(OutputsModel$Exp[i_pdt] > 0, "#10B510", "#FF0303"), border = NA)
635
636
637
638
639
640
641
642
643
644
645
646
647
    # 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)
648
    text(x = 50, y = xy_min_EXPO[2]+Param[3]*fact_res/3, labels = "Exp.\nstore",
649
         cex = 1.4, pos = 4)
650
651
652
653
654
655
656
657
658
659
660
    points(x = 180, y = xy_min_EXPO[2]+Param[3]*fact_res/3+05, pch = 43, # +
         cex = 2.0, col = "#10B510")
    points(x = 178, y = xy_min_EXPO[2]+Param[3]*fact_res/3-35, pch = 95, # -
         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))
    }
661

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

748