In utah-dwq/wasteloadR: Tools for developing wasteloads

Knit to PDF

Note - to knit to PDF, you will need to install a version of latex in RStudio and phantomjs (only need to do this once). To do this run tinytex::install_tinytex() & `webshot::install_phantomjs()` in the R console. You can then knit to PDF clicking the dropdown menu beside the knit button in the above menu bar.

Packages

Install wasteloadR & wqTools if necessary.

devtools::install_github('utah-dwq/wqTools', upgrade = c('never')) 
devtools::install_github('utah-dwq/wasteloadR', upgrade = c('never'))
#Note- recommend updating other packages independently

library(dataRetrieval)
library(wqTools)
library(wasteloadR)
library(dplyr)
library(plotly)
library(knitr)
library(writexl)

User inputs

Site number 10171000 is Jordan River at 1700 S. You can find appropriate gauge sites via USGS NWIS or with wqTools::findSites(). Parameter code 00060 = discharge in cfs.

You can customize your query date range here, or query the full record and subset later after examining the daily values data.

Site <- "10171000" # Enter one USGS station number
startDate <- "1989-10-01" # Enter start date of gauge data download

Read USGS daily flow values and gauge info

gauge_data=dataRetrieval::readNWISdv(siteNumbers=Site,
  startDate=startDate, parameterCd='00060') %>% 
    dplyr::rename(discharge_cfs=X_00060_00003, discharge_qa_code=X_00060_00003_cd)

gauge_info <- readNWISsite(Site)

SiteName <- gauge_info$station_nm[1]

Export raw gauge data

write_xlsx(list(gauge_info=gauge_info, gauge_data=gauge_data), path="raw_gauge_data.xlsx")

Assign seasons and water year

Note - The seasons() function defines seasons as winter = c(1, 2, 3), spring = c(4, 5, 6), summer = c(7, 8, 9), fall = c(10, 11, 12) by default, but can be custom modified. See the irrigation/non-irrigation section for an example of custom season definitions.

gauge_data=gauge_data %>% 
  mutate(
    month=lubridate::month(Date), 
    season=seasons(Date), 
    wateryear=waterYear(Date)
  )

Plots of daily discharge

Based on a visual inspection, you may want to subset the data down to a specific time period, make custom season assignments, or take other QAQC actions.

plot_ly(data=gauge_data, x=~Date, y=~discharge_cfs, type='scatter', mode='lines') %>% 
  layout(
    xaxis = list(title = ''),
    yaxis = list(title = 'Discharge (cfs)')
    ) %>%   
    config(displaylogo = FALSE,
        modeBarButtonsToRemove = c(
            'sendDataToCloud',
            'hoverClosestCartesian',
            'hoverCompareCartesian',
            'lasso2d',
            'select2d'
        )
    )

plot_ly() %>% 
  add_trace(data=gauge_data, type='box', x=~wateryear, y=~discharge_cfs) %>% 
  layout(title = "",
        yaxis = list(title = 'Discharge (cfs)'),
        xaxis = list(side = 'left', title = 'Water Year', showline=TRUE)
    ) %>% 
    config(displaylogo = FALSE,
        modeBarButtonsToRemove = c(
            'sendDataToCloud',
            'hoverClosestCartesian',
            'hoverCompareCartesian',
            'lasso2d',
            'select2d'
        )
    )

plot_ly() %>% 
  add_trace(data=gauge_data, type='box', x=~season, y=~discharge_cfs) %>% 
  layout(title = "",
        yaxis = list(title = 'Discharge (cfs)'),
        xaxis = list(side = 'left', title = 'Season', showline=TRUE)
    ) %>% 
    config(displaylogo = FALSE,
        modeBarButtonsToRemove = c(
            'sendDataToCloud',
            'hoverClosestCartesian',
            'hoverCompareCartesian',
            'lasso2d',
            'select2d'
        )
    )

plot_ly() %>% 
  add_trace(data=gauge_data, type='box', x=~month, y=~discharge_cfs) %>% 
  layout(title = "",
        yaxis = list(title = 'Discharge (cfs)'),
        xaxis = list(side = 'left', title = 'Month', showline=TRUE)
    ) %>% 
    config(displaylogo = FALSE,
        modeBarButtonsToRemove = c(
            'sendDataToCloud',
            'hoverClosestCartesian',
            'hoverCompareCartesian',
            'lasso2d',
            'select2d'
        )
    )

Example: subset to selected years

SubStartDate <- "1989-10-01" # Enter start date of gauge data for use in low flow analysis
gauge_data=subset(gauge_data, Date>='1999-10-01')

nQy calculations

The nQy function is contained within the wasteloadR R package{target="_blank"}. It calculates two types of nQy values. One from a fitted probability of annual n-day minima, and one from the 1/n-th quantile of of annual n-day minima. See the help file with ?nQy() for more information.

Annual 7Q10

gauge_7Q10 <- nQy(data=gauge_data, n=7, y=10, date_col='Date', 
                     q_col='discharge_cfs', plot_fit=F, min_days=328)

gauge_7Q10_annualnday <- gauge_7Q10$annual_mins_fit
kable(gauge_7Q10_annualnday, caption=paste0(SiteName, " Annual nday Timeseries"), digits=3)

Examine annual 7Q10 fit plot

Note: by default, nQy() will display a fitted probability plot of n-day annual minima. The function also exports this plot as the object 'fit_plot'. You can use plot_fit=F when running the argument to turn off automatic plotting and examine plots manually. The fit_plot is a ggplot2 object and can be modified via normal ggplot methods. For example, here we add a custom title to the plot.

gauge_7Q10$fit_plot + 
  ggtitle(paste0(SiteName, ' Probability Fit Plot')) +
  theme(legend.position="bottom")

Seasonal 7Q10s

To calculate seasonal or other custom nQy values, subset your data to desired season and feed to nQy function. Use an appropriate min_days setting. In this case, we've divided the year into 4 seasons, and want our min_days to be approximately 90% of a season (365/4*0.9 = 82).

gauge_7Q10_winter=nQy(data=subset(gauge_data, season=='winter'), n=7, y=10,
  date_col='Date', q_col='discharge_cfs', plot_fit=F, min_days=82)
gauge_7Q10_spring=nQy(data=subset(gauge_data, season=='spring'), n=7, y=10,
  date_col='Date', q_col='discharge_cfs', plot_fit=F, min_days=82)
gauge_7Q10_summer=nQy(data=subset(gauge_data, season=='summer'), n=7, y=10,
  date_col='Date', q_col='discharge_cfs', plot_fit=F, min_days=82)
gauge_7Q10_fall=nQy(data=subset(gauge_data, season=='fall'), n=7, y=10,
  date_col='Date', q_col='discharge_cfs', plot_fit=F, min_days=82)

Other examples: custom seasons and other nQys

Irrigation/non-irrigation

Use the same approach as for seasonal values above. Modify the seasons and min_days as appropriate.

### User Input Required ###
irrMonth <- c(5,6,7,8,9) # Enter months of irrigation season
nonirrMonth <- c(1,2,3,4,10,11,12) # Enter months of non-irrigation season
# Minimum number of days in season with data to include year in analysis
irrminDays <- 120
# Minimum number of days in season with data to include year in analysis
nonirrminDays <- 160 

gauge_data=gauge_data %>%
  mutate(irrigation=seasons(Date, seasons=list('y'=irrMonth, 'n'=nonirrMonth)))
guage_irr=nQy(data=subset(gauge_data, irrigation=='y'), n=7, y=10, date_col='Date',
  q_col='discharge_cfs', min_days=irrminDays)
gauge_7Q10_nonirr=nQy(data=subset(gauge_data, irrigation=='n'), n=7, y=10,
  date_col='Date', q_col='discharge_cfs', min_days=nonirrminDays)

Other nQy types (e.g. 1Q10)

gauge_1Q10=nQy(data=gauge_data, n=1, y=10, date_col='Date', q_col='discharge_cfs')

Grab sample alternative approach

This example shows an alternative option to calculate grab sample flow measurement summary statistics in cases where an appropriate local USGS gauge is not available to calculate low flow statistics. To do this we'll query flow measurements at DWQ site locations from the Water Quality Portal, then group & summarize the data to calculate summary statistics for potential use as critical flow values.

Read grab sample data

Use Utah DWQ sites on the Logan River as an example.

grab_flow=readWQP(type="result",
  siteid=c("UTAHDWQ_WQX-4905250","UTAHDWQ_WQX-4905290", "UTAHDWQ_WQX-4905310"),
  CharacteristicName="Flow")
grab_flow$ActivityStartDate=as.Date(grab_flow$ActivityStartDate)

grab_flow_sites=readWQP(type="sites",
  siteid=c("UTAHDWQ_WQX-4905250","UTAHDWQ_WQX-4905290", "UTAHDWQ_WQX-4905310"))

grab_flow=merge(grab_flow, grab_flow_sites)

# Check on units
kable(table(grab_flow$CharacteristicName, grab_flow$ResultMeasure.MeasureUnitCode))

Summarize grab sample flow

In this example, we group by site code, parameter, and season to calculate summary statistics for each grouping. You can add or drop grouping columns to calculate stats on different groups (e.g. remove the site code column group to pool data from multiple sites, or remove season to get annual stats).

grab_flow_summary=grab_flow %>% 
  mutate(
    season=seasons(ActivityStartDate), 
  ) %>%
  group_by(MonitoringLocationIdentifier, MonitoringLocationName, 
    CharacteristicName, season) %>%
  summarize(
    count=n(),
    min=min(ResultMeasureValue),
    pctile10=quantile(ResultMeasureValue,0.1), 
    .groups="drop"
  )


kable(grab_flow_summary)

utah-dwq/wasteloadR documentation built on Oct. 19, 2024, 11:04 p.m.