watertools: Function to Obtain and Assess Water Quality Data

## Functions to plot flow and load duration curves and calculate load reductions

## quiets concerns of R CMD check re: the .'s that appear in pipelines
if(getRversion() >= "2.15.1") utils::globalVariables(c("."))


# functions ---------------------------------------------------------------

# Geometric Mean ==========================================================

#' Geometric mean
#'
#' Generic function calculates the geometric mean
#'
#' @param x numeric vector of values
#' @param na.rm a logical value indicating weather \code{NA} values should
#' be stripped before the computation occurs. Defaults to \code{FALSE}
#' @param trim the fraction (o to 0.5) of observations to be trimmed from
#' each end of \code{x} before the geometric mean is computed. Values of
#' trim outside that range are taken as the nearest endpoint.
#' @param ... further arguments passed to or from other methods.
#' @export
#'
#' @examples x <- rlnorm(50, log(100), log(2))
#' xm <- gm_mean(x)
gm_mean = function(x, na.rm = TRUE, trim = 0, ...){
  exp(mean(log(x, ...), na.rm = na.rm, trim = trim, ...))
}


# FDC =====================================================================

# take a dataframe of daily mean streamflow values and generate the flow
# duration curve

# input should be a dataframe with
# column 0 = date in yyyy-mm-dd format
# column 2 = daily mean flow in cfs
#' Flow duration curve
#'
#' Calculates the flow exceedance probabilities from an input dataframe with daily mean streamflow values.
#'
#' @param data dataframe with first column dates in \code{yyyy-mm-dd} format, second column is assumed daily mean flow in cfs
#' @param filter a logical value indicating if zero values should be stripped after calculating the flow exceedance probabilities, useful for plotting the dataframe
#' @param ... further arguments passed to or from other methods
#' @import dplyr
#'
#' @return dataframe of flow measurements with column of FDC values
#' @export
fdc <- function(data, filter, ...){

  if(filter == TRUE) { #if filter = TRUE, calculate the cum dist, then remove the zeros. Handy for plotting
    d <- data %>%
      arrange(-.[[2]]) %>%
      mutate(cumdist = cume_dist(-.[[2]])) %>%
      filter(.[[2]] > 0)
  } else { #filter == FALSE, don't remove zero's, needed for calculating the flow exceedance probabilities
    d <- data %>%
      arrange(-.[[2]]) %>%
      mutate(cumdist = cume_dist(-.[[2]]))
  }

  d

}

# flowExceedance ==========================================================

# using the flow duration curve, identify the flows at desired flow
# exceedance values. Useful for LDC calculations and TMDLs.
# use a list of probs to return a df with flows at the midpoint of
# each exceedance category as desired by user
#' flowExceedance
#'
#' returns a dataframe with flows at specified flow exceedance probabilities
#'
#' @param data dataframe, assumed to be the output from \code{fdc()}.
#' @param probs value or vector of probabilities that flow should be computed.
#'
#' @examples \dontrun{#generate sample data
#' library(tidyr)
#' library(dplyr)
#'
#' df <- data.frame(
#'   date = seq(as.Date("2016-01-01"), as.Date("2016-12-31"), by = "day"),
#'   flow = rlnorm(366, 500, 75))
#'
#' dff <- fdc(df, filter = FALSE)
#' flowExceedance(dff, c(0.05,0.25,0.50,0.75,0.95))}
#' @export
flowExceedance <- function(data, probs){
  data = data$flow
  data.frame(exceedance = (100-probs*100), flow = stats::quantile(data, probs, type = 5))
}

# plotFDC =================================================================

## using fdc, plot the fdc in ggplot
#' plotFdc
#'
#' generates a flow duration curve as a \code{ggplot} object
#'
#' @param fdc expects output from \code{fdc()}
#' @param flowExceedance expects output from \code{flowExceedance()}
#' @param breaks optional argument specifying breaks used to to plot the flow duration curve
#' @import ggplot2
#'
#' @return \code{ggplot} object
#' @export
plotFdc <- function(fdc, flowExceedance, breaks){

  fdc$exceedance <- fdc$cumdist * 100

  p <- (ggplot(NULL) +
          geom_line(data = fdc, aes_(x = as.name(names(fdc)[4]),
                                     y = as.name(names(fdc)[2]))) +
          geom_point(data = flowExceedance,
                     aes_(x = as.name(names(flowExceedance)[1]),
                          y = as.name(names(flowExceedance)[2])),
                     shape = 7, color = "red") +
          scale_y_log10() +
          ylab("Flow (cfs)") + xlab("Percent of Days Flow Exceeded"))

  if(missing(breaks)) {
    p
  } else {
    p + scale_x_continuous(breaks = breaks)
  }
}


# ldc ---------------------------------------------------------------------

## calculate a load duration curve using the output from fdc, optionally join a dataframe of water quality measurements
##

#' ldc
#'
#' Calculates the load duration curve using the output from \code{fdc()}. Current built in conversions are for primary contact bacteria (126 MPN/100mL)
#'
#' @param fdc expects output from \code{fdc()}
#' @param conversion required character string. Valid inputs include: \code{"ecoli"}
#' @param observed optional dataframe of observed daily water quality measurement data, with first column in \code{yyyy-mm-dd} format.
#' @import dplyr
#' @export
ldc <- function(fdc, conversion, observed){

  if(conversion == "ecoli"){
    f <- (126/100) * 28316.8 * 86400 ###for bacteria cubic feet/sec * MPN/100mL * mL/cubic feet * sec/day = mpn/day
    g <- 100 #concentration per unit 100mL for bacteria
    l <- 28316.8 * 86400
  } else {
    f <- NA ## eventually add other conversions
    g <- 1
    l <- NA
  }

  if(missing(observed)) {
    fdc %>%
      mutate(loadcurve = .[[2]] * f)
  } else {
    fdc %>%
      mutate(loadcurve = .[[2]] * f) %>%
      left_join(observed) %>%
      mutate(measload = (.[[5]]/g) * .[[2]] * l)
  }

}

# plotLDC ----------------------------------------------------------

# plots the LDC using the output from ldc()

#' plotLdc
#'
#' @param ldc expects output from \code{ldc()}
#' @param breaks optional argumen specifying breaks used to plot the load duration curve
#' @import ggplot2
#' @return \code{ggplot} object
#' @export
plotLdc <- function(ldc, breaks){

  ldc$exceedance <- ldc$cumdist * 100

  p <- (ggplot(NULL) +
          geom_line(data = ldc, aes_(x = as.name(names(ldc)[7]),
                                     y = as.name(names(ldc)[4]))) +
          geom_point(data = ldc,
                     aes_(x = as.name(names(ldc)[7]),
                          y = as.name(names(ldc)[6])),
                     color = "red") +
          scale_y_log10() +
          ylab("Load Unit") + xlab("Percent of Days Load Exceeded"))

  if(missing(breaks)) {
    p
  } else {
    p + scale_x_continuous(breaks = breaks)
  }
}