EX 5. Streamflow analysis: Annual water balance

In this exercise, you will: - read USGS streamflow data into R
- calculate monthly discharge
- calculate runoff coefficients

0. Load the libraries you’ll need.

The USGS has a website dedicated to R: https://owi.usgs.gov/R/ There are packages for streamflow and water quality data retrieval and analysis. Here, we install and use the dataRetrieval package from the USGS.

install.packages("dataRetrieval", 
                 repos=c("http://owi.usgs.gov/R",
                         getOption("repos")))

library(zoo)
library(raster)
library(rgdal)
library(hydroTSM)
library(dataRetrieval)  # USGS package that gets streamflow data direct from the USGS website

See what dataRetrieval can do, using “vignette”:

vignette("dataRetrieval",package = "dataRetrieval")

II. Load and summarize rainfall data

A. Load the daily time series for Poway station.

See Exercise 4. Subset to only post-Oct-1955 (see Ex4).

B. Summarize preciptation data to annual values by water year and plot

See Exercise 4… Look at your data. Is WaterYear correct?

Aggregate the data to annual (water year) total precipitation (see Ex 4) and plot to make sure it worked:

Make sure your column names are “WaterYear” and “Rainfall.mm”.

II. Load and summarize streamflow data

A. Load streamflow at the San Diego River, Mast Blvd.

1. Specify what station you’re interested in.

  # Try the code below with the site.code here, then use the site code for your watershed.
site.code = "11023340"  #  The USGS streamgage code for Los Penasquitos

2. See what data is available.

readNWISsite(site.code)  # Note:  doing readNWISsite("11022480") gives the same result.
what.data = whatNWISdata(siteNumber = site.code)
what.data[1:10,]  # just look a first 10 records

3. Read in daily discharge data direct from the USGS website.

parameter.code = "00060"  # this is the code for stream discharge.
start.date = ""  # Blanks get all of the data
end.date = ""
x = readNWISdv(site.code,parameter.code,start.date,end.date)
head(x)

##   agency_cd  site_no       Date X_00060_00003 X_00060_00003_cd
## 1      USGS 11023340 1964-10-01           0.5                A
## 2      USGS 11023340 1964-10-02           0.5                A
## 3      USGS 11023340 1964-10-03           0.5                A
## 4      USGS 11023340 1964-10-04           0.6                A
## 5      USGS 11023340 1964-10-05           0.6                A
## 6      USGS 11023340 1964-10-06           0.6                A

4. Rename the columns of x.

The names for the discharge and QA columns aren’t very nice, so rename them:

names(x)[c(4,5)] = c("Q.ft.s","QA")
head(x)

##   agency_cd  site_no       Date Q.ft.s QA
## 1      USGS 11023340 1964-10-01    0.5  A
## 2      USGS 11023340 1964-10-02    0.5  A
## 3      USGS 11023340 1964-10-03    0.5  A
## 4      USGS 11023340 1964-10-04    0.6  A
## 5      USGS 11023340 1964-10-05    0.6  A
## 6      USGS 11023340 1964-10-06    0.6  A

5. Tabulate the quality infomation.

The quality codes can be found by googling “usgs streamflow data quality flags”, which points you to http://waterdata.usgs.gov/usa/nwis/uv?codes_help.

Scroll down the webiste to “Daily Value Qualification Code” to see what the quality flags mean.

table(x$QA)  #  Creates a table of the counts of each quality flag value.
#  Based on the quality flags and this table, can you use all of the data?

B. Calculate the depth of annual discharge in mm/yr.

1. Use the aggregate function as you did in Exercise 4 to calculate total rainfall by water year.

Be sure to rename the columns of the data frame that has the annual mean values. ** Also be sure to convert the year and mean to numeric, if they are factors. **

  # What from and to dates should you use for your dataset?
breaks <- seq(from=as.Date("1965-10-01"), to=as.Date("2023-10-01"), by="year")  
years.breaks = as.numeric(format(breaks,"%Y"))
labels.wy = years.breaks[2:length(breaks)]
x$WaterYear <- cut(x$Date, breaks,labels=labels.wy)

data.avail.by.wy = aggregate(x$Q.ft.s,by=list(x$WaterYear),FUN=length)
annual.mean.cfs = aggregate(x$Q.ft.s,by=list(x$WaterYear),FUN=mean)
names(annual.mean.cfs) = c("WYear","MeanQcfs")

# The output of aggregate is sometimes a factor, so convert both Year and MeanQ to numeric.
annual.mean.cfs$WYear = as.numeric(as.character(annual.mean.cfs$WYear))
annual.mean.cfs$MeanQcfs = as.numeric(as.character(annual.mean.cfs$MeanQcfs))
head(annual.mean.cfs)

##   WYear MeanQcfs
## 1  1966 8.196438
## 2  1967 6.814000
## 3  1968 1.462350
## 4  1969 6.770712
## 5  1970 2.001315
## 6  1971 2.317260

2. Calculate annual runoff in mm/year.

For the calculation, you need the drainage area of your watershed.

site.data = readNWISsite(site.code)
drain.area.mi2 = site.data$drain_area_va
drain.area.mi2

## [1] 42.1

Now calculate annual discharge from your watershed in mm by doing units conversions on annual.mean.cfs. (write out the units conversion first on paper!). I’m not showing you the code…you have to write it.

Call this new annual discharge data frame (in mm), and call it “Q.ann.df”. Use the data.frame command to make sure the result is a data.frame class. Make sure it has column names of “WaterYear” and “Q.ann.mm”.

III. Annual rainfall-runoff relationships.

A. Plot annual rainfall vs annual runoff:

First, you need to merge the datasets, so you have a data frame with only the years that have data for both.

Q.P.merge = merge(Q.ann.df,x.annual.wy,by.x="WaterYear",by.y="WaterYear")
#  Print Q.P.merge to make sure it worked.

plot(Q.P.merge$Rainfall.mm,Q.P.merge$Q.ann.mm,xlab="Annual precipitation, mm",ylab="Annual runoff, mm")
#  Add a 1:1 line to see if runoff is ever greater than precipitation.
abline(0,1)  #  Run ?abline to see how to use it.

B. Calculate the annual runoff coefficient (Q/P) for each water year and plot its timeseries.

I’m not showing he code to do this… There’s one value that looks odd…is that value realistic? Why or why not? Look back at Ex 4 and the number of days with rainfall values.

Exclude the year with low rainfall due to missing data. You could remove the whole year, or just set the value to NA in the time series.

Q.P.merge$Rainfall.mm[Q.P.merge$Rainfall.mm<50]=NA
# And with the Q.P time series:
Q.P[is.na(Q.P.merge$Rainfall.mm)] = NA

And replot:

plot(Q.P.merge$WaterYear,Q.P,xlab="Year",ylab="Annual Q/P",type="l")

C. Multipanel plots.

In order to clarify the relationship between P, and Q/P, plot both on a single, multipanel plot.

par(mfrow=c(2,1))  # mfrow creates a multipanel plot.  c(3,1) makes the plot have 3 rows and 1 column.
par(mar=c(0,0,0,0),oma=c(4,4,1,1))  # mar is the margins between each plot.  oma is the outer margins of the plot.
# First panel.
plot(Q.P.merge$WaterYear,Q.P.merge$Rainfall.mm,xaxt="n",type="l")
  #  xaxt="n" means don't plot the x-axis tics or labels.
axis(side=1,labels=FALSE)  #  side=1 means the bottom, side=2 means the left side.
mtext("Annual P, mm",side=2,line=2.5) # adds text to the left side (y-axis).

# Second panel.
plot(Q.P.merge$WaterYear,Q.P,xlab="Year",ylab="Annual Q/P",type="l")
mtext("Annual Q/P",side=2,line=2)

D. Plotting multiple colors for points in a series.

Plot Q/P vs P, and color the points by the year

# Generate a vector where the text is "black"
wy = Q.P.merge$WaterYear
colvec = rep("black",times=length(wy))

# Choose some break years, and assign the colors to years in those intervals.
colvec[(wy>1990) & (wy<=2000)] = "grey"  # 1991-2000 will be colored grey
colvec[(wy>2000)] = "white"              # 2000-2015 will be colored white

plot(Q.P.merge$Rainfall.mm,Q.P,xlab="Precipitation, mm",ylab="Q/P",bg=colvec,col="black",pch=22)
#  "bg" is the argument for the fill color.
#  pch is the symbol type

#  try ?pch to see the options.

legend("bottomright",c("1965-1990","1991-2000","2001-2022"),pch=22,pt.bg=c("black","grey","white"))

  # legend uses the argument "pt.bg" for the background color of the points

E. Plotting multiple colors/symbols: Method using “points”

A second way to plot different series as different point symbols is to use “points”.

#  First, add Q/P to the original dataframe that has all the data.  This makes sure that Q/P is attached to each of the smaller dataframes you'll make in the next step.
Q.P.merge$Q.P = Q.P 
# Next, separate the one dataframe into 3 separate dataframes by the year.
Q.P.merge.pre.1990 = Q.P.merge[Q.P.merge$WaterYear <= 1990,]
Q.P.merge.1990.2000 = Q.P.merge[(Q.P.merge$WaterYear > 1990) & (Q.P.merge$WaterYear <= 2000),]
Q.P.merge.post.2000 = Q.P.merge[Q.P.merge$WaterYear > 2000,]

#  Find the full range of y values in the original dataframe.
y.value.range = range(Q.P.merge$Q.P,na.rm=TRUE)
x.value.range = range(Q.P.merge$Rainfall.mm,na.rm=TRUE)

#  Now, plot the first series using "plot"
plot(Q.P.merge.pre.1990$Rainfall.mm,Q.P.merge.pre.1990$Q.P,col="black",pch=22,xlim=x.value.range,ylim=y.value.range)

# Now, add additional series using "points"
points(Q.P.merge.1990.2000$Rainfall.mm,Q.P.merge.1990.2000$Q.P,col="grey",pch=19)
points(Q.P.merge.post.2000$Rainfall.mm,Q.P.merge.post.2000$Q.P,col="blue",pch=12)

# Add a legend
legend("bottomright",c("1965-1990","1991-2000","2001-2022"),pch=c(22,19,12),col=c("black","grey","blue"))

Plotting tip: let’s turn the y-axis numbers so they are not turned 90 degrees:

plot(Q.P.merge.pre.1990$Rainfall.mm,Q.P.merge.pre.1990$Q.P,col="black",pch=22,xlim=x.value.range,ylim=y.value.range,las=1)

# Now, add additional series using "points"
points(Q.P.merge.1990.2000$Rainfall.mm,Q.P.merge.1990.2000$Q.P,col="grey",pch=19)
points(Q.P.merge.post.2000$Rainfall.mm,Q.P.merge.post.2000$Q.P,col="blue",pch=12)

# Add a legend
legend("bottomright",c("1965-1990","1991-2000","2001-2022"),pch=c(22,19,12),col=c("black","grey","blue"))

For the HW, you’ll make the axis labels nice…!

You’ll also plot rainfall vs runoff:

plot(Q.P.merge$Rainfall.mm,Q.P.merge$Q.ann.mm,xlab="Annual Precipitation, mm",ylab="Annual Q, mm",bg=colvec,col="black",pch=22)
#  "bg" is the argument for the fill color.
#  pch is the symbol type

#  try ?pch to see the options.

legend("bottomright",c("1965-1990","1991-2000","2001-2022"),pch=22,pt.bg=c("black","grey","white"))