In TxWRI/ldc: Calculate and Plot Pollutant Load Duration Curves

knitr::opts_chunk$set(
    fig.path = "man/figures/README-",
    message = FALSE,
    warning = FALSE,
    collapse = TRUE,
    comment = "#>",
    dev = "ragg_png",
    dpi = 300,
    out.width = "100%"
)

ldc


ldc provides automated and fairly opinionated functions for generating pollutant load duration curves (LDCs) in freshwater streams. Due to the automated nature, there isn't much ability to adjust methodology or customize the generated LDCs since much of the calculation is abstracted away from the user.

ldc has three major functions:

• calc_ldc takes and input dataset of matched flow and pollutant concentrations to generate a table with exceedance probabilities grouped by user specified break points.

• summ_ldc uses the output from calc_ldc to generate a summary dataframe

• draw_ldc uses the output from both functions to generate a LDC figure as a ggplot object.

Installation

ldc is currently on Github. First install the remotes package then install ldc from Github:

remotes::install_github("TxWRI/ldc")

Example

An example using the data in the ldc package is shown below.

Setup and format data

library(ldc)
library(dplyr)
library(units)
library(ggplot2)

## this will calculate a ldc for indicator bacteria

## ldc uses the unit package to facilitate unit conversions
## we need to make the cfu unit first, since it isn't included 
## in the units package
install_unit("cfu")

## format the data for use in ldc
tres_palacios <- as_tibble(tres_palacios) |>
  ## flow must have units, here is is in cfs
  mutate(Flow = set_units(Flow, "ft^3/s"))|>
  ## pollutant concentration must have units
  mutate(Indicator_Bacteria = set_units(Indicator_Bacteria, "cfu/100mL"))

tres_palacios

Calculate exceedance probability

## specify the allowable concentration
allowable_concentration <- 126
## set the units
units(allowable_concentration) <- "cfu/100mL"

## calculate the exceedance probabilities along with
## allowable pollutant loads and measured pollutant loads
## at given probabilities
df_ldc <- calc_ldc(tres_palacios, 
                   Q = Flow, 
                   C = Indicator_Bacteria, 
                   allowable_concentration = allowable_concentration)

df_ldc

Summarize data

df_sum <- summ_ldc(df_ldc, 
                   Q = Flow, 
                   C = Indicator_Bacteria, 
                   Exceedance = P_Exceedance,
                   groups = Flow_Category,
                   method = "geomean")
df_sum

Plot LDC

draw_ldc(df_ldc, 
         df_sum, 
         y_lab = expression(paste(italic("E. coli"))),
         label_nudge_y = log10(1000)) + 
  scale_y_log10() +
  annotation_logticks(sides = "l") +
  theme_bw() +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        legend.direction = "vertical",
        panel.grid = element_blank())

Units

ldc relies on the units package to facilitate unit conversions and tracking of units across variables. This is handy if we want to transform units on the fly. In the above summary table, median daily flow volume is reported in units of 100mL/day. This isn't a logical unit to communicate, lets change it to million. gallons/day.

df_sum |>
  mutate(Median_Daily_Flow_Volume = set_units(Median_Daily_Flow_Volume, "1E6gallons/day")) -> df_sum
df_sum

cfu/day is a really big number. We can convert that to billion cfu/day.

df_sum |>
  mutate(Median_Flow_Load = set_units(Median_Flow_Load, "1E9cfu/day")) -> df_sum
df_sum

If we want to plot these, we also need to convert the df_ldc variables to matching units.

df_ldc |>
  mutate(Daily_Load = set_units(Daily_Load, "1E9cfu/day"),
         Allowable_Daily_Load = set_units(Allowable_Daily_Load, "1E9cfu/day")) -> df_ldc

Updated units will carry over to the plot:

draw_ldc(df_ldc, 
         df_sum, 
         y_lab = expression(paste(italic("E. coli"))),
         label_nudge_y = log10(1000)) + 
  scale_y_log10() +
  annotation_logticks(sides = "l") +
  theme_bw() +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        legend.direction = "vertical",
        panel.grid = element_blank())

TxWRI/ldc documentation built on Feb. 13, 2022, 9:22 a.m.