plot.CalGR.R 4.49 KB
Newer Older
1
plot.CalGR <- function(x, xlab = NULL, ylab =  NULL, main = NULL, which = c("perf", "iter", "ts"), ...) {
2

3
  if (!inherits(x, "CalGR")) {
unknown's avatar
unknown committed
4
5
    stop("Non convenient data for x argument. Must be of class \"CalGR\"")
  }
6
7
  if (! any(which %in% c("perf", "iter", "ts"))) {
    stop("Non convenient data for which argument. Must be of class \"perf\", \"iter\" or \"ts\"")
8
  }
9

10
11
  nbParamX <- .TypeModelGR(x)$NbParam #as.numeric(gsub("\\D", "", x$TypeModel))
  nbParamC <- ifelse(.TypeModelGR(x)$CemaNeige, 2, 0)
unknown's avatar
unknown committed
12
13
  nbParam <- nbParamX + nbParamC
  nmParam <- c(sprintf("X%i", 1:nbParamX), sprintf("C%i", seq_len(nbParamC)))
14

15
16
  opar <- par(no.readonly = TRUE)
  on.exit(par(opar))
17

18
  if (any(which[1L] %in% c("perf"))) {
19
    plot(x$OutputsModel, Qobs = x$Qobs, ...)
20
  }
21

22
  if (any(which[1L] %in% c("iter"))) {
unknown's avatar
unknown committed
23
24
25
26
27
28
29
30
31
32
33
34
35
    layout.list <- list(matrix(c(1:2), ncol = 2),
                        matrix(c(1:3, 3), ncol = 2),
                        matrix(NA),
                        matrix(c(1:5, 5), ncol = 3),
                        matrix(c(1:5, 0, 6, 6), ncol = 4),
                        matrix(c(1:7, 7), ncol = 4),
                        matrix(c(1:5, 0, 6:8, 8), ncol = 5),
                        matrix(c(1:9, 9), ncol = 5))
    layout.mat <- layout.list[[nbParam]]
    layout.w <- c(rep(1, ncol(layout.mat)-1), 2)
    layout.h <- rep(1, nrow(layout.mat))
    layout(mat = layout.mat, widths = layout.w, heights = layout.h)
  }
36
37
38
39
  # if (which[1L] %in%  c("both")) {
  #   layout.mat <- rbind(layout.mat, max(layout.mat)+1, max(layout.mat)+2)
  #   layout(mat = layout.mat, widths = layout.w, heights = layout.h)
  # }
40
  if (which[1L] %in% c("ts")) {
unknown's avatar
unknown committed
41
    layout(mat = matrix(1:2), widths = c(1, 2), heights = c(1, 2))
42
  }
43

44
  if (any(which[1L] %in% c("iter"))) {
45
46
    ParamLab <- data.frame(Name  = c(sprintf("X%i", 1:6), sprintf("C%i", 1:2)),
                           Label = c("prod. store capacity [mm]",
47
                                     "intercatchment exch. coef. [mm/TimeUnit]",
48
                                     "routing store capacity [mm]",
49
50
51
52
                                     "UH time constant [TimeUnit]",
                                     "intercatchment exch. threshold [-]",
                                     "coef. for emptying exp. store [mm]",
                                     "weight for snowpack thermal state [-]",
53
                                     "degree-day melt coef. [mm/degC/TimeUnit]"))
54
    ParamLab$Label <- gsub("TimeUnit", substr(.TypeModelGR(x)$TimeUnit, 1, 1), ParamLab$Label)
55
56


57
    par(mar = c(2.5, 3.5, 2.0, 1.0), mgp = c(2.0, 0.6, 0.0), oma = c(0, 0, 4, 0))
unknown's avatar
unknown committed
58
    for (i in seq_len(ncol(x$OutputsCalib$HistParamR))) {
59
60
      plot(x$OutputsCalib$HistParamR[, i],
           col = ifelse(grepl("X", nmParam[i]), "darkorchid3", "darkolivegreen3"),
unknown's avatar
unknown committed
61
           type = "b", pch = 19, cex = 0.8,
62
           xlab = "", ylab = ParamLab[ParamLab$Name == nmParam[i], "Label"],
63
           main = nmParam[i])
unknown's avatar
unknown committed
64
    }
65
66
    plot(x$OutputsCalib$HistCrit,
         col = "tomato",
unknown's avatar
unknown committed
67
         type = "b", pch = 19, cex = 1.0,
68
         xlab = "", ylab = "",
69
         main = gsub(".*_", "", x$CalCrit))
70

71
72
73
    mainIter <- "Evolution of parameters and efficiency criterion\nduring the iterations of the steepest-descent step"
    title(main = ifelse(is.null(main), mainIter, main),
          outer = TRUE, line = 0.6, cex.main = 1.6)
unknown's avatar
unknown committed
74
  }
75

76
  if (any(which[1L] %in%  c("ts"))) {
77
78
79
80
    if (is.null(xlab)) {
      xlab <- "Time"
    }
    if (is.null(ylab)) {
81
      ylab  <- paste0(c("flow [mm/", "precip. [mm/"), .TypeModelGR(x)$TimeUnit, "]")
82
83
84
85
    } else {
      if (length(ylab) < 2) {
        ylab <- c(ylab, "")
      }
86
    }
unknown's avatar
unknown committed
87
88
89
90
    data <- data.frame(DatesR = x$OutputsModel$DatesR,
                       Precip = x$OutputsModel$Precip,
                       Qobs   = x$Qobs,
                       Qsim   = x$OutputsModel$Qsim)
91

unknown's avatar
unknown committed
92
93
    par(mar = c(0.1, 4, 4, 1), xaxt = "n")
    plot(Precip ~ DatesR, data = data, type = "h", col = "royalblue",
94
         xlab = "", ylab = ylab[2L], main = main, ylim = rev(range(data$Precip)))
95

unknown's avatar
unknown committed
96
    par(mar = c(5, 4, 0.1, 1), xaxt = "s")
97
    matplot(data$DatesR, data[, c("Qobs", "Qsim")], type = "l", col = c("black", "orangered"), lty = 1,
98
         xlab = xlab, ylab = ylab[1L], main = "", xaxt = "n")
unknown's avatar
unknown committed
99
    axis.POSIXct(side = 1, x = data$DatesR)
100

unknown's avatar
unknown committed
101
102
    # if (plot.na) {
    #   axis(side = 1, at = as.POSIXct(data$DatesR[is.na(data$Qobs)]),
103
104
    #        labels = FALSE, lwd.ticks = 3, col.ticks = col.na, tck = 0.025, lend = "butt")
    #
unknown's avatar
unknown committed
105
106
107
    #   legend("topright", legend = "NA", pch = 15, col = col.na, bty = "n", cex = 0.8)
    # }
  }
108
  # box()
109
110

}