<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/seasonal-analysis-in-egret/"></div></div><div><a data-url="https://www.r-bloggers.com/seasonal-analysis-in-egret/" data-counturl="https://www.r-bloggers.com/seasonal-analysis-in-egret/" data-count="horizontal" data-text="Seasonal Analysis in EGRET" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://owi.usgs.gov/blog/seasonal-analysis/"> The USGS OWI blog </a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<h1>Introduction</h1>
<p>This document describes how to obtain seasonal information from the R package <a target="_blank" rel="nofollow" href="https://cran.r-project.org/package=EGRET">EGRET</a> – Exploration and Graphics for RivEr Trends: An R-package for the analysis of long-term changes in water quality and streamflow, including the water-quality method Weighted Regressions on Time, Discharge, and Season (WRTDS)</p>
<p>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.</p>
<h1>Getting started</h1>
<p>First, you need to have installed and loaded the <code>EGRET</code> package. Then, you’ll need need to create an <code>eList</code> object. See the <code>EGRET</code> vignette or user guide <a target="_blank" rel="nofollow" href="http://pubs.usgs.gov/tm/04/a10/">here</a> for more information on how to use your own water quality data in <code>EGRET</code>. Once the <code>eList</code> object has been created, run the <code>modelEstimation</code> function to create a WRTDS model.</p>
<p>For this post, we will use the Choptank River as an example. There is an example data set included in <code>EGRET</code> with measured nitrate in the Choptank River. Here is an example of how to load that example data:</p>
<pre><code>library(EGRET)
eList &lt;- Choptank_eList
</code></pre>
<p>To start with lets look at the results calculated for complete water years.</p>
<pre><code>tableResults(eList)
</code></pre>
<pre><code>## 
##    Choptank River 
##    Inorganic nitrogen (nitrate and nitrite)
##    Water Year 
## 
##    Year   Discharge    Conc    FN_Conc     Flux    FN_Flux
##              cms            mg/L             10^6 kg/yr 
## 
##    1980      4.25     0.949     1.003    0.1154     0.106
##    1981      2.22     1.035     0.999    0.0675     0.108
##    1982      3.05     1.036     0.993    0.0985     0.110
##    1983      4.99     1.007     0.993    0.1329     0.112
....
##    2006      3.59     1.382     1.362    0.1409     0.147
##    2007      4.28     1.408     1.382    0.1593     0.149
##    2008      2.56     1.477     1.401    0.1008     0.149
##    2009      3.68     1.409     1.419    0.1328     0.149
##    2010      7.19     1.323     1.438    0.2236     0.149
##    2011      5.24     1.438     1.457    0.1554     0.148
</code></pre>
<p>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 0.149 10<sup>6</sup> 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. <em>Note</em>: some lines (1984-2005) were removed from this blog post simply to save space.</p>
<p>The next step is to establish what season you are interested in looking at. All functions in <code>EGRET</code> 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 <code>setPA</code>. The <code>setPA</code> function has two arguments:</p>
<ul>
<li><code>paStart</code> is the number of the calendar month that is the start of the season.</li>
<li><code>paLong</code> is the length of the season in months (it can be any number from 1 to 12).</li>
</ul>
<p>For example let’s say we want to consider the winter, defined here as December through February, the starting month (<code>paStart</code>) would be 12, and the length (<code>paLong</code>) would be 3:</p>
<pre><code>eList &lt;- setPA(eList, paStart = 12, paLong = 3)
</code></pre>
<p>Now lets view our results when we are focusing on the winter season.</p>
<pre><code>tableResults(eList)
</code></pre>
<pre><code>## 
##    Choptank River 
##    Inorganic nitrogen (nitrate and nitrite)
##    Season Consisting of Dec Jan Feb 
## 
##    Year   Discharge    Conc    FN_Conc     Flux    FN_Flux
##              cms            mg/L             10^6 kg/yr 
## 
##    1980     4.220      1.10      1.11    0.1403     0.156
##    1981     1.960      1.20      1.13    0.0735     0.161
##    1982     5.057      1.16      1.15    0.1764     0.165
##    1983     3.504      1.21      1.16    0.1310     0.169
....
##    2006     5.843      1.51      1.51    0.2675     0.214
##    2007     5.417      1.55      1.53    0.2419     0.216
##    2008     2.436      1.62      1.55    0.1217     0.217
##    2009     2.711      1.66      1.56    0.1396     0.216
##    2010    13.779      1.26      1.57    0.4161     0.215
##    2011     3.369      1.66      1.57    0.1669     0.215
</code></pre>
<p>Note that now the estimated flow normalized flux is 0.215 10<sup>6</sup> 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.</p>
<p>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 0.149 10<sup>6</sup> kg. But for the winter season we would compute the mass as 0.215 * 90 / 365 to give us a result of 0.053 10<sup>6</sup> kg. To get this result we took the annual rate (0.215) 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 (90) which is the sum of the length of those three months (31 + 31 + 28) to simply get a total mass.</p>
<p>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).</p>
<p>Next we will look at descriptions of change.</p>
<h1>Seasonal Changes</h1>
<p>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.</p>
<pre><code>eList &lt;- setPA(eList, paStart = 10, paLong = 12)
tableChange(eList, yearPoints = c(1990,2010))
</code></pre>
<pre><code>## 
##    Choptank River 
##    Inorganic nitrogen (nitrate and nitrite)
##    Water Year 
## 
##            Concentration trends
##    time span       change     slope    change     slope
##                      mg/L   mg/L/yr        %       %/yr
## 
##  1990  to  2010      0.31     0.016        28       1.4
## 
## 
##                  Flux Trends
##    time span          change        slope       change        slope
##                   10^6 kg/yr    10^6 kg/yr /yr      %         %/yr
##  1990  to  2010         0.02      0.00099           15         0.76
</code></pre>
<pre><code>eList &lt;- setPA(eList, paStart = 12, paLong = 3)
tableChange(eList, yearPoints = c(1990,2010))
</code></pre>
<pre><code>## 
##    Choptank River 
##    Inorganic nitrogen (nitrate and nitrite)
##    Season Consisting of Dec Jan Feb 
## 
##            Concentration trends
##    time span       change     slope    change     slope
##                      mg/L   mg/L/yr        %       %/yr
## 
##  1990  to  2010      0.24     0.012        18      0.89
## 
## 
##                  Flux Trends
##    time span          change        slope       change        slope
##                   10^6 kg/yr    10^6 kg/yr /yr      %         %/yr
##  1990  to  2010         0.02        0.001           10         0.51
</code></pre>
<p>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 +0.020 10<sup>6</sup> kg/year (it went from 0.129 to 0.149) and then, when we look at the winter season the change is also +0.020 10<sup>6</sup> kg/yr (from 0.195 to 0.215). 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.</p>
<p>The results show above also express the change as a slope (either 0.001 or 0.00099 virtually identical to each other) and these are simply the change results divided by the number of years. The next entry in the <code>tableChange</code> 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 + 15 % over the 20 years, while the winter season shows an increase of 10 %. That is because the same amount of increase 0.02 10<sup>6</sup> kg / year is being compared to a the smaller number (1990 flow normalized annual flux of 0.129 10<sup>6</sup> kg/year) in the first table and compared to a larger number (seasonal flux of 0.195 10<sup>6</sup> 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.</p>
<h1>Seasonal Load Fraction</h1>
<p>Next, we can think about the seasonal load fraction.</p>
<p>You will need to read in two new function called <code>setupSeasons</code> and <code>setupYearsPlus</code> 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 <a target="_blank" rel="nofollow" href="https://cran.r-project.org/package=dplyr"><code>dplyr</code></a>, a package that is useful for general data exploration and manipulation.</p>
<pre><code>library(dplyr)

setupSeasons &lt;- function(eList){
  Daily &lt;- eList$Daily
  
  SeasonResults &lt;- setupYearsPlus(Daily, 
                                  paLong = eList$INFO$paLong, 
                                  paStart = eList$INFO$paStart)
  AnnualResults &lt;- setupYearsPlus(Daily, 
                                  paLong = 12, 
                                  paStart = eList$INFO$paStart) %>%
    filter(Counts >= 365) #To make sure we are getting full years
  
  divideBy &lt;- 1000000
  
  annualPctResults &lt;- AnnualResults %>%
    mutate(FluxYear = Flux*Counts/divideBy,
           FNFluxYear = FNFlux*Counts/divideBy) %>%
    select(FluxYear, FNFluxYear, Year)
  
  seasonPctResults &lt;- SeasonResults %>%
    mutate(FluxSeason = Flux*Counts/divideBy,
           FNFluxSeason = FNFlux*Counts/divideBy) %>%
    left_join(annualPctResults, by=&quot;Year&quot;) %>%
    mutate(pctFlux = 100*FluxSeason/FluxYear,
           pctFNFlux = 100*FNFluxSeason/FNFluxYear) %>%
    select(-Q, -Conc, -Flux, -FNFlux, -FNConc, -Counts) %>%
    rename(seasonStart = paStart,
           seasonLong = paLong)
  
  return(seasonPctResults)
}

setupYearsPlus &lt;- function (Daily, paLong = 12, paStart = 10){

  monthsToUse &lt;- seq(paStart, length=paLong)
  monthsToUse[monthsToUse > 12] &lt;- monthsToUse[monthsToUse > 12] - 12
  
  crossesYear &lt;- paLong + (paStart - 1) > 12
  
  AnnualResults &lt;- Daily %>%
    mutate(Year =  as.integer(format(Date, &quot;%Y&quot;))) %>%
    filter(Month %in% monthsToUse) %>%
    mutate(Year = if(crossesYear){
      ifelse(Month >= paStart, Year + 1, Year)
    } else {
      Year
    }) %>%
    group_by(Year) %>%
      summarise(DecYear = mean(DecYear, na.rm = TRUE),
                Q = mean(Q, na.rm = TRUE),
                Conc = mean(ConcDay, na.rm = TRUE),
                Flux = mean(FluxDay, na.rm = TRUE),
                FNConc = mean(FNConc, na.rm = TRUE),
                FNFlux = mean(FNFlux, na.rm = TRUE),
                Counts = sum(!is.na(ConcDay))) %>%
    mutate(paLong = paLong,
           paStart = paStart) 
      
  return(AnnualResults)
  
}
</code></pre>
<p>Simply use the loaded <code>eList</code> to calculate these seasonal load fractions. Let’s go back to the winter season (Dec-Feb):</p>
<pre><code>eList &lt;- setPA(eList, paStart = 12, paLong = 3)
seasonPctResults &lt;- setupSeasons(eList)
</code></pre>
<h1>Looking at your results</h1>
<p>What you now have is a data frame called <code>seasonPctResults</code>. The columns it contains are the following:</p>
<table>
<colgroup>
<col width="19%">
<col width="73%">
</colgroup>
<thead>
<tr>
<th>variable</th>
<th>Definition</th>
</tr>
</thead>
<tbody>
<tr>
<td>DecYear</td>
<td>Decimal Year of the mid-date of the season</td>
</tr>
<tr>
<td>Year</td>
<td>Calendary Year of mid-date of the year</td>
</tr>
<tr>
<td>FluxYear</td>
<td>Estimated flux for the year in millions of kg</td>
</tr>
<tr>
<td>FNFluxYear</td>
<td>Flow Normalized flux for the year in millions of kg</td>
</tr>
<tr>
<td>FluxSeason</td>
<td>Estimated flux for the season in millions of kg</td>
</tr>
<tr>
<td>FNFluxSeason</td>
<td>Flow Normalized flux for the season in millions of kg</td>
</tr>
<tr>
<td>pctFlux</td>
<td>Season flux as a percentage of Annual Flux</td>
</tr>
<tr>
<td>pctFNFlux</td>
<td>FlowNormalized Seasonal Flux as a percent of Flow Normalized Annual Flux</td>
</tr>
<tr>
<td>seasonLong</td>
<td>Length of the Season in Months</td>
</tr>
<tr>
<td>seasonStart</td>
<td>Starting Month of the Season, 1=January</td>
</tr>
</tbody>
</table>
<h1>Plotting the time series</h1>
<p>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 <code>EGRET</code> function <code>genericEGRETDotPlot</code> to automatically to pick some plotting styles that are consistent with other <code>EGRET</code> plots.</p>
<pre><code>plotTitle = paste(&quot;Seasonal Flux as a Percent of Annual Flux\n&quot;,
                  eList$INFO$shortName, eList$INFO$paramShortName,
                  &quot;\nSolid line is percentage of flow normalized flux&quot;) 
par(mar=c(5,6,4,2) + 0.1,mgp=c(3,1,0))
plot(seasonPctResults$DecYear, seasonPctResults$pctFlux,pch=20,
     yaxs=&quot;i&quot;,ylim = c(0,100),las=1,tck=.01,
     xlab=&quot;Year&quot;,ylab=&quot;Percentage of Annual Flux&quot;,
     main=plotTitle,cex=1.5)
lines(seasonPctResults$DecYear,seasonPctResults$pctFNFlux,col=&quot;green&quot;,lwd=2)
axis(3, labels = FALSE,tck=.01)
axis(4, labels = FALSE,tck=.01)
</code></pre>
<p><img title="Seasonal flux as a percentage of annual flux." alt="Seasonal flux as a percentage of annual flux." src="https://i2.wp.com/owi.usgs.gov/blog/static/seasonal-analysis/unnamed-chunk-8-1.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>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.</p>
<h1>Computing averages over a period of years</h1>
<p>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.</p>
<p>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.</p>
<ul>
<li>
<p>Filter the data frame <code>seasonPctResults</code> for the years 2000 – 2010.</p>
</li>
<li>
<p>Now we can compute the sum of the annual fluxes 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.</p>
</li>
</ul>
<pre><code>years00_10 &lt;- filter(seasonPctResults, Year >= 2000, Year &lt;= 2010)

sumYears &lt;- sum(years00_10$FluxYear)
 
sumSeasons &lt;- sum(years00_10$FluxSeason)

avePct &lt;- 100 * sumSeasons / sumYears
</code></pre>
<p>The total flux for all years in the period of interest in millions of kg is <code>sumYears</code> = 1.7297595.</p>
<p>The total seasonal flux for all years of the period of interest in millions of kg is <code>sumSeasons</code> = 0.6099614.</p>
<p>The percentage of the total flux for the years 2000 through 2010 that was transported in the winter months is <code>avePct</code> = 35.2627852.</p>
<p>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.</p>
<h2>Questions</h2>
<ul>
<li>
<p>Robert M. Hirsch <a target="_blank" rel="nofollow" href="mailto:rhirsch@usgs.gov"><i></i></a> <a target="_blank" rel="nofollow" href="https://scholar.google.com/citations?user=Jt5I-0gAAAAJ"><i></i></a> <a target="_blank" rel="nofollow" href="https://www.researchgate.net/profile/Robert_Hirsch3"><i></i></a> <a target="_blank" rel="nofollow" href="https://www.usgs.gov/staff-profiles/robert-hirsch"><i></i></a></p>
</li>
<li>
<p>Laura DeCicco <a target="_blank" rel="nofollow" href="mailto:ldecicco@usgs.gov"><i></i></a> <a target="_blank" rel="nofollow" href="https://twitter.com/DeCiccoDonk"><i></i></a> <a target="_blank" rel="nofollow" href="https://github.com/ldecicco-usgs"><i></i></a> <a target="_blank" rel="nofollow" href="https://scholar.google.com/citations?hl=en&user=jXd0feEAAAAJ"><i></i></a> <a target="_blank" rel="nofollow" href="https://www.researchgate.net/profile/Laura_De_Cicco"><i></i></a> <a target="_blank" rel="nofollow" href="https://www.usgs.gov/staff-profiles/laura-decicco"><i></i></a></p>
</li>
</ul>
<p>Information on USGS-R packages used in this post:</p>
<table>
<colgroup>
<col width="19%">
<col width="73%">
</colgroup>
<tbody>
<tr>
<td><a target="_blank" rel="nofollow" href="https://github.com/USGS-R/EGRET"><img alt="USGS-R image icon" src="https://i1.wp.com/owi.usgs.gov/blog/images/USGS_R.png?w=456&ssl=1" data-recalc-dims="1"></a></td>
<td><a target="_blank" rel="nofollow" href="https://github.com/USGS-R/EGRET/issues">EGRET</a>: Exploration and Graphics for RivEr Trends: An R-package for the analysis of long-term changes in water quality and streamflow, including the water-quality method Weighted Regressions on Time, Discharge, and Season (WRTDS)</td>
</tr>
<tr>
<td><a target="_blank" rel="nofollow" href="https://github.com/USGS-R/dataRetrieval"><img alt="USGS-R image icon" src="https://i1.wp.com/owi.usgs.gov/blog/images/USGS_R.png?w=456&ssl=1" data-recalc-dims="1"></a></td>
<td><a target="_blank" rel="nofollow" href="https://github.com/USGS-R/dataRetrieval/issues">dataRetrieval</a>: This R package is designed to obtain USGS or EPA water quality sample data, streamflow data, and metadata directly from web services.</td>
</tr>
</tbody>
</table>
<p><div><div><div><div data-width="55" data-show-faces="false" data-send="true" data-layout="box_count" data-height="62" data-href="https://www.r-bloggers.com/seasonal-analysis-in-egret/"></div></div><div><a data-url="https://www.r-bloggers.com/seasonal-analysis-in-egret/" data-counturl="https://www.r-bloggers.com/seasonal-analysis-in-egret/" data-count="vertical" data-text="Seasonal Analysis in EGRET" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://owi.usgs.gov/blog/seasonal-analysis/"> The USGS OWI blog </a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/?s=data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/?s=Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/?s=map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/?s=animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/?s=RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/?s=sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/?s=LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/?s=SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/?s=eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/?s=git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/?s=hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/?s=Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/?s=PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/?s=time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/?s=trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/Ka23_weSWh4"></a></div>