refactor: index time series before using it

- remove all specific indexations
Refs #112

R/plot.OutputsModel.R

@@ -15,12 +15,18 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
   })
 
 
-
   OutputsModel <- x
+
+  ## index time series
   if (!is.null(IndPeriod_Plot)) {
+    if (length(IndPeriod_Plot) == 0) {
+      IndPeriod_Plot <- seq_along(OutputsModel$DatesR)
+    }
+    IndPeriod_Plot <- seq_along(IndPeriod_Plot)
     OutputsModel <- .IndexOutputsModel(OutputsModel, IndPeriod_Plot)
     Qobs <- Qobs[IndPeriod_Plot]
-    IndPeriod_Plot <- seq_along(IndPeriod_Plot)
+  } else {
+    IndPeriod_Plot <- seq_along(OutputsModel$DatesR)
   }
 
 
@@ -191,9 +197,6 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     #   stop("the time step of the model inputs could not be found")
     # }
   }
-  if (length(IndPeriod_Plot) == 0) {
-    IndPeriod_Plot <- 1:length(OutputsModel$DatesR)
-  }
   if (inherits(OutputsModel, "CemaNeige")) {
     NLayers <- length(OutputsModel$CemaNeigeLayers)
   }
@@ -209,12 +212,12 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
   }
   BOOL_QobsZero <- FALSE
   if (BOOL_Qobs) {
-    SelectQobsNotZero <- round(Qobs[IndPeriod_Plot], 4) != 0
+    SelectQobsNotZero <- round(Qobs, 4) != 0
     BOOL_QobsZero <- sum(!SelectQobsNotZero, na.rm = TRUE) > 0
   }
   BOOL_QsimZero <- FALSE
   if (BOOL_Qsim) {
-    SelectQsimNotZero <- round(OutputsModel$Qsim[IndPeriod_Plot], 4) != 0
+    SelectQsimNotZero <- round(OutputsModel$Qsim, 4) != 0
     BOOL_QsimZero <- sum(!SelectQsimNotZero, na.rm = TRUE) > 0
   }
   if ( BOOL_Qobs & !BOOL_Qsim) {
@@ -310,7 +313,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
 
 
 
-    Xaxis <- as.POSIXct(OutputsModel$DatesR[IndPeriod_Plot])
+    Xaxis <- as.POSIXct(OutputsModel$DatesR)
 
     if (!is.null(BasinArea)) {
       Factor_UNIT_M3S <- switch(NameTS,
@@ -340,7 +343,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     mar <- c(3, 5, 1, 5)
 
     par(new = FALSE, mar = mar)
-    ylim1 <- range(OutputsModel$Precip[IndPeriod_Plot], na.rm = TRUE)
+    ylim1 <- range(OutputsModel$Precip, na.rm = TRUE)
     ylim2 <- ylim1 * c(1.0, 1.1)
     ylim2 <- rev(ylim2)
 
@@ -348,7 +351,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     if (NameTS %in% c("month", "year")) {
       lwdP <- lwd * 2
     }
-    plot(Xaxis, OutputsModel$Precip[IndPeriod_Plot],
+    plot(Xaxis, OutputsModel$Precip,
          type = "h", xaxt = "n", yaxt = "n", yaxs = "i", ylim = ylim2,
          col = "royalblue", lwd = lwdP * lwdk, lend = 1,
          xlab = "", ylab = "", ...)
@@ -359,7 +362,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
       legend("bottomright", legend = c("solid","liquid"),
              col = c("lightblue", "royalblue"), lty = c(1, 1), lwd = c(lwd, lwd),
              bty = "o", bg = bg, box.col = bg, cex = cex.leg)
-      points(Xaxis, PsolLayerMean[IndPeriod_Plot],
+      points(Xaxis, PsolLayerMean,
            type = "h", xaxt = "n", yaxt = "n", yaxs = "i", ylim = ylim2,
            col = "lightblue", lwd = lwdP * lwdk, lend = 1,
            xlab = "", ylab = "", ...)
@@ -378,23 +381,23 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
 
     par(new = FALSE, mar = mar)
     if (!BOOLPLOT_ActuEvap) {
-      ylim1 <- range(OutputsModel$PotEvap[IndPeriod_Plot], na.rm = TRUE)
+      ylim1 <- range(OutputsModel$PotEvap, na.rm = TRUE)
       xlabE <- "pot. evap."
     } else {
-      ylim1 <- range(c(OutputsModel$PotEvap[IndPeriod_Plot],
-                       OutputsModel$AE[IndPeriod_Plot]),
+      ylim1 <- range(c(OutputsModel$PotEvap,
+                       OutputsModel$AE),
                      na.rm = TRUE)
       xlabE <- "evap."
     }
     ylim2 <- ylim1 #* c(1.0, 1.1)
 
 
-    plot(Xaxis, OutputsModel$PotEvap[IndPeriod_Plot],
+    plot(Xaxis, OutputsModel$PotEvap,
          type = "l", xaxt = "n", yaxt = "n", ylim = ylim2,
          col = "green3", lwd = lwd * lwdk,
          xlab = "", ylab = "", ...)
     if (BOOLPLOT_ActuEvap) {
-      lines(Xaxis, OutputsModel$AE[IndPeriod_Plot],
+      lines(Xaxis, OutputsModel$AE,
            type = "l", xaxt = "n", yaxt = "n", ylim = ylim2,
            col = "green4", lwd = lwd * lwdk, lty = 3)
       legend("topright", legend = c("pot.", "actu."), col = c("green3", "green4"),
@@ -416,10 +419,10 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     mar <- c(3, 5, 1, 5)
 
     par(new = FALSE, mar = mar)
-    ylim1 <- range(OutputsModel$AE[IndPeriod_Plot], na.rm = TRUE)
+    ylim1 <- range(OutputsModel$AE, na.rm = TRUE)
     ylim2 <- ylim1 #* c(1.0, 1.1)
 
-    plot(Xaxis, OutputsModel$AE[IndPeriod_Plot],
+    plot(Xaxis, OutputsModel$AE,
          type = "l", xaxt = "n", yaxt = "n", ylim = ylim2,
          col = "green4", lwd = lwd * lwdk,
          xlab = "", ylab = "", ...)
@@ -448,12 +451,12 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
         SnowPackLayerMean <- SnowPackLayerMean + OutputsModel$CemaNeigeLayers[[iLayer]]$Temp/NLayers
       }
     }
-    plot(Xaxis, SnowPackLayerMean[IndPeriod_Plot], type = "n", ylim = ylim1, xlab = "", ylab = "", xaxt = "n", yaxt = "n", ...)
+    plot(Xaxis, SnowPackLayerMean, type = "n", ylim = ylim1, xlab = "", ylab = "", xaxt = "n", yaxt = "n", ...)
     for (iLayer in 1:NLayers) {
-      lines(Xaxis, OutputsModel$CemaNeigeLayers[[iLayer]]$Temp[IndPeriod_Plot], lty = 3, col = "orchid", lwd = lwd * lwdk * 0.8)
+      lines(Xaxis, OutputsModel$CemaNeigeLayers[[iLayer]]$Temp, lty = 3, col = "orchid", lwd = lwd * lwdk * 0.8)
     }
     abline(h = 0, col = "grey", lty = 2)
-    lines(Xaxis, SnowPackLayerMean[IndPeriod_Plot], type = "l", lwd = lwd * lwdk *1.0, col = "darkorchid4")
+    lines(Xaxis, SnowPackLayerMean, type = "l", lwd = lwd * lwdk *1.0, col = "darkorchid4")
     axis(side = 2, at = pretty(ylim1), labels = pretty(ylim1), cex.axis = cex.axis, ...)
 
     mtext(side = 2, expression(paste("temp. [", degree, "C]"), sep = ""),  padj = 0.2, line = line, cex = cex.lab)
@@ -483,9 +486,9 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
         SnowPackLayerMean <- SnowPackLayerMean + OutputsModel$CemaNeigeLayers[[iLayer]]$SnowPack/NLayers
       }
     }
-    plot(Xaxis, SnowPackLayerMean[IndPeriod_Plot], type = "l", ylim = ylim1, lwd = lwd * lwdk *1.2, col = "royalblue", xlab = "", ylab = "", xaxt = "n", yaxt = "n", ...)
+    plot(Xaxis, SnowPackLayerMean, type = "l", ylim = ylim1, lwd = lwd * lwdk *1.2, col = "royalblue", xlab = "", ylab = "", xaxt = "n", yaxt = "n", ...)
     for (iLayer in 1:NLayers) {
-      lines(Xaxis, OutputsModel$CemaNeigeLayers[[iLayer]]$SnowPack[IndPeriod_Plot], lty = 3, col = "royalblue", lwd = lwd * lwdk *0.8)
+      lines(Xaxis, OutputsModel$CemaNeigeLayers[[iLayer]]$SnowPack, lty = 3, col = "royalblue", lwd = lwd * lwdk *0.8)
     }
     axis(side = 2, at = pretty(ylim1), labels = pretty(ylim1), cex.axis = cex.axis, ...)
 
@@ -515,21 +518,21 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     DATA3[!SelectQsimNotZero] <- mean(OutputsModel$Qsim, na.rm = TRUE) / 10000
     DATA3 <- log(DATA3)
 
-    ylim1 <- range(DATA3[IndPeriod_Plot], na.rm = TRUE)
+    ylim1 <- range(DATA3, na.rm = TRUE)
     if (BOOL_Qobs) {
-      ylim1 <- range(c(ylim1, DATA2[IndPeriod_Plot]), na.rm = TRUE)
+      ylim1 <- range(c(ylim1, DATA2), na.rm = TRUE)
     }
     ylim2 <- c(ylim1[1], 1.1*ylim1[2])
     plot(Xaxis, rep(NA, length(Xaxis)), type = "n", ylim = ylim2, xlab = "", ylab = "", xaxt = "n", yaxt = "n", ...)
     txtleg <- NULL
     colleg <- NULL
     if (BOOL_Qobs) {
-      lines(Xaxis, DATA2[IndPeriod_Plot], lwd = lwd * lwdk , lty = 1, col = par("fg"))
+      lines(Xaxis, DATA2, lwd = lwd * lwdk , lty = 1, col = par("fg"))
       txtleg <- c(txtleg, "observed")
       colleg <- c(colleg, par("fg"))
     }
     if (BOOL_Qsim) {
-      lines(Xaxis, DATA3[IndPeriod_Plot], lwd = lwd * lwdk , lty = 1, col = "orangered")
+      lines(Xaxis, DATA3, lwd = lwd * lwdk , lty = 1, col = "orangered")
       txtleg <- c(txtleg, "simulated")
       colleg <- c(colleg, "orangered")
     }
@@ -549,21 +552,21 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     kPlot <- kPlot + 1
     mar <- c(3, 5, 1, 5)
     par(new = FALSE, mar = mar)
-    ylim1 <- range(OutputsModel$Qsim[IndPeriod_Plot], na.rm = TRUE)
+    ylim1 <- range(OutputsModel$Qsim, na.rm = TRUE)
     if (BOOL_Qobs) {
-      ylim1 <- range(c(ylim1, Qobs[IndPeriod_Plot]), na.rm = TRUE)
+      ylim1 <- range(c(ylim1, Qobs), na.rm = TRUE)
     }
     ylim2 <- c(ylim1[1], 1.1*ylim1[2])
     plot(Xaxis, rep(NA, length(Xaxis)), type = "n", ylim = ylim2, xlab = "", ylab = "", xaxt = "n", yaxt = "n", ...)
     txtleg <- NULL
     colleg <- NULL
     if (BOOL_Qobs) {
-      lines(Xaxis, Qobs[IndPeriod_Plot], lwd = lwd * lwdk , lty = 1, col = par("fg"))
+      lines(Xaxis, Qobs, lwd = lwd * lwdk , lty = 1, col = par("fg"))
       txtleg <- c(txtleg, "observed")
       colleg <- c(colleg, par("fg"))
     }
     if (BOOL_Qsim) {
-      lines(Xaxis, OutputsModel$Qsim[IndPeriod_Plot], lwd = lwd * lwdk , lty = 1, col = "orangered")
+      lines(Xaxis, OutputsModel$Qsim, lwd = lwd * lwdk , lty = 1, col = "orangered")
       txtleg <- c(txtleg, "simulated")
       colleg <- c(colleg, "orangered")
     }
@@ -586,9 +589,9 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     mar <- c(3, 5, 1, 5)
 
     if (log_scale) {
-      errorQ <- log(OutputsModel$Qsim[IndPeriod_Plot]) - log(Qobs[IndPeriod_Plot])
+      errorQ <- log(OutputsModel$Qsim) - log(Qobs)
     } else {
-      errorQ <- OutputsModel$Qsim[IndPeriod_Plot] - Qobs[IndPeriod_Plot]
+      errorQ <- OutputsModel$Qsim - Qobs
     }
     par(new = FALSE, mar = mar)
     ylim1 <- range(errorQ[SelectNotZero], na.rm = TRUE)
@@ -634,18 +637,18 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     } else {
       ## Data_formating_as_table
       DataModel <- as.data.frame(matrix(as.numeric(NA), nrow = length(IndPeriod_Plot), ncol = 5))
-      DataModel[, 1] <- as.numeric(format(OutputsModel$DatesR[IndPeriod_Plot], format = "%Y%m%d%H"))
+      DataModel[, 1] <- as.numeric(format(OutputsModel$DatesR, format = "%Y%m%d%H"))
       if (BOOL_Pobs) {
-        DataModel[, 2] <- OutputsModel$Precip[IndPeriod_Plot]
+        DataModel[, 2] <- OutputsModel$Precip
       }
       if (BOOL_Psol) {
-        DataModel[, 3] <- PsolLayerMean[IndPeriod_Plot]
+        DataModel[, 3] <- PsolLayerMean
       }
       if (BOOL_Qobs) {
-        DataModel[, 4] <- Qobs[IndPeriod_Plot]
+        DataModel[, 4] <- Qobs
       }
       if (BOOL_Qsim) {
-        DataModel[, 5] <- OutputsModel$Qsim[IndPeriod_Plot]
+        DataModel[, 5] <- OutputsModel$Qsim
       }
       colnames(DataModel) <- c("Dates", "Precip", "Psol", "Qobs", "Qsim")
       TxtDatesDataModel <- formatC(DataModel$Dates, format = "d", width = 8, flag = "0")
@@ -791,15 +794,15 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     xlim <- c(0, 1)
     if ( BOOL_Qobs & !BOOL_Qsim) {
       # SelectNotZero <- SelectQobsNotZero
-      ylim <- range(log(Qobs[IndPeriod_Plot][SelectNotZero]), na.rm = TRUE)
+      ylim <- range(log(Qobs[SelectNotZero]), na.rm = TRUE)
     }
     if (!BOOL_Qobs &  BOOL_Qsim) {
       # SelectNotZero <- SelectQsimNotZero
-      ylim <- range(log(OutputsModel$Qsim[IndPeriod_Plot][SelectNotZero]), na.rm = TRUE)
+      ylim <- range(log(OutputsModel$Qsim[SelectNotZero]), na.rm = TRUE)
     }
     if ( BOOL_Qobs &  BOOL_Qsim) {
       # SelectNotZero <- SelectQobsNotZero & SelectQsimNotZero
-      ylim <- range(log(c(Qobs[IndPeriod_Plot][SelectNotZero], OutputsModel$Qsim[IndPeriod_Plot][SelectNotZero])), na.rm = TRUE)
+      ylim <- range(log(c(Qobs[SelectNotZero], OutputsModel$Qsim[SelectNotZero])), na.rm = TRUE)
     }
     SelectNotZero <- ifelse(is.na(SelectNotZero), FALSE, SelectNotZero)
     if (any(SelectNotZero)) {
@@ -814,7 +817,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
       txtleg <- NULL
       colleg <- NULL
       if (BOOL_Qobs) {
-        DATA2 <- log(Qobs[IndPeriod_Plot][SelectNotZero])
+        DATA2 <- log(Qobs[SelectNotZero])
         SeqQuant <- seq(0, 1, by = 1/(length(DATA2)))
         Quant <- as.numeric(quantile(DATA2, SeqQuant, na.rm = TRUE))
         Fn <- ecdf(DATA2)
@@ -827,7 +830,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
         colleg <- c(colleg, par("fg"))
       }
       if (BOOL_Qsim) {
-        DATA2 <- log(OutputsModel$Qsim[IndPeriod_Plot][SelectNotZero])
+        DATA2 <- log(OutputsModel$Qsim[SelectNotZero])
         SeqQuant <- seq(0, 1, by = 1/(length(DATA2)))
         Quant <- as.numeric(quantile(DATA2, SeqQuant, na.rm = TRUE))
         Fn <- ecdf(DATA2)
@@ -862,9 +865,9 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     mar <- c(6, 5, 1, 5)
     par(new = FALSE, mar = mar)
     if (any(SelectNotZero)) {
-      ylim <- log(range(c(Qobs[IndPeriod_Plot][SelectQobsNotZero & SelectQsimNotZero], OutputsModel$Qsim[IndPeriod_Plot][SelectQobsNotZero & SelectQsimNotZero]), na.rm = TRUE))
-      plot(log(Qobs[IndPeriod_Plot][SelectQobsNotZero & SelectQsimNotZero]),
-           log(OutputsModel$Qsim[IndPeriod_Plot][SelectQobsNotZero & SelectQsimNotZero]),
+      ylim <- log(range(c(Qobs[SelectQobsNotZero & SelectQsimNotZero], OutputsModel$Qsim[SelectQobsNotZero & SelectQsimNotZero]), na.rm = TRUE))
+      plot(log(Qobs[SelectQobsNotZero & SelectQsimNotZero]),
+           log(OutputsModel$Qsim[SelectQobsNotZero & SelectQsimNotZero]),
            type = "p", pch = 1, cex = 0.9, col = par("fg"), lwd = lwd,
            xlim = ylim, ylim = ylim, xaxt = "n", yaxt = "n", xlab = "", ylab = "", ...)
       abline(a = 0, b = 1, col = "royalblue", lwd = lwd)