R/utility_functions.R

In ChrisBogner/ViscousFlow: Analysis of tracer experiments with the viscous flow approach by Peter Germann

#' Plot the tracer breakthrough curve
#'
#' @param tracer_data tibble or data.frame. Tracer breakthrough data. The first column must be time, the second the tracer concentration.
#' @param time_threshold numeric. Time of the breakthrough (i.e. start of the concentration increase)
#' @param xlab_tag string or expression for the x label. Default is Time (sec).
#' @param ylab_tag string or expression for the y label. Default is c/c0.
#' @param file_name string. The name for the pdf file to save the figure. Default is NULL.
#' @param save_pdf logical. Should de figure be saved to a pdf file. Default is FALSE.
#'
#' @return ggplot object
#' @export
#'
#' @examples
#' data(tracer)
#' plot_tracer(tracer_data = tracer, time_threshold = 40000)
plot_tracer <- function(tracer_data, time_threshold = NULL,
                        xlab_tag = 'Time (sec)', ylab_tag = expression(c/c[0]),
                        file_name = NULL,
                        save_pdf = FALSE) {

  g <- ggplot2::ggplot(data = tracer_data,
                       ggplot2::aes_string(x = colnames(tracer_data)[1],
                             y = colnames(tracer_data)[2])) +
    ggplot2::geom_line() +
    ggplot2::geom_point() +
    ggplot2::xlab(xlab_tag) +
    ggplot2::ylab(ylab_tag)
  if(!is.null(time_threshold)) {
    g <- g + ggplot2::geom_vline(xintercept = time_threshold, colour = 'blue', lty = 2)
  }
  else g
  if(save_pdf) ggplot2::ggsave(file_name)
  g
}

###########################

#' Extracts data from scales
#'
#' Extracts weight data from ascii files generated by HTerm 0.8.1beta
#'
#' @param data_orig the original data, the format is "hh:mm:ss.sss:\n weight unit
#' @param weight_unit character, weight unit transmitted by the scale: g or kg
#' @param do_plot logical. Should the extracted data be plotted. Default is TRUE.
#'
#' @return a tibble. First column time in sec, second column weight in unit
#' @export
#'
#' @examples
extract_data <- function(data_orig, weight_unit, do_plot = T) {
  time_xs <- data_orig %>%
    dplyr::slice(seq(1, dim(data_orig)[1], 2)) %>%
    dplyr::pull(var = 1) %>%
    stringr::str_sub(start = 1, end = 8) %>%
    lubridate::hms()

  time_rel <- lubridate::time_length(time_xs - time_xs[1], unit = 'seconds')

  dat <- data_orig %>%
    dplyr::slice(seq(2, dim(data_orig)[1], 2)) %>%
    dplyr::pull(var = 1) %>%
    stringr::str_remove_all('\\s') %>%
    stringr::str_sub(start = 1, end = -(stringr::str_length(weight_unit)+1)) %>%
    as.numeric()
  res <- tibble::tibble('time_rel_sec' = time_rel, dat)
  colnames(res)[2] <- paste0('weight_', weight_unit)

  g <- ggplot2::ggplot(res, ggplot2::aes_string(x = 'time_rel_sec',
                                                y = colnames(res)[2])) +
    ggplot2::geom_line()

  if(do_plot) print(g)
  res
}

##################################

#' Calculates the flow rate from cumulative drainage
#'
#' The relative weights are calculated by taking the first derivative using the function D1 from library sfsmisc. Then, they are transformed into a flow rate. Both, a smoothed version of the derivative and the trivial differences are returned.
#'
#' @param scale_s tibble. First column time in sec, second column cumulative weight of the drained water
#' @param rho_H20 numeric. Density of water, default is 0.99704 g/mL or g/cm^3
#' @param D numeric. Diameter of the soil column in m.
#' @param delta_out numeric. Output sampling time interval, default is 30 sec. Must equal at least the length of the original sampling interval.
#' @param convert_mm_h logical. Should the flow rate be converted to mm/h. Default is TRUE.
#'
#' @return tibble. First column is time in sec, second is the flow rate in mm/h
#' @export
#'
#' @examples
calculate_drainage_flow <- function(scale_s, rho_H20 = 0.99704, D, delta_out = 30, convert_mm_h = TRUE) {

  delta_t <- scale_s %>%
    dplyr::pull(var = 1) %>%
    diff() %>%
    mean() %>%
    round()

   sampling_time <- scale_s %>%
     dplyr::slice(seq(1, dim(scale_s)[1], by = delta_out/delta_t)) %>%
     dplyr::pull(var = 1)

  res <- scale_s %>%
    dplyr::pull(var = 2) %>%
    diff()

  res_clean <- ifelse(abs(res) >= stats::quantile(res, probs = 0.99), NA, res)
  res_clean <- ifelse(res_clean < 0, NA, res_clean)
  res_clean <- zoo::na.approx(res_clean, na.rm = FALSE)
  res_clean <- ifelse(is.na(res_clean), 0, res_clean)
  res_clean <- c(0,cumsum(res_clean))

  # resample at sampling_time and join
  dat_clean <- tibble::tibble(x = unlist(scale_s[,1]), y = res_clean)
  sampling_join <- tibble::tibble(x = sampling_time)
  dat_deriv <- dplyr::inner_join(sampling_join, dat_clean, by = 'x')


  # smoothed derivative
  flow_rate_smoothed <- sfsmisc::D1ss(x = dat_deriv$x,
                                      y = dat_deriv$y,
                                      xout = dat_deriv$x,
                                      spar.offset = 0)

  # trivial differences
  flow_rate <- sfsmisc::D1tr(x = dat_deriv$x,
                             y = dat_deriv$y)

  flow_rate <- ifelse(flow_rate < 0, 0, flow_rate)
  flow_rate_smoothed <- ifelse(flow_rate_smoothed < 0, 0, flow_rate_smoothed)

  if(convert_mm_h) {
    flow_rate <- flow_rate * 3600 /(rho_H20 * 1000 * (D/2)^2 * pi)
    flow_rate_smoothed <- flow_rate_smoothed * 3600 /(rho_H20 * 1000 * (D/2)^2 * pi)
  }

  tibble::tibble(sec = sampling_time, flow_rate, flow_rate_smoothed)
}