Utils.R 32.5 KB
Newer Older
1
.onAttach <- function(libname, pkgname) {
2
  if (packageVersion("htmlwidgets") <= "1.5.2") {
3
    base::packageStartupMessage("\n---------------------------\n")
4
    base::packageStartupMessage("This version of 'airGRteaching' is designed to work with 'htmlwidgets' >= 1.5.2.9000 (troubles with 'dygraphs')")
5
6
7
8
9
10
    base::packageStartupMessage("Install the latest version of 'htmlwidgets' from GitHub with the following command lines:")  
    base::packageStartupMessage("\tinstall.packages(\"remotes\")\n\tremotes::install_github(\"ramnathv/htmlwidgets\")")
    base::packageStartupMessage("\n---------------------------\n")
  }
}

11
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
30
31
32
33
34
35
36
37
38
39
40
41
## =================================================================================
## function to test if a remote url exists
## =================================================================================

.CheckUrl <- function(url, timeout = 2) {
  con <- url(description = url)
  check <- suppressWarnings(try(open.connection(con = con, open = "rt", timeout = t), silent =  TRUE)[1])
  suppressWarnings(try(close.connection(con), silent = TRUE))
  is.null(check)
}




42
43
44
45
## =================================================================================
## function to compute the start and stop id of equal values in a vector
## =================================================================================

46
.StartStop <- function(x, FUN) {
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
  naQ_rle <- rle(FUN(x))
  naQ_ide <- cumsum(naQ_rle$lengths)[naQ_rle$values]   + 1
  naQ_ids <- naQ_ide - naQ_rle$lengths[naQ_rle$values] - 1
  idNA   <- data.frame(start = naQ_ids, stop = naQ_ide)
  idNA$start <- ifelse(idNA$start < 1         , 1       , idNA$start)
  idNA$stop  <- ifelse(idNA$stop  > length(x), length(x), idNA$stop )
  idNA
}




## =================================================================================
## function for drawing several shadows of dygraphic regions simultaneously
## =================================================================================

63
.DyShadingMulti <- function(dygraph, ts, idStart, IdStop, ...) {
64
65
66
67
68
69
70
71
72
73
  for (i in seq_along(idStart)) {
    dygraph <- dygraphs::dyShading(dygraph = dygraph,
                                   from    = as.character(ts)[idStart[i]],
                                   to      = as.character(ts)[IdStop[i]],
                                   ...)
  }
  dygraph
}


74
## =================================================================================
75
## function to manage the model units
76
77
## =================================================================================

78
.TypeModelGR <- function(x) {
79
  
80
81
  if (!is.list(x)) {
    x <- list(TypeModel = x)
82
  }
83
  if (any(class(x) %in% c("PrepGR", "CalGR", "SimGR")) || names(x) %in% "TypeModel") {
84
    x <- x$TypeModel
85
86
  }
  
87
  StrName    <- "(.*)(GR)(\\d{1})(\\D{1})"
88
  NameModel  <- gsub(StrName, "\\2\\3\\4", x)
89
  TimeUnitFR <- gsub(StrName, "\\4", x)
90
91
92
93
  # 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)
94
95
  NbParam    <- gsub(StrName, "\\3", x)
  isCN       <- grepl("CemaNeige"  , x)
96
  
97
98
99
  res <- list(TypeModel = x, NameModel = NameModel, CemaNeige = isCN, 
              NbParam = as.numeric(NbParam),
              TimeUnit = TimeUnit, TimeLag = TimeLag)
100
  return(res)
101
102
103
}


104
105


106
## =================================================================================
107
## function to plot the gr models diagrams (only GR4J and GR5J)
108
## =================================================================================
109

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

  if (HydroModel == "GR6J") {
    
    # Triche pour la taille du reservoire exponentiel
    tmp_triche   <- 0#80
    
682
683
684
685
686
    # 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)

687
688
    # Reservoir exponentiel
    rect(xleft = xy_min_EXPO[1], xright = xy_min_EXPO[1]+base_res,
689
         ybottom = xy_min_EXPO[2], ytop = xy_min_EXPO[2]+logExpIpdt*fact_resExp+tmp_triche,
690
         col = ifelse(OutputsModel$Exp[i_pdt] > 0, "#10B510", "#FF0303"), border = NA)
691
692
693
    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],
694
             y0 = xy_min_EXPO[2], y1 = xy_min_EXPO[2]+logExpMax*fact_resExp+tmp_triche)
695
    segments(x0 = xy_min_EXPO[1]+base_res, x1 = xy_min_EXPO[1]+base_res,
696
             y0 = xy_min_EXPO[2], y1 = xy_min_EXPO[2]+logExpMax*fact_resExp+tmp_triche)
697
    segments(x0 = xy_min_EXPO[1], x1 = xy_min_EXPO[1],
698
             y0 = xy_min_EXPO[2], y1 = xy_min_EXPO[2]-logExpMax*fact_resExp-tmp_triche)
699
    segments(x0 = xy_min_EXPO[1]+base_res, x1 = xy_min_EXPO[1]+base_res,
700
             y0 = xy_min_EXPO[2], y1 = xy_min_EXPO[2]-logExpMax*fact_resExp-tmp_triche)
701
702
    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, # +
703
         cex = 2.0, col = "#10B510")
704
    points(x = 178, y = xy_min_EXPO[2]-20, pch = 95, # -
705
706
707
708
709
710
711
712
         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))
    }
713

714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
    # 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))
733
734
735
  }
  
  
736
737
738
739
740
741
742
743
  # --------------------------------------------------------------------------------
  # 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)
744
745
  pch <- ifelse(OutputsModel$AExch1[i_pdt] < 0, tri_R, tri_L)
  points(x = xy_min_ROUT[1]+base_res+130, y = y_Ech_Q9,
746
747
         type = "p", pch = pch, col = col_P,
         cex = cex_tri(OutputsModel$AExch1[i_pdt], fact = fact_triangle, max = cex_max_poly))
748
  
749
750
751
752
  # 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)
753
  pch <- ifelse(OutputsModel$AExch2[i_pdt] < 0, tri_R, tri_L)
754
  points(x = xy_Q1[1]+100, y = y_Ech_Q1,
755
         type = "p", pch = pch, col = col_P,
756
757
758
759
760
761
762
763
764
765
766
767
         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)) 
  }
768
  
769
  if (HydroModel%in% c("GR4J", "GR6J")) {
770
771
772
773
    
    # --------------------------------------------------------------------------------
    # UH 1 & 2 PLOT
    # --------------------------------------------------------------------------------
774
775
    
    tmpUH1 <- PR_mat_UH1 %*% UH1[seq_len(PR_mat_UH1_lg)] * 0.75 # 0.75 au lieu de 0.90 pour reduire dif. visuelle
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
    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))
    
  }
  
  
793
  if (HydroModel == "GR5J") {
794
795
796
797
798
799
800
801
802
803
804
805
806
807
    
    # --------------------------------------------------------------------------------
    # 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))
    
  }
808
  
809
810
}

811