R/inversion.R

In transfR: Transfer of Hydrograph from Gauged to Ungauged Catchments

#' Estimate net rainfall by inversion
#'
#' Estimate net rainfall by inverse modelling, where the model is a convolution between net rainfall
#' and a unit hydrograph in order to simulate discharge.
#' @name inversion
#' @param Qobs discharge vector or object of class \code{transfR}. If no unit is provided,
#' \code{Qobs} is assumed to be in [mm/h]
#' @param uh unit hydrograph vector
#' @param deltat time step of the time series. If no unit is provided, \code{deltat} is assumed to be in [min]
#' @param RnAp a priori estimate of net rainfall. If no unit is provided, \code{RnAp} is assumed to be in [mm/h]
#' @param Bd parameter used to maintain a minimum value of standard deviation for low discharge values.
#' If no unit is provided, \code{Bd} is assumed to be in [mm/h]
#' @param Dd decorrelation time of discharge errors. If no unit is provided, \code{Dd} is assumed to be
#' in [h]
#' @param Bp parameter used to maintain a minimum value of standard deviation for low net rainfall values.
#' If no unit is provided, \code{Bp} is assumed to be in [mm/h]
#' @param Tp decorrelation time of net rainfall errors. If no unit is provided, \code{Tp} is assumed to
#' be in [h]
#' @param Ad parameter equivalent to the coefficient of variation of the discharge measurement error. If
#' no unit is provided, \code{Ad} is assumed to be dimensionless
#' @param Ap parameter equivalent to the coefficient of variation of the net rainfall error. If no unit
#' is provided, \code{Ap} is assumed to be dimensionless
#' @param warmup length of the warmup period. If no unit is provided, \code{warmup} is assumed to be in [days]
#' @param cooldown length of the period removed at the end of the simulation. If no unit is provided,
#' \code{cooldown} is assumed to be in [days]
#' @param dosplit logical, if true the inversion is performed by
#' subperiods of length defined by \code{split}
#' @param split length of the subperiods if dosplit is true. If no unit is provided, \code{split} is assumed to be
#' in [days]
#' @param fixedpar logical; if FALSE, Ap and Ad are calibrated dynamically according to the coefficient of variation of
#' RnAp and Qobs respectively (see details)
#' @param parallel logical indicating whether the computation should be parallelised.
#' Subperiods are parallelised for a single time series, whereas catchments are
#' parallelised for a \code{transfR} object
#' @param cores the number of cores to use for parallel execution if \code{parallel} is TRUE. If not specified, the number of cores is set to the value of \code{parallel::detectCores()}
#' @param cluster an optional parallel cluster to reuse
#' @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{Qobs}. If \code{Qobs} is a transfR object,
#' it is returned with the computed results added as new attributes.
#' @import sf stars doParallel foreach
#' @importFrom units set_units as_units drop_units
#' @importFrom stats na.omit sd
#' @details Given a convolution between the unit hydrograph (\code{uh}) and net rainfall that simulates
#' streamflow at the outlet (\code{Qobs}), where net rainfall is the only unknown variable, this function estimates
#' net rainfall by inversion \insertCite{Tarantola1982,Menke1989,Boudhraa2018}{transfR}. It requires an
#' a priori estimate of net rainfall (which can be obtained with \link{rapriori}) and a description
#' of errors in discharge (\code{Ad}, \code{Bd}, \code{Dd}) and net rainfall (\code{Ap},
#' \code{Bp}, \code{Tp}); these errors are assumed to be Gaussian and unbiased. Default values of these parameters
#' are taken from \insertCite{deLavenne2016;textual}{transfR}. If \code{fixedpar} is deactivated, \code{Ap}
#' is estimated at 20% of the coefficient of variation of RnAp, and \code{Ad} is estimated at 5% of the coefficient
#' of variation of Qobs.
#'
#' It is recommended to use \code{warmup} and \code{cooldown} periods to reduce oscillations
#' caused by inversion.
#'
#' If \code{Qobs} is a \code{transfR} object, results are stored in a new space-time attribute
#' called "RnInv".
#' @seealso
#' \link{rapriori}
#' @references
#' \insertRef{Boudhraa2018}{transfR}
#'
#' \insertRef{deLavenne2016}{transfR}
#'
#' \insertRef{Menke1989}{transfR}
#'
#' \insertRef{Tarantola1982}{transfR}
#' @examples
#' \donttest{data(Oudon)
#' icatch <- 1 # Catchment index
#' itime <- 1:1000 # Using the first values for a quicker example
#' Qobs <- Oudon$obs[["Qobs"]][itime,icatch]
#' Qspec <- units::set_units(Qobs/st_area(st_geometry(Oudon$obs)[icatch]), "mm/h")
#' deltat <- units::set_units(1, "h")
#' uc <- velocity(hl = Oudon$hl[[icatch]])
#' uh <- uh(hl = Oudon$hl[[icatch]], uc = uc, deltat = units::set_units(1,"h"))$prob
#' RnAp <- rapriori(Qobs = Qspec, lagtime = lagtime(hl = Oudon$hl[[icatch]], uc = uc),
#' deltat = deltat)
#' RnInv <- inversion(Qobs = Qspec, RnAp = RnAp, uh = uh, deltat = deltat, parallel = TRUE, cores=2)}
#' @export
inversion <- function(Qobs,...) UseMethod("inversion")

#' @name inversion
#' @export
inversion.default <- function(Qobs, uh, RnAp, deltat, ...){
  #--- Assumed units
  Qobs <- units::set_units(Qobs,"mm/h")
  uh <- units::set_units(uh,1)
  RnAp <- units::set_units(RnAp,"mm/h")
  deltat <- units::set_units(deltat,"min")
  #--- Inversion
  inversion.units(Qobs = Qobs, uh = uh, RnAp = RnAp, deltat = deltat, ...)
}

#' @name inversion
#' @export
inversion.units <- function(Qobs, uh, RnAp, deltat, Bd = 0.01, Dd = 1, Bp = 0.001, Tp = 20, Ad = 0.01,
                            Ap = 0.9, warmup = 10, cooldown = 8, dosplit = TRUE, split = 30,
                            fixedpar = TRUE, parallel = FALSE, cores = NULL, cluster = NULL, ...){

  if(parallel && is.null(cluster)){
    if(missing(cores) || is.null(cores)) cores <- parallel::detectCores()
    cluster <- parallel::makeCluster(cores)
    on.exit(parallel::stopCluster(cluster))
  }

  out_units <- units(RnAp)

  #--- Define some parameter values according to inputs characteristics
  if(!fixedpar){
    if(missing(Ad)) Ad <- units::drop_units(stats::sd(Qobs,na.rm = T)/mean(Qobs,na.rm = T)*0.05)
    if(missing(Ap)) Ap <- units::drop_units(stats::sd(RnAp,na.rm = T)/mean(RnAp,na.rm = T)*0.2)
  }

  #--- Inversion parameters
  Ad <- units::set_units(Ad,1)
  Bd <- units::set_units(Bd,"mm/h")
  Ap <- units::set_units(Ap,1)
  Bp <- units::set_units(Bp,"mm/h")
  Dd <- units::set_units(Dd,"h")
  Tp <- units::set_units(Tp,"h")

  #--- Warmup and cooldown periods
  deltat   <- units::set_units(deltat,"min")
  warmup   <- units::set_units(warmup,"days")
  cooldown <- units::set_units(cooldown,"days")
  split    <- units::set_units(split,"days")
  npdt_warmup   <- ceiling(units::drop_units(warmup/deltat))
  npdt_cooldown <- ceiling(units::drop_units(cooldown/deltat))
  npdt_split    <- ceiling(units::drop_units(split/deltat))

  #--- Check units
  if(length(units(Qobs/units::set_units(1,"m"))$numerator)!=0) stop("Qobs should be a specific discharge. Check units.")
  if(length(units(RnAp/units::set_units(1,"m"))$numerator)!=0) stop("RnAp should be a specific discharge. Check units.")
  if(!(length(units(Ad)$numerator)==0 & length(units(Ad)$denominator)==0)) stop("Ad should be without units.")
  if(!(length(units(Bd/units::set_units(1,"mm/h"))$numerator)==0 & length(units(Bd/units::set_units(1,"mm/h"))$denominator)==0)) stop("Bd should be a specific discharge. Check units.")
  if(!(length(units(Ap)$numerator)==0 & length(units(Ap)$denominator)==0)) stop("Ap should be without units.")
  if(!(length(units(Bp/units::set_units(1,"mm/h"))$numerator)==0 & length(units(Bp/units::set_units(1,"mm/h"))$denominator)==0)) stop("Bp should be a specific discharge. Check units.")
  if(!(length(units(Dd/units::set_units(1,"h"))$numerator)==0 & length(units(Dd/units::set_units(1,"h"))$denominator)==0)) stop("Dd should be a time. Check units.")
  if(!(length(units(Tp/units::set_units(1,"h"))$numerator)==0 & length(units(Tp/units::set_units(1,"h"))$denominator)==0)) stop("Tp should be a time. Check units.")

  Rn <- inversion_numeric(
    Qobs = units::drop_units(units::set_units(Qobs, "mm/h")),
    uh = units::drop_units(units::set_units(uh, 1)),
    RnAp = units::drop_units(units::set_units(RnAp, "mm/h")),
    deltat = units::drop_units(units::set_units(deltat, "min")),
    Bd = units::drop_units(units::set_units(Bd, "mm/h")),
    Dd = units::drop_units(units::set_units(Dd, "h")),
    Bp = units::drop_units(units::set_units(Bp, "mm/h")),
    Tp = units::drop_units(units::set_units(Tp, "h")),
    Ad = units::drop_units(units::set_units(Ad, 1)),
    Ap = units::drop_units(units::set_units(Ap, 1)),
    warmup = units::drop_units(units::set_units(warmup, "days")),
    cooldown = units::drop_units(units::set_units(cooldown, "days")),
    dosplit = dosplit,
    split = units::drop_units(units::set_units(split, "days")),
    fixedpar = TRUE,
    parallel = parallel,
    cores = cores,
    cluster = cluster,
    ...
  )
  units(Rn) <- units::as_units("mm/h")
  units::set_units(Rn, out_units, mode = "standard")
}

inversion_numeric <- function(Qobs, uh, RnAp, deltat, Bd = 0.01, Dd = 1, Bp = 0.001, Tp = 20, Ad = 0.01,
                              Ap = 0.9, warmup = 10, cooldown = 8, dosplit = TRUE, split = 30,
                              fixedpar = TRUE, parallel = FALSE, cores = NULL, cluster = NULL, ...){

  #--- Define some parameter values according to inputs characteristics
  if(!fixedpar){
    if(missing(Ad)) Ad <- stats::sd(Qobs,na.rm = T)/mean(Qobs,na.rm = T)*0.05
    if(missing(Ap)) Ap <- stats::sd(RnAp,na.rm = T)/mean(RnAp,na.rm = T)*0.2
  }

  #--- Warmup and cooldown periods
  npdt_warmup   <- ceiling(warmup*24*60/deltat)
  npdt_cooldown <- ceiling(cooldown*24*60/deltat)
  npdt_split    <- ceiling(split*24*60/deltat)

  #--- Managing periods with NA
  bna <- !is.na(RnAp)&!is.na(Qobs)
  nna <- which(bna)
  if(length(nna)<=(npdt_warmup+npdt_cooldown)) return(rep(NA,length(Qobs)))

  #--- Detect if several periods without NA and run individual inversions for each period
  period <- as.numeric(bna[1])
  for(i in 1:(length(bna)-1)) period[i+1] <- period[i]+abs(bna[i]-bna[i+1])
  if(max(period)>2){
    Rn <- NA
    for(p in unique(period[bna])) Rn[period==p] <- inversion_numeric(Qobs = Qobs[period==p], uh = uh, RnAp = RnAp[period==p], deltat = deltat, Bd = Bd, Dd = Dd, Bp = Bp, Tp = Tp, Ad = Ad, Ap = Ap, warmup = warmup, cooldown = cooldown, dosplit = dosplit, split = split, fixedpar = TRUE, parallel = parallel, cores = cores, cluster = cluster)
    return(Rn)
  }

  #--- Checking
  if(any((nna[2:length(nna)]-nna[2:length(nna)-1])>1)) stop("Time series must be continuous. NA values can be handled only at the beginning or the end of a time series.")
  # if(sum(!bna)>0) warning(paste0(sum(!bna)," NA values at the beginning or end of Qobs or RnAp times series. This part of the time series will not be used."))
  if(warmup<10) warning("Warmup period is short (< 10 days); you might observe oscillations in the Rn time series.")
  if(cooldown<8) warning("Cooldown period is short (< 8 days); you might observe oscillations in the Rn time series.")
  if(npdt_warmup<(5*length(uh))) warning("Warmup period might be too short; you might observe oscillations in the Rn time series.")
  if(npdt_cooldown<(4*length(uh))) warning("Cooldown period might be too short; you might observe oscillations in the Rn time series.")

  #--- Vector dimension
  Q <- Qobs[nna]
  R <- RnAp[nna]
  npdtR<-length(R)
  npdtQ<-length(Q)
  qstr <- format(Q, scientific = FALSE, trim = TRUE, drop0trailing = TRUE, decimal.mark = ".")
  ndec <- max(nchar(sub("^[^.]*\\.?", "", qstr)))

  #--- Default values
  ZRn <- 0
  NRn <- NA

  if(dosplit&(length(Q)>npdt_split)){
    if((npdt_warmup+npdt_cooldown)>=npdt_split) stop("The entire subperiod is used for initialisation: (warmup+cooldown)>=split. Increase the value of 'split' or decrease the values of 'warmup' and 'cooldown'.")
    cuts_start <- seq(from=1,to=(length(Q)+npdt_split-npdt_warmup-npdt_cooldown),by=(npdt_split-npdt_warmup-npdt_cooldown))
    cuts_end   <- cuts_start+npdt_split-1
    cuts_end[cuts_end>length(Q)] <- length(Q)
    cuts_start <- cuts_start[!duplicated(cuts_end)]
    cuts_end   <- cuts_end[!duplicated(cuts_end)]
    if(parallel){
      if(is.null(cluster)) stop("A cluster must be supplied for parallel inversion.")
      doParallel::registerDoParallel(cl=cluster)
      tmp <- foreach::"%dopar%"(foreach::foreach(i = 1:length(cuts_start), .combine='c',
                                                 .export = c("inversion_numeric", "uh2mat")),
                                na.omit(inversion_numeric(Qobs = Q[cuts_start[i]:cuts_end[i]], uh = uh, RnAp = R[cuts_start[i]:cuts_end[i]], deltat = deltat, Bd = Bd, Dd = Dd, Bp = Bp, Tp = Tp, Ad = Ad, Ap = Ap, warmup = warmup, cooldown = cooldown, dosplit=FALSE, split = split, fixedpar = TRUE, parallel = FALSE)))
      Rn <- c(rep(NRn,npdt_warmup),tmp,rep(NRn,npdt_cooldown))
    }else{
      Rn <- rep(NRn,length(Q))
      for(i in 1:length(cuts_start)){
        tmp <- inversion_numeric(Qobs = Q[cuts_start[i]:cuts_end[i]], uh = uh, RnAp = R[cuts_start[i]:cuts_end[i]], deltat = deltat, Bd = Bd, Dd = Dd, Bp = Bp, Tp = Tp, Ad = Ad, Ap = Ap, warmup = warmup, cooldown = cooldown, dosplit=FALSE, split = split, fixedpar = TRUE, parallel = FALSE)
        Rn[cuts_start[i]:cuts_end[i]][!is.na(tmp)] <- tmp[!is.na(tmp)]
      }
    }
    # If NA in Qobs or RnAp
    FullRn <- rep(NRn,length(Qobs))
    FullRn[nna] <- Rn
    return(FullRn)
  }

  #--- Transfer function matrix
  M <- uh2mat(uh,nrow=npdtQ,ncol=npdtR)
  tM <- Matrix::t(M)

  #------------------------------------
  # Covariance matrix on runoff
  #------------------------------------
  #--- Standard deviation (confidence in measurement)
  TQobs = Ad * Q + Bd
  #--- Matrix computation
  u<-npdtQ-1
  v<-0:u
  Mtmp1<-matrix(rep(0:u,npdtQ),npdtQ,npdtQ)
  Mtmp2<-matrix(rep(0:-u,npdtQ),npdtQ,npdtQ,byrow=TRUE)
  MTQobs3<-Mtmp1+Mtmp2

  CovQ_tmp <- -0.5*((abs(MTQobs3)*deltat/60)/Dd)^2
  CovQ<-outer(TQobs,TQobs)*exp(CovQ_tmp)
  #------------------------------------
  # Covariance matrix on net rainfall
  #------------------------------------
  #--- Standard deviation (confidence in a priori rainfall)
  TRnAp = Ap * R + Bp
  #--- Matrix computation
  u<-npdtR-1
  v<-0:u
  Mtmp1<-matrix(rep(0:u,npdtR),npdtR,npdtR)
  Mtmp2<-matrix(rep(0:-u,npdtR),npdtR,npdtR,byrow=TRUE)
  MTRnAp3<-Mtmp1+Mtmp2

  CovAp_tmp <- -0.5*((abs(MTRnAp3)*deltat/60)/Tp)^2
  CovAp<-outer(TRnAp,TRnAp)*exp(CovAp_tmp)

  #------------------------------------
  # Solution: net rainfall estimation
  #------------------------------------
  # Same equation as before; M is stored sparsely to avoid multiplying zeros:
  # Rn <- as.vector(R + CovAp%*%t(M) %*% solve((M %*% CovAp %*% t(M) + CovQ)) %*% (Q - M%*%R))
  A <- as.matrix(M %*% CovAp %*% tM) + CovQ
  res <- Q - as.vector(M %*% R)
  Rn <- as.vector(R + CovAp %*% as.vector(tM %*% solve(A, res)))
  #--- Removing negative Rn, warmup period and cooldown period
  Rn[Rn<ZRn] <- ZRn
  if(npdt_warmup>0)   Rn[1:npdt_warmup] <- NRn
  if(npdt_cooldown>0) Rn[(length(Rn)-npdt_cooldown+1):length(Rn)] <- NRn
  # If NA in Qobs or RnAp
  FullRn <- rep(NRn,length(Qobs))
  FullRn[nna] <- round(Rn,ndec)
  return(FullRn)
}

uh2mat <- function(uh,nrow,ncol){
  k <- length(uh)
  j <- rep(seq_len(ncol), each = k)
  i <- j + rep(seq_len(k) - 1, times = ncol)
  x <- rep(uh, times = ncol)
  keep <- i <= nrow
  Matrix::sparseMatrix(i = i[keep], j = j[keep], x = x[keep], dims = c(nrow,ncol))
}


#' @name inversion
#' @export
inversion.transfR <- function(Qobs, verbose=TRUE, parallel=FALSE, cores=NULL, cluster=NULL, ...){
  object <- Qobs
  if(!"deltat"%in%names(object)) stop("Time step attribute (deltat) is missing in transfR object. See as_transfr().")
  if(!"Qobs"%in%names(object$st)) stop("Discharge observations (Qobs) are missing from the spatio-temporal arrays (st) of the transfR object. See as_transfr().")
  if(!"RnAp"%in%names(object$st)) stop("The a priori estimate of net rainfall (RnAp) is missing from the spatio-temporal arrays (st) of the transfR object. See rapriori().")
  ncatchments <- dim(object$st)[2]
  if(parallel && is.null(cluster)){
    if(missing(cores) || is.null(cores)) cores <- parallel::detectCores()
    cores <- min(cores,ncatchments)
    cluster <- parallel::makeCluster(cores)
    on.exit(parallel::stopCluster(cluster))
  }
  if(object$deltat==units::set_units(1,"day")){out_units <- units::as_units("mm/d")}else{out_units <- units::as_units("mm/h")}
  object$st$RnInv <- NA
  if(parallel){
    if(verbose) cat("Computing inversion for",ncatchments,"catchments in parallel\n")
    jobs <- vector("list",ncatchments)
    for(i in seq_len(ncatchments)){
      Qspec <- object$st$Qobs[,i]/st_area(st_geometry(object$st)[i])
      jobs[[i]] <- list(Qobs = Qspec, uh = object$uh[[i]], RnAp = object$st$RnAp[,i],
                        deltat = object$deltat)
    }
    doParallel::registerDoParallel(cl=cluster)
    dots <- list(...)
    job <- NULL # Define the foreach iterator for R CMD check
    results <- foreach::"%dopar%"(
      foreach::foreach(job = jobs, .packages = "transfR"),
      do.call(inversion, c(job, dots, list(parallel = FALSE)))
    )
    for(i in seq_len(ncatchments)) object$st$RnInv[,i] <- results[[i]]
  }else{
    for(i in seq_len(ncatchments)){
      if(verbose) progress("Computing inversion for catchment ",i,ncatchments)
      Qspec <- object$st$Qobs[,i]/st_area(st_geometry(object$st)[i])
      object$st$RnInv[,i] <- inversion(Qobs = Qspec, uh = object$uh[[i]],
                                       RnAp = object$st$RnAp[,i], deltat = object$deltat,
                                       parallel = FALSE, ...)
    }
  }
  object$st[["RnInv"]] <- units::set_units(object$st[["RnInv"]],out_units,mode="standard") # could not find a way to keep the units provided by inversion
  return(object)
}