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_ldctakes and input dataset of matched flow and pollutant concentrations to generate a table with exceedance probabilities grouped by user specified break points. -
summ_ldcuses the output fromcalc_ldcto generate a summary dataframe -
draw_ldcuses the output from both functions to generate a LDC figure as a ggplot object.
ldc is currently on Github. First install the remotes package then install ldc from Github:
remotes::install_github("TxWRI/ldc")An example using the data in the ldc package is shown below.
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
#> # A tibble: 7,671 x 4
#> site_no Date Flow Indicator_Bacteria
#> <chr> <date> [ft^3/s] [cfu/100mL]
#> 1 08162600 2000-01-01 0.84 NA
#> 2 08162600 2000-01-02 3 NA
#> 3 08162600 2000-01-03 3.4 NA
#> 4 08162600 2000-01-04 2.6 NA
#> 5 08162600 2000-01-05 1.6 NA
#> 6 08162600 2000-01-06 3.2 NA
#> 7 08162600 2000-01-07 11 NA
#> 8 08162600 2000-01-08 17 NA
#> 9 08162600 2000-01-09 22 NA
#> 10 08162600 2000-01-10 18 NA
#> # ... with 7,661 more rowsCalculate 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
#> # A tibble: 7,671 x 9
#> site_no Date Flow Indicator_Bacteria Daily_Flow_Volume Daily_Load
#> <chr> <date> [ft^3/s] [cfu/100mL] [100mL/d] [cfu/d]
#> 1 08162600 2000-08-18 0.22 NA 5382466. NA
#> 2 08162600 2000-03-07 0.42 NA 10275617. NA
#> 3 08162600 2000-03-06 0.6 NA 14679453. NA
#> 4 08162600 2000-08-22 0.75 NA 18349317. NA
#> 5 08162600 2000-03-08 0.78 NA 19083289. NA
#> 6 08162600 2000-08-17 0.78 NA 19083289. NA
#> 7 08162600 2000-09-11 0.8 NA 19572604. NA
#> 8 08162600 2000-01-01 0.84 NA 20551235. NA
#> 9 08162600 2000-08-21 0.93 NA 22753153. NA
#> 10 08162600 2000-03-05 1 NA 24465755. NA
#> # ... with 7,661 more rows, and 3 more variables: Allowable_Daily_Load [cfu/d],
#> # P_Exceedance <dbl>, Flow_Category <fct>Summarize data
df_sum <- summ_ldc(df_ldc,
Q = Flow,
C = Indicator_Bacteria,
Exceedance = P_Exceedance,
groups = Flow_Category,
method = "geomean")
df_sum
#> # A tibble: 3 x 6
#> Flow_Category Median_Flow Median_P Geomean_C Median_Daily_Fl~ Median_Flow_Load
#> <fct> [ft^3/s] <dbl> [cfu/100mL] [100mL/d] [cfu/d]
#> 1 High Flows 58.9 0.200 259. 1441032996. 373699594298.
#> 2 Medium Flows 16.5 0.600 83.6 403684965. 33764600228.
#> 3 Low Flows 8.66 0.900 67.1 211995771. 14229395256.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())
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
#> # A tibble: 3 x 6
#> Flow_Category Median_Flow Median_P Geomean_C Median_Daily_Fl~ Median_Flow_Load
#> <fct> [ft^3/s] <dbl> [cfu/100mL] [1E6gallons/d] [cfu/d]
#> 1 High Flows 58.9 0.200 259. 38.1 373699594298.
#> 2 Medium Flows 16.5 0.600 83.6 10.7 33764600228.
#> 3 Low Flows 8.66 0.900 67.1 5.60 14229395256.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
#> # A tibble: 3 x 6
#> Flow_Category Median_Flow Median_P Geomean_C Median_Daily_Fl~ Median_Flow_Load
#> <fct> [ft^3/s] <dbl> [cfu/100mL] [1E6gallons/d] [1E9cfu/d]
#> 1 High Flows 58.9 0.200 259. 38.1 374.
#> 2 Medium Flows 16.5 0.600 83.6 10.7 33.8
#> 3 Low Flows 8.66 0.900 67.1 5.60 14.2If 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_ldcUpdated 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())