feat: Add model LRQ

NAMESPACE

@@ -56,6 +56,7 @@ export(RunModel_GR6J)
 export(RunModel_Lag)
 export(RunModel_LLR)
 export(RunModel_LBLR)
+export(RunModel_LRQ)
 export(SeriesAggreg)
 export(TransfoParam)
 export(TransfoParam_CemaNeige)

R/RunModel_LRQ.R

+RunModel_LRQ <- function(InputsModel, RunOptions, Param, QcontribDown) {
+  NParam <- 2
+  #RunModel_LRQ
+  ## argument check
+  if (!inherits(InputsModel, "InputsModel")) {
+    stop("'InputsModel' must be of class 'InputsModel'")
+  }
+  if (!inherits(InputsModel, "SD")) {
+    stop("'InputsModel' must be of class 'SD'")
+  }
+  if (!inherits(RunOptions, "RunOptions")) {
+    stop("'RunOptions' must be of class 'RunOptions'")
+  }
+  if (!is.vector(Param) | !is.numeric(Param)) {
+    stop("'Param' must be a numeric vector")
+  }
+  if (sum(!is.na(Param)) != NParam) {
+    stop(paste("'Param' must be a vector of length", NParam, "and contain no NA"))
+  }
+  if (inherits(QcontribDown, "OutputsModel")) {
+    if (is.null(QcontribDown$Qsim)) {
+      stop("'QcontribDown' should contain a key 'Qsim' containing the output of the runoff of the downstream subcatchment")
+    }
+    if (length(QcontribDown$Qsim) != length(RunOptions$IndPeriod_Run)) {
+      stop("Time series Qsim in 'QcontribDown' should have the same length as 'RunOptions$IndPeriod_Run'")
+    }
+    if (!identical(RunOptions$IndPeriod_WarmUp, 0L) && !identical(RunOptions$Outputs_Sim, RunOptions$Outputs_Cal)) {
+      # This test is not necessary during calibration but usefull in other cases because
+      # WarmUpQsim is then used for downstream sub-basins because of the delay in Qupstream
+      if (is.null(QcontribDown$RunOptions$WarmUpQsim) ||
+          length(QcontribDown$RunOptions$WarmUpQsim) != length(RunOptions$IndPeriod_WarmUp)) {
+        stop("Time series WarmUpQsim in 'QcontribDown' should have the same length as 'RunOptions$IndPeriod_WarmUp'")
+      }
+    }
+  } else if (is.vector(QcontribDown) && is.numeric(QcontribDown)) {
+    if (length(QcontribDown) != length(RunOptions$IndPeriod_Run)) {
+      stop("'QcontribDown' should have the same length as 'RunOptions$IndPeriod_Run'")
+    }
+  } else {
+    stop("'QcontribDown' must be a numeric vector or a 'OutputsModel' object")
+  }
+
+  # data set up
+
+  NbUpBasins <- length(InputsModel$LengthHydro)
+  if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
+    RunOptions$IndPeriod_WarmUp <- NULL
+  }
+  IndPeriod1   <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
+  LInputSeries <- as.integer(length(IndPeriod1))
+  IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
+
+  if (inherits(QcontribDown, "OutputsModel")) {
+    OutputsModel <- QcontribDown
+    if (is.null(OutputsModel$RunOptions$WarmUpQsim)) {
+      OutputsModel$RunOptions$WarmUpQsim <- rep(NA, length(RunOptions$IndPeriod_WarmUp))
+    }
+    QsimDown <- c(OutputsModel$RunOptions$WarmUpQsim, OutputsModel$Qsim)
+  } else if (is.vector(QcontribDown) && is.numeric(QcontribDown)) {
+    OutputsModel <- list()
+    class(OutputsModel) <- c("OutputsModel", class(RunOptions)[-1])
+    QsimDown <- c(rep(NA, length(RunOptions$IndPeriod_WarmUp)),
+                  QcontribDown)
+  }
+
+
+  ## Parameters set up
+  TParam <- Param[1L]
+  KParam <- Param[2L]
+
+  PT <-  InputsModel$LengthHydro * 1e3 / TParam / RunOptions$FeatFUN_MOD$TimeStep
+  PK <- KParam * PT
+
+
+  ## Lag set up
+
+  Lag <-  function(t, Qroute, Qupstream){
+    Q <-  c(Qupstream[1:(1 + floor(t))], Qroute[1:(length(Qroute) - (floor(t) + 1))])
+    return(Q)
+  }
+  ## set up initial states
+  IniSD <- RunOptions$IniStates[grep("SD", names(RunOptions$IniStates))]
+  if (length(IniSD) > 0) {
+    if (sum(floor(PT)) + NbUpBasins != length(IniSD)) {
+      stop(
+        sprintf(
+          "SD initial states has a length of %i and a length of %i is required",
+          length(IniSD),
+          sum(floor(PT)) + NbUpBasins
+        )
+      )
+    }
+    IniStates <- lapply(seq_len(NbUpBasins), function(x) {
+      iStart <- 1
+      if (x > 1) {
+        iStart <- iStart + sum(floor(PT[1:x - 1]) + 1)
+      }
+      as.vector(IniSD[iStart:(iStart + PT[x])])
+    })
+  } else {
+    IniStates <- lapply(
+      seq_len(NbUpBasins),
+      function(iUpBasins) {
+        iWarmUp <- seq.int(
+          from = max(1, IndPeriod1[1] - floor(PT[iUpBasins]) - 1),
+          to   = max(1, IndPeriod1[1] - 1)
+        )
+        ini <- InputsModel$Qupstream[iWarmUp, iUpBasins]
+        if (length(ini) != floor(PT[iUpBasins] + 1)) {
+          # If warm-up period is not enough long complete beginning with first value
+          ini <- c(rep(ini[1], floor(PT[iUpBasins] + 1) - length(ini)), ini)
+        }
+        return(as.vector(ini))
+      }
+    )
+  }
+  ## Lag & route model computation
+  Qsim_m3 <- matrix(0
+                    , nrow = length(IndPeriod1), ncol = NbUpBasins)
+  QsimDown_input <- matrix(QsimDown *
+                       InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3, ncol = 1)
+
+
+
+  for (upstream_basin in seq_len(NbUpBasins)) {
+    Qupstream_m3 <- c(IniStates[[upstream_basin]],
+                      InputsModel$Qupstream[IndPeriod1, upstream_basin])
+    Qroute <- Qupstream_m3
+
+    for (q in seq_along(Qroute)[1:(length(Qroute) - 1)]){
+      Qroute[q + 1] <- ((Qroute[q]^0.5 + (((Qupstream_m3[q] + Qupstream_m3[q + 1])/2)/PK[upstream_basin])) /
+                          (1 + ((Qroute[q]^0.5) / PK[upstream_basin])))^2
+    }
+
+    Qllr <- Lag(PT[upstream_basin], Qroute, Qupstream_m3)[IndPeriod1]
+    Qsim_m3[, upstream_basin] <- Qllr
+  }
+
+  Qsim_m3 <- cbind(QsimDown_input, Qsim_m3)
+  Qsim_m3 <- rowSums(Qsim_m3)
+
+  ## OutputsModel
+
+  if ("Qsim_m3" %in% RunOptions$Outputs_Sim) {
+    OutputsModel$Qsim_m3 <- Qsim_m3[IndPeriod2]
+  }
+
+  if ("Qsim" %in% RunOptions$Outputs_Sim) {
+    # Convert back Qsim to mm
+    OutputsModel$Qsim <- Qsim_m3[IndPeriod2] / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
+    # message("Qsim: ", paste(OutputsModel$Qsim, collapse = ", "))
+  }
+
+  if ("QsimDown" %in% RunOptions$Outputs_Sim) {
+    # Convert back Qsim to mm
+    OutputsModel$QsimDown <- QsimDown[IndPeriod2]
+  }
+
+  # Warning for negative flows or NAs only in extended outputs
+  if (length(RunOptions$Outputs_Sim) > 2) {
+    if (any(OutputsModel$Qsim[!is.na(OutputsModel$Qsim)] < 0)) {
+      warning(length(which(OutputsModel$Qsim < 0)), " time steps with negative flow, set to zero.")
+      OutputsModel$Qsim[OutputsModel$Qsim < 0] <- 0
+    }
+    # Warning for NAs
+    if (any(is.na(OutputsModel$Qsim))) {
+      warning(length(which(is.na(OutputsModel$Qsim))), " time steps with NA values")
+    }
+  }
+
+  if ("StateEnd" %in% RunOptions$Outputs_Sim) {
+    SD <- lapply(seq(NbUpBasins), function(x) {
+      lastTS <- RunOptions$IndPeriod_Run[length(RunOptions$IndPeriod_Run)]
+      InputsModel$Qupstream[(lastTS - floor(PT[x])):lastTS, x]
+    })
+    if (is.null(OutputsModel$StateEnd)) {
+      OutputsModel$StateEnd <- CreateIniStates(RunModel_LLR, InputsModel, SD = SD)
+    } else {
+      OutputsModel$StateEnd$SD <- SD
+    }
+    # message("StateEnd: ", paste(OutputsModel$StateEnd$SD, collapse = ", "))
+  }
+  if ("WarmUpQsim" %in% RunOptions$Outputs_Sim) {
+    OutputsModel$RunOptions$WarmUpQsim <- Qsim_m3[seq_len(length(RunOptions$IndPeriod_WarmUp))] / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
+  }
+
+  if ("Param" %in% RunOptions$Outputs_Sim) {
+    OutputsModel$RunOptions$Param <- c(Param, OutputsModel$RunOptions$Param)
+  }
+
+  class(OutputsModel) <- c(class(OutputsModel), "SD")
+
+  return(OutputsModel)
+}

R/TransfoParam_LRQ.R

+TransfoParam_LRQ <- function(ParamIn, Direction) {
+
+  ## number of model parameters
+  NParam <- 2L
+
+  ## check arguments
+  isVecParamIn <- is.vector(ParamIn)
+  if (isVecParamIn) {
+    ParamIn <- matrix(ParamIn, nrow = 1)
+  }
+  if (!inherits(ParamIn, "matrix")) {
+    stop("'ParamIn' must be of class 'matrix'")
+  }
+  if (!inherits(Direction, "character") | length(Direction) != 1 | any(!Direction %in% c("RT", "TR"))) {
+    stop("'Direction' must be a character vector of length 1 equal to 'RT' or 'TR'")
+  }
+  if (ncol(ParamIn) != NParam) {
+    stop(sprintf("the Lag & Route model requires %i parameters", NParam))
+  }
+
+
+  ## transformation
+  if (Direction == "TR") {
+    ParamOut <- ParamIn
+    ParamOut[, 1] <- 20 * (ParamIn[, 1] + 10) / 20.0
+    ParamOut[, 2] <- 20 * (ParamIn[, 2] + 10) / 20.0
+  }
+  if (Direction == "RT") {
+    ParamOut <- ParamIn
+    ParamOut[, 1] <- ParamIn[, 1] * 20.0 / 20 - 10
+    ParamOut[, 2] <- ParamIn[, 2] * 20.0 / 20 - 10
+  }
+
+
+  ## check output
+  if (isVecParamIn) {
+    ParamOut <- as.vector(ParamOut)
+  }
+
+
+  return(ParamOut)
+
+}

R/Utils.R

@@ -83,6 +83,9 @@
       }else if(identical(RunModel_LBLR, FUN_SD)){
         res$NParamSD <- FeatMod$NbParam[16]
         res$FUN_SD <- "LBLR"
+      }else if (identical(RunModel_LRQ, FUN_SD)){
+        res$NParamSD <- FeatMod$NbParam[17]
+        res$FUN_SD <- "LRQ"
       }
     }
     return(res)

R/UtilsCalibOptions.R

@@ -33,7 +33,10 @@
       FUN_ROUT <- TransfoParam_LLR
     } else if (identical(RunModel_LBLR, FUN_SD)){
       FUN_ROUT <- TransfoParam_LBLR
+    } else if (identical(RunModel_LRQ, FUN_SD)){
+      FUN_ROUT <- TransfoParam_LRQ
     }
+
   }
 
     ## set FUN_TRANSFO

inst/modelsFeatures/FeatModelsGR.csv

@@ -15,3 +15,4 @@ CemaNeigeGR5H;CemaNeigeGR5H;7;hourly;NA;GR;airGR
 Lag;Lag;1;NA;NA;SD;airGR
 LLR;LLR;2;NA;NA;SD;airGR
 LBLR;LBLR;2;NA;NA;SD;airGR
+LRQ;LRQ;2;NA;NA;SD;airGR