Introduction

This document describes how to obtain seasonal information from the R package EGRET. 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 the seasonal information for a single year, or averages over several years, or in terms of flow normalized fluxes.

Getting started

First, you need to have installed and loaded the EGRET package. Then, you'll need need to create an eList object. See the EGRET vignette or user guide here for more information on how to use your own water quality data in EGRET. Once the eList object has been created, run the modelEstimation function to create a WRTDS model.

library(EGRET)
library(dplyr)
eList <- Choptank_eList

For this post, we will use the r eList$INFO$shortName, measuring r eList$INFO$paramShortName as an example.

To start with lets look at the results calculated for complete water years.

water_year <- tableResults(eList)
water_year_2010_flux <- water_year$`FN Flux [10^6kg/yr]`[water_year$Year == 2010]

Looking at the last column of these results we see that, for example, the flow normalized flux in water year 2010 is estimated to be r water_year_2010_flux 106 kg/year. Now, let's say we had a particular interest in the winter season which we define here as the months of December, January, and February.

The next step is to establish what season you are interested in looking at. All functions in EGRET can be done on the "water year" (Oct-Sept), the calendar year (Jan-Dec), or any set of sequential months. To define what period of analysis (PA) to use, there is a function setPA. The setPA function has two arguments:

For example let's say we want to consider the winter, defined here as December through February, the starting month (paStart) would be 12, and the length (paLong) would be 3:

eList <- setPA(eList, paStart = 12, paLong = 3)

Now lets view our results when we are focusing on the winter season.

winter <- tableResults(eList)
winter_2010_flux <- winter$`FN Flux [10^6kg/yr]`[winter$Year == 2010]

filter_setPA <- function(Daily, paStart, paLong){

  monthsToUse <- seq(from = paStart, 
                     by = 1,
                     length.out = paLong)
  monthsToUse[monthsToUse > 12] <- monthsToUse[monthsToUse > 12] - 12 

  Daily_filtered <- Daily %>% 
    filter(Month %in% monthsToUse) %>% 
    mutate(YearIndex = trunc(DecYear))

  crossesYear <- paLong + (paStart - 1) > 12

  if(crossesYear){
    Daily_filtered$YearIndex[Daily_filtered$Month >= paStart] <-
      Daily_filtered$YearIndex[Daily_filtered$Month >= paStart] + 1
  }

  get_years <- Daily_filtered %>% 
    group_by(YearIndex) %>% 
    mutate(Year = trunc(mean(DecYear))) %>% 
    ungroup() 


  return(get_years)

}

number_of_days <- eList$Daily %>% 
  filter_setPA(paStart = eList$INFO$paStart,
               paLong = eList$INFO$paLong) %>% 
  group_by(Year) %>% 
  summarise(n_days = sum(!is.na(Q)))

number_of_days_2010 <- number_of_days$n_days[number_of_days$Year == 2010]

Note that now the estimated flow normalized flux is r winter_2010_flux 106 kg/year. Let's take a moment to think about that in relation to the previous result. What it is saying is that the flux (mass per unit time) is greater during this three month period that the average flux for the entire water year. That makes sense, some seasons will have higher than average fluxes and other seasons will have lower than average fluxes.

Now, it might be that we want to think about these results, not in terms of a flux but in terms of mass alone. In the first result (water year) the mass for the year is r water_year_2010_flux 106 kg. But for the winter season we would compute the mass as r winter_2010_flux * r number_of_days_2010 / 365 to give us a result of r signif(winter_2010_flux*number_of_days_2010/365, digits = 3) 106 kg. To get this result we took the annual rate (r winter_2010_flux) and divided by the number of days in the year (365) to get a rate in mass per day and then multiplied by the number of days in the season (r number_of_days_2010) which is the sum of the length of those three months to simply get a total mass.

We have taken pains here to run through this calculation because users have sometimes been confused, wondering how the flux for a part of the year can be greater than the flux for the whole year. Depending on the ultimate objective of the analysis one might want to present the seasonal results in either of these ways (mass or mass per unit time).

Next we will look at descriptions of change.

Seasonal Changes

Let's use the tableChange function to explore the change from 1990 to 2010. We will do it first for the full water year and then for the winter season.

eList <- setPA(eList, paStart = 12, paLong = 12)
change_annual_flux <- tableChangeSingle(eList, flux = TRUE,
                                   yearPoints = c(1990,2010))
change_annual_conc <- tableChangeSingle(eList, flux = FALSE,
                                   yearPoints = c(1990,2010))

eList <- setPA(eList, paStart = 12, paLong = 3)
change_season_flux <- tableChangeSingle(eList, flux = TRUE,
                                   yearPoints = c(1990,2010))
change_season_conc <- tableChangeSingle(eList, flux = FALSE,
                                   yearPoints = c(1990,2010))

What we see is a fairly large difference between the concentration trends for the winter as compared to the whole year (concentration rose more steeply for the winter than it did for the year on average). But what we are going to focus on here is the trend in flux. The first thing to note is that the change from 1990 to 2010 is identical for the winter season and the year as a whole. That is, the change in the flow normalized flux for the full water year is r change_annual_flux[["change [10^6kg/yr]"]] 106 kg/year (it went from r water_year[["FN Flux [10^6kg/yr]"]][water_year$Year == 1990] to r water_year[["FN Flux [10^6kg/yr]"]][water_year$Year == 2010]) and then, when we look at the winter season the change is also r change_season_flux[["change [10^6kg/yr]"]] 106 kg/yr (from r winter[["FN Flux [10^6kg/yr]"]][winter$Year == 1990] to r winter[["FN Flux [10^6kg/yr]"]][winter$Year == 2010]). So the change for this season over the 20 year period is essentially the same as the change for the entire water year. In other words, the results tell us that although the change in flux (mass per unit time) for the winter is the same as for the full year, the change in the mass for the winter season is about 25% of the change for the full year (because the winter consists of about 25% of the days in the year). Thus, we can conclude that the winter change is not atypical of the changes for the other parts of the year.

The results shown above also express the change as a slope (either r format(change_annual_flux[["slope [10^6kg/yr/yr]"]]) or r format(change_season_flux[["slope [10^6kg/yr/yr]"]]) virtually identical to each other) and these are simply the change results divided by the number of years. The next entry in the tableChange output is for change expressed as %. Here we see a big difference between the winter and the whole year. The whole year shows an increase of r change_annual_flux[["change[percent]"]] % over the 20 years, while the winter season shows an increase of r change_season_flux[["change[percent]"]] %. That is because the same amount of increase r change_season_flux[["change [10^6kg/yr]"]] 106 kg / year is being compared to a the smaller number (1990 flow normalized annual flux of r water_year[["FN Flux [10^6kg/yr]"]][water_year$Year == 1990] 106 kg/year) in the first table and compared to a larger number (seasonal flux of r winter[["FN Flux [10^6kg/yr]"]][winter$Year == 1990] 106 kg/year) in the second table. So, even though the change is focused equally in the winter and non-winter months, the percentage change for the winter is smaller than the percentage change for the whole year.

Seasonal Load Fraction

Next, we can think about the seasonal load fraction.

You will need to read in two new functions 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. The functions use the package dplyr, a package that is useful for general data exploration and manipulation.

setupSeasons <- function(eList){

  Daily <- eList$Daily

  paLong <- eList$INFO$paLong
  paStart <- eList$INFO$paStart

  Daily_season <- filter_setPA(Daily,
                               paStart = paStart,
                               paLong = paLong)
  Daily_annual <- filter_setPA(Daily,
                               paStart = paStart,
                               paLong = 12)

  #Cleanup units:
  divideBy <- 1000000
  divideBy <- 1

  # Convert flux to kg/year
  unit_scale <- fluxConst[[9]]@unitFactor / divideBy

  SeasonResults <- Daily_season %>% 
      group_by(Year) %>% 
      summarize(DecYear = mean(DecYear),
                Counts = sum(!is.na(ConcDay)),
                Flux = mean(FluxDay) * unit_scale,
                FNFlux = mean(FNFlux) * unit_scale)

  AnnualResults <- Daily_annual %>%
    group_by(Year) %>% 
    summarize(AnnualCounts = sum(!is.na(ConcDay)),
              AnnualFlux = mean(FluxDay) * unit_scale,
              AnnualFNFlux = mean(FNFlux) * unit_scale) %>%
    filter(AnnualCounts >= 365)


  seasonPctResults <- SeasonResults %>%
    select(DecYear, Year, Flux, FNFlux, Counts) %>% 
    left_join(select(AnnualResults,
                     Year, AnnualCounts,
                     AnnualFlux, AnnualFNFlux), by="Year") %>%
      filter(!is.na(AnnualFlux)) %>%
    mutate(pctFlux = 100*Flux*Counts/(AnnualFlux*AnnualCounts),
           pctFNFlux = 100*FNFlux*Counts/(AnnualFNFlux*AnnualCounts)) %>%
    rename(SeasonFlux = Flux,
           SeasonFNFlux = FNFlux)

  return(seasonPctResults)
}

Simply use the loaded eList to calculate these seasonal load fractions. Let's go back to the winter season (Dec-Feb):

eList <- setPA(eList, paStart = 12, paLong = 3)
seasonPctResults <- setupSeasons(eList)

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|Calendar Year of mid-date of the year| |AnnualFlux|Estimated flux for the water year in millions of kg| |AnnualFNFlux|Flow Normalized flux for the water year in millions of kg| |SeasonFlux|Estimated flux for the season in millions of kg| |SeasonFNFlux|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|

Plotting the time series

We can make a graph showing the percentage flux (estimated annual and flow normalized). Note, this workflow uses base R plotting functions. You could also use the EGRET function genericEGRETDotPlot to automatically to pick some plotting styles that are consistent with other EGRET plots.

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") 
par(mar=c(5,6,4,2) + 0.1,mgp=c(3,1,0))
plot(seasonPctResults$DecYear, seasonPctResults$pctFlux,pch=20,
     yaxs="i",ylim = c(0,100),las=1,tck=.01,
     xlab="Year",ylab="Percentage of Annual Flux",
     main=plotTitle,cex=1.5)
lines(seasonPctResults$DecYear,seasonPctResults$pctFNFlux,col="green",lwd=2)
axis(3, labels = FALSE,tck=.01)
axis(4, labels = FALSE,tck=.01)

We can interpret this example graph as follows. The winter flux of r eList$INFO$paramShortName 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.

Keep in mind, the way we are defining "year" is what year the ending year of the period of anaylsis fell. So, for this analysis, the full 2010 "year" is from Dec. 2009 through the end of November 2010.

years00_10 <- filter(seasonPctResults, Year >= 2000, Year <= 2010)

sumYears <- sum(years00_10$AnnualFlux)

sumSeasons <- sum(years00_10$SeasonFlux * years00_10$Counts /
                    years00_10$AnnualCounts)

avePct <- 100 * sumSeasons / sumYears

The total flux for all years in the period of interest in millions of kg is sumYears = r signif(sumYears, digits = 3).

The total seasonal flux for all years of the period of interest in millions of kg is sumSeasons = r signif(sumSeasons, digits = 3).

The percentage of the total flux for the years 2000 through 2010 that was transported in the winter months is avePct = r signif(avePct, digits = 3).

This can be determined for any set of years simply by changing Year values in the dplyr::filter function. So, for the years 1990-1999:

years90_99 <- filter(seasonPctResults, Year >= 1990, Year <= 1999)

c("sumYears" = sum(years90_99$AnnualFlux),
"sumSeasons" = sum(years90_99$SeasonFlux * years90_99$Counts / years90_99$AnnualCounts),
"avePct" = 100 * sumSeasons / sumYears)


USGS-R/EGRET documentation built on Feb. 9, 2024, 5:30 p.m.