In USGS-R/EGRETextra: Exploration and Graphics for RivEr Trends (EGRET) Extra Code

library(rmarkdown)
options(continue=" ")
options(width=60)
library(knitr)
library(EGRET)

Introduction

This document describes how to obtain information from EGRET results that describe the seasonal distribution of fluxes. For example, we might want to know the fraction of the load that takes place in the winter season (say that is December, January, and February). We can look at it for a single year, we can look at averages of it over several years, or we can look at it in terms of flow normalized fluxes.

Getting started

First, you need to have loaded the EGRET package and you need to have run the modelEstimation function and as a result of that, have an eList object.

Next, you will need to read in two new function called setupSeasons and setupYearsPlus designed for this purpose. You can copy them from here and paste them into your workspace (all as a single copy and paste) or you can create an .R file from them that you will source each time you want to use them.

setupSeasons <- function(localDaily, paLong, paStart){
  SeasonResults <- setupYearsPlus(localDaily, paLong = paLong, paStart = paStart)
  AnnualResults <- setupYearsPlus(localDaily, paLong = 12, paStart = paStart)
  numYears <- length(AnnualResults$DecYear)
  divide <- 1000000
  DecYear <- AnnualResults$DecYear
  Year <- trunc(DecYear)
  FluxYear <- AnnualResults$Flux*AnnualResults$Counts/divide
  FNFluxYear <- AnnualResults$FNFlux*AnnualResults$Counts/divide
  FluxSeason <- SeasonResults$Flux[1:numYears]*SeasonResults$Counts[1:numYears]/divide
  FNFluxSeason <- SeasonResults$FNFlux[1:numYears]*SeasonResults$Counts[1:numYears]/divide
  pctFlux <- ifelse(is.na(FluxYear)|is.na(FluxSeason),NA,100*FluxSeason/FluxYear)
  pctFNFlux <- ifelse(is.na(FNFluxYear)|is.na(FNFluxSeason),NA,100*FNFluxSeason/FNFluxYear)
  seasonPctResults <- data.frame(DecYear,Year,FluxYear,FNFluxYear,FluxSeason,FNFluxSeason,pctFlux,pctFNFlux)
  seasonLong <- rep(paLong,numYears)
  seasonStart <- rep(paStart,numYears)
  seasonPctResults <- data.frame(seasonPctResults,seasonLong,seasonStart)
  return(seasonPctResults)
}

setupYearsPlus <- function (localDaily, paLong = 12, paStart = 10){

# This is an augmented version of setupYears 
#  that also returns the number of good days in each year or season

  numDays <- length(localDaily$MonthSeq)
  firstMonthSeq <- localDaily$MonthSeq[1]
  lastMonthSeq <- localDaily$MonthSeq[numDays]
  Starts <- seq(paStart, lastMonthSeq, 12)
  Ends <- Starts + paLong - 1
  StartEndSeq <- data.frame(Starts, Ends)
  StartEndSeq <- StartEndSeq[(StartEndSeq$Starts >= firstMonthSeq) & 
                               (StartEndSeq$Ends <= lastMonthSeq), ]
  firstMonth <- StartEndSeq[1, 1]
  numYears <- length(StartEndSeq$Starts)
  DecYear <- rep(NA, numYears)
  Q <- rep(NA, numYears)
  Conc <- rep(NA, numYears)
  Flux <- rep(NA, numYears)
  FNConc <- rep(NA, numYears)
  FNFlux <- rep(NA, numYears)
  Counts <- rep(NA, numYears)
  for (i in 1:numYears) {
    startMonth <- (i - 1) * 12 + firstMonth
    stopMonth <- startMonth + paLong - 1
    DailyYear <- localDaily[which(localDaily$MonthSeq %in% 
                                    startMonth:stopMonth), ]
    counter <- ifelse(is.na(DailyYear$ConcDay), 0, 1)
    if (length(counter) > 0) {
      good <- (sum(counter) > 25)
    }
    else {
      good <- FALSE
    }
    DecYear[i] <- mean(DailyYear$DecYear)
    Q[i] <- mean(DailyYear$Q)
    if (good) {
      Conc[i] <- mean(DailyYear$ConcDay, na.rm = TRUE)
      Flux[i] <- mean(DailyYear$FluxDay, na.rm = TRUE)
      FNConc[i] <- mean(DailyYear$FNConc, na.rm = TRUE)
      FNFlux[i] <- mean(DailyYear$FNFlux, na.rm = TRUE)
      Counts[i] <- sum(counter)
    }
  }
  PeriodStart <- rep(paStart, numYears)
  PeriodLong <- rep(paLong, numYears)
  AnnualResults <- data.frame(DecYear, Q, Conc, Flux, FNConc, 
                              FNFlux, PeriodLong, PeriodStart, Counts)
  return(AnnualResults)
}

The next step is to establish what season you are interested in looking at. We do this by specifying paStart and paLong.

paStart is the number of the calendar month that is the start of the season.
paLong is the length of the season in months (it can be any number from 1 to 12).

For example lets say we want to consider the winter, defined here as December through February. This code we would use is. This is written with the example data set Choptank_eList, which comes out of the EGRET package. In running this script you would delete the line eList <- Choptank_eList and enter the values of paLong and paStart that you wish to use.

library(EGRET)
eList <- Choptank_eList
Daily <- eList$Daily
seasonPctResults <- setupSeasons(Daily, paLong = 3, paStart = 12)

Looking at your results

What you now have is a data frame called seasonPctResults. The columns it contains are the following:

|variable| Definition| |:----|:----| |DecYear|Decimal Year of the mid-date of the season| |Year|Calendary Year of mid-date of the year| |FluxYear|Estimated flux for the year in millions of kg| |FNFluxYear|Flow Normalized flux for the year in millions of kg| |FluxSeason|Estimated flux for the season in millions of kg| |FNFluxSeason|Flow Normalized flux for the season in millions of kg| |pctFlux|Season flux as a percentage of Annual Flux| |pctFNFlux|FlowNormalized Seasonal Flux as a percent of Flow Normalized Annual Flux| |seasonLong|Length of the Season in Months| |seasonStart|Starting Month of the Season, 1=January

You can just print it out as a simple table:

seasonPctResults

Plotting the time series

We can make a graph showing the percentage flux (estimated annual and flow normalized)

nYears <- length(seasonPctResults$DecYear)
xlim <- c(seasonPctResults$DecYear[1]-1,seasonPctResults$DecYear[nYears]+1)
xTicks <- pretty(xlim)
ylim <- c(0,100)
yTicks <- seq(0,100,10)
plotTitle = paste("Seasonal Flux as a Percent of Annual Flux\n",eList$INFO$shortName,eList$INFO$paramShortName,"\nSolid line is percentage of flow normalized flux") 
genericEGRETDotPlot(seasonPctResults$DecYear,seasonPctResults$pctFlux,xlim=xlim,ylim=ylim,xTicks=xTicks,yTicks=yTicks,xaxs="i",yaxs="i",xlab="Year",ylab="Percentage of Annual Flux",plotTitle=plotTitle,xDate=TRUE,cex=1.5)
par(new=TRUE)
genericEGRETDotPlot(seasonPctResults$DecYear,seasonPctResults$pctFNFlux,xlim=xlim,ylim=ylim,xTicks=xTicks,yTicks=yTicks,xaxs="i",yaxs="i",xlab="",ylab="",plotTitle=plotTitle,xDate=TRUE,cex=1.5,type="l",col="green",lwd=2)

We can interpret this example graph as follows. The winter flux of nitrate fluctuates a good deal from year to year. From a low of around 10% to a high of around 60% but the mean percentage hasn't changed much over the years. It is around 35% of the annual total flux.

Computing averages over a period of years

Let's say we wanted to answer the question, what percentage of the annual total flux moved in the winter season during the years 2000 through 2010. We can answer that question with a simple set of calculations.

• First we need to look at the list of annual values that we printed out above and find the index numbers for the two years specified. The year 2000 is number 21 on the list and 2010 is number 31 on the list.

• Now we can compute the sum of the annual fluxs for those years and the sum of the seasonal fluxes for those years, and then get our answer by taking the ratio and multiplying by 100.

```r sumYears <- sum(seasonPctResults$FluxYear[21:31])

This is the total flux for all years

in the period of interest in millions of kg

sumYears

sumSeasons <- sum(seasonPctResults$FluxSeason[21:31])

This is the total seasonal flux for all years

of the period of interest in millions of kg

sumSeasons

avePct <- 100 * sumSeasons / sumYears

This is the percentage of the total flux for the

period of interest that was transported during the season of interest

avePct ````` This is the percentage of the total flux for the years 2000 through 2010 that was transported in the winter months.

This can be determined for any set of years simply by changing the two numbers inside the brackets to the index numbers of the first and last years of interest.

Disclaimer

Software created by USGS employees along with contractors and grantees (unless specific stipulations are made in a contract or grant award) are to be released as Public Domain and free of copyright or license. Contributions of software components such as specific algorithms to existing software licensed through a third party are encouraged, but those contributions should be annotated as freely available in the Public Domain wherever possible. If USGS software uses existing licensed components, those licenses must be adhered to and redistributed.

Although this software has been used by the U.S. Geological Survey (USGS), no warranty, expressed or implied, is made by the USGS or the U.S. Government as to accuracy and functionality, nor shall the fact of distribution constitute any such warranty, and no responsibility is assumed by the USGS in connection therewith.

USGS-R/EGRETextra documentation built on May 9, 2019, 6:08 p.m.