Flow and Salinity Data Sets

Purpose

The purpose of this vignette is to cover the basics for creating seasonally detrended flow and seasonally detrended salinity data sets. These data sets are used in support of applying gam4. In gam4, the dependent variable is modeled with a non-linear term as a function of year; a seasonality term includes an interaction term which allows seasonality to vary over the period of record; and includes a hydrologic term that allows for factoring wet/dry conditions in the model.

The seasonally detrended flow and salinity data sets are created using the detrended.flow and detrended.salinity functions included in baytrends. Refer to the help for each of these functions (e.g., ??detrended.flow or ??detrended.salinity) for the specific computational steps involved with each function.

Load baytrends

The first step is to load the baytrends library. Loading baytrends will result in loading multiple additional libraries as well.

library(baytrends)
#> Loading:baytrends v2.0.7
#> 
#> Attaching package: 'baytrends'
#> The following object is masked from 'package:stats':
#> 
#>     nobs

Seasonally Detrended Flow

The function, detrended.flow, returns a list of seasonally detrended flow and companion statistics. It relies on USGS’ dataRetrieval package to retrieve daily flow data.

It is the user’s responsibility to save the list that is returned from the detrended.flow function as flow.detrended for integration with baytrends.

For purposes of the Chesapeake Bay Program, it is expected that the user would typically identify all nine of the River Input Monitoring (RIM) stations. The nine RIM stations for the Chesapeake Bay have been loaded as a data frame in baytrends, usgsGages.

Table 1. Chesapeake Bay River Input Monitoring Stations

The below code chunk is an example of retrieving and seasonally detrending daily flow data for all nine of the RIM stations from 1983-2017. To retrieve data for other USGS gages, the variables, usgsGagesID and siteName would be assigned with vectors that contain the list of USGS gage IDs and their site names. The vector of gage IDs are used to retrieve data while the vector of siteNames is used for diagnostic plots created by the function, detrend.flow.

To best match up with water quality data, we recommend retrieving flow data for one year prior to the first year of water quality data. This allows for creating a time-averaged flow data set and not loose the first few months of water quality data due to lack of matching flow data. Data retrievals should also be made in light of the time needed by the USGS to review and approve their flow records.

Information about the settings for dvAvgWinSel, dvAvgWgtSel, dvAvgSidesSel, lowess.f, span, and max.fill are described in the help file for this function.

library(baytrends)
# Define Function Inputs
usgsGageID    <- usgsGages$usgsGageID # all RIM stations
siteName      <- usgsGages$siteName   # all RIM stations
yearStart     <- 1983
yearEnd       <- 2017
dvAvgWinSel   <- c(1, 5, 10, 15, 20, 30, 40, 50, 60, 90, 120, 150, 180, 210)
dvAvgWgtSel   <- "uniform"
dvAvgSidesSel <- 1
lowess.f      <- 0.2

# Run detrended.flow function
flow.detrended <- detrended.flow(usgsGageID, siteName, yearStart, yearEnd
                                 , dvAvgWinSel, dvAvgWgtSel, dvAvgSidesSel
                                 , lowess.f)

# Save list to data file separate use
save(flow.detrended, file='mySeasonallyDetrendedFlow.rda')

01578310-Susquehanna River at Conowingo, MD

Figure 1. Daily Flow [cms] as a Function of Day of Year. (Solid red line is the predicted GAM model using the model: gam(LogQ ~ s(doy,bs=‘cc’)).)

Figure 2. 1-day-smoothed Seasonally-adjusted Log Flow Residuals as a Function of Day of Year.

Figure 3. 10-day-smoothed Seasonally-adjusted Log Flow Residuals as a Function of Day of Year.

01578310-Susquehanna River at Conowingo, MD (Avg. Per.: d1)

Figure: d1: Mean, standard deviation, and number of observations as a Function of Day of Year.

01578310-Susquehanna River at Conowingo, MD (Avg. Per.: d10)

Figure: d10: Mean, standard deviation, and number of observations as a Function of Day of Year.

Seasonally Detrended Salinity

The function, detrended.salinity, returns a list of seasonally detrended salinity and companion statistics.

It is the user’s responsibility to save the list that is returned from the detrended.salinity function as salinity.detrended for integration with baytrends.

The detrended.salinity evaluates a user-supplied data set. baytrends includes a sample salinity data set, sal, that corresponds to the water quality data stored in dataCensored. Table 2 and Table 3 display some sample data and the structure of sal.

Table 2. Example Salinity Data

Table 3. Data Structure Associated with sal

#> 'data.frame':    51092 obs. of  4 variables:
#>  $ station : chr  "LE3.6" "LE3.6" "LE3.6" "LE3.6" ...
#>  $ date    : POSIXct, format: "1984-06-27" "1984-06-27" ...
#>  $ layer   : chr  "S" "M" "M" "B" ...
#>  $ salinity: num  11.7 12.1 12.6 13.4 21.6 ...

The below code chunk is an example for analyzing all of the data contained in the sample data frame, sal. As currently developed, detrended.salinity will average salinity in the surface (‘S’) and above pycnocline (‘AP’) layers to create an average salinity value for the surface and above pycnocline (SAP) layer. detrended.salinity will also average salinity in the bottom (‘B’) and below pycnocline (‘AP’) layers to create an average salinity value for the bottom and below pycnocline (BAP) layer.

The use and settings for dvAvgWinSel, lowess.f, minObs, and minObs.sd are described in the help file. The detrened.salinity function will assume the user wants to evaluate the salinity for all of the stations in sal and for all the years of available data.

# Define Function Inputs
df.sal      <- sal
dvAvgWinSel <- 30
lowess.f    <- 0.2
minObs      <- 40
minObs.sd   <- 10
                 
# Run Function
salinity.detrended <- detrended.salinity(df.sal, dvAvgWinSel, 
                                 lowess.f, minObs, minObs.sd)   

# Save list to data file separate use
save(salinity.detrended, file='mySeasonallyDetrendedSalinity.rda')

Example GAM Analysis

It is now possible to run gam4 using the results of the above analyses. The below code chunk will run the gam4 for surface total nitrogen at station CB5.4.

dfr       <- analysisOrganizeData(dataCensored, reports = NA)
df        <- dfr[["df"]]
analySpec <- dfr[["analySpec"]]

# down select analyzed GAM
analySpec$gamModels <- analySpec$gamModels[5]

stat = 'CB5.4'; dep = 'tn'; layer = 'S';
gamResult <- gamTest(df = df, dep = dep, stat = stat, layer = layer
                     , analySpec = analySpec
                     , salinity.detrended = salinity.detrended)
#> Warning: Detrended flow data (flow.detrended) not passed as argument in gamTest,
#> but models included in analySpec include those with 'flw_sal' term. Consider
#> passing detrended flow data or reducing models specified in analySpec.