R/convolution.R

In transfR: Transfer of Hydrograph from Gauged to Ungauged Catchments

#' Convolution of net rainfall with unit hydrograph
#'
#' Simulate the discharge by a convolution between the unit hydrograph and the net rainfall.
#' @name convolution
#' @param Rn net rainfall vector or an object of class \code{transfR}
#' @param uh unit hydrograph vector
#' @param Rcol name of the space-time attribute for the net rainfall simulation in the \code{transfR} object
#' @param Qcol name of the space-time attribute for the discharge in the \code{transfR} object
#' @param continuous logical indicating whether simulated values should be omitted when a time step might be influenced
#' by past or future rainfall (according to the time span of the unit hydrograph)
#' @param save_donor logical indicating if additional discharge simulations should be computed using the
#' net rainfall of each individual donor catchment instead of just the weighted average net rainfall. This
#' requires \code{save_donor} to have been TRUE when using \link{mixr}
#' @param verbose logical indicating if information messages should be written to the console
#' @param ... further arguments passed to or from other methods
#' @return An object of the same class as \code{Rn}. If \code{Rn} is a transfR object,
#' it is returned with the computed results added as new attributes.
#' @examples
#' data(Oudon)
#' icatch <- 1
#' uc <- velocity(hl = Oudon$hl[[icatch]])
#' uh <- uh(hl = Oudon$hl[[icatch]], uc = uc, deltat = units::set_units(1,"h"))$prob
#' Rn <- units::set_units(c(1,5,2), "mm/h")
#' Qsim <- convolution(Rn = Rn, uh = uh)
#' @import sf stars
#' @importFrom units set_units drop_units
#' @export
convolution = function(Rn, ...) UseMethod("convolution")

#' @name convolution
#' @export
convolution.default <- function(Rn, uh, continuous=FALSE, ...){
  Q <- rep(0,length(Rn)+length(uh))
  Q <- .Call("c_convolution", rn=Rn, uh=uh)
  if(continuous){
    Q[1:(length(uh)-1)] <- NA
    return(Q[1:length(Rn)])
  }else{return(Q)}
}

#' @name convolution
#' @export
convolution.units <- function(Rn, uh, ...){
  Q <- convolution.default(Rn=units::drop_units(Rn), uh=units::drop_units(uh), ...)
  units(Q) <- units(Rn)
  return(Q)
}


#' @name convolution
#' @export
convolution.transfR = function(Rn, Rcol="RnSim", Qcol="Qsim", save_donor=FALSE, verbose=TRUE, ...){
  object <- Rn
  if(!"st"%in%names(object)) stop("No spatio-temporal array (st) of class 'stars' was found in the transfR object. See as_transfr().")
  if(!Rcol%in%names(object$st)) stop(paste0("Attribute '",Rcol,"' not found in the stars object.\nUse mixr() first to simulate net rainfall from neighbouring catchments, or select an existing attribute with the 'Rcol' argument."))
  Qsim <- matrix(NA,nrow=dim(object$st)[1],ncol=dim(object$st)[2])
  dim(Qsim) <- dim(object$st)
  for(i in 1:dim(object$st)[2]){
    if(verbose) progress("Computing simulated discharge for catchment ",i,dim(object$st)[2])
    iQsim <- convolution(Rn=object$st[[Rcol]][,i],uh=object$uh[[i]],continuous=TRUE)
    iQsim <- iQsim*st_area(st_geometry(object$st)[i])
    iQsim <- units::set_units(iQsim,"m^3/s")
    Qsim[,i] <- iQsim
  }
  object$st[[Qcol]] <- Qsim
  if(save_donor){
    iRcol <- grep("RnDonor",names(object$st))
    if(length(iRcol)==0) stop("save_donor is TRUE, but no 'RnDonor' attributes were found. Also use save_donor=TRUE with mixr() to retain each donor's Rn.")
    if(verbose) cat("\n")
    for(i in iRcol){
      if(verbose) progress("Additional discharge simulation using the single donor ",which(i==iRcol),length(iRcol),"")
      object <- convolution(Rn = object, Rcol = names(object$st)[i], Qcol = gsub("RnDonor","QsimDonor",names(object$st)[i]),
                            continuous = TRUE, save_donor = FALSE, verbose = FALSE, ...)
    }
    #Re-order columns just for readability of numerous attributes
    n <- names(object$st)
    icol <- c(grep("Qsim",n),grep("RnSim",n),grep("RnDonor",n),grep("Idonor",n),grep("Wdonor",n))
    object$st <- object$st[c(icol,(1:length(n))[!(1:length(n))%in%icol])]
  }
  units(object$st[[Qcol]]) <- units(iQsim) # could not find a way to keep the units
  return(object)
}