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")

A. Save your extracted DAYMET time series.

First, run the DAYMET extraction for your watershed (Exercise 3), and then save the DAYMET daily precipitation series in your personal class directory, so you don’t have to re-run the DAYMET loop every time:

outdir = << your class directory >>
write.table(prcp.all,file=paste0(outdir,"DAYMET_mywatershed.txt"))  # where mywatershed is the name of your watershed
  #  write.table is used for data frames

B. Read back in your DAYMET series.

Close R (don’t save the workspace), and restart R. Now, reload your DAYMET precipitation data for the lower SDR:

  # Change "outdir", which is where you wrote your DAYMET rainfall file.
outdir = "G:/mydocuments/SDSU_classes/past/2016_2017_100_576_574/2016F_100_576/GEOG576/hw/data/"
  # Change the filename "DAYMET_lowerSDR.txt" to your DAYMET file name.
daymet = read.table(paste0(outdir,"DAYMET_lowerSDR.txt"))
#  and look at the data:
head(daymet)
##         Date Precip.mm
## 1 1980-01-01         0
## 2 1980-01-02         0
## 3 1980-01-03         0
## 4 1980-01-04         0
## 5 1980-01-05         0
## 6 1980-01-06         0

What class is your Date column?

class(daymet$Date)

Date was read as a “factor”, which messes up subsequent analysis, so reformat Date as a Date class:

daymet$Date = as.Date(daymet$Date,format="%Y-%m-%d")

And confirm that it worked as expected:

head(daymet)
class(daymet$Date)

# Also see what class the Precip.mm column is in--will that class work?
class(daymet$Precip.mm)

C. Summarize DAYMET data to annual values by water year and plot.

1. Generate the WaterYear (as in Exercise 2, IV.C).
a. Here is a simplified version of Ex 2 IVC. Note the use of “from” and “to” in the seq function.
breaks <- seq(from=as.Date("1980-10-01"), to=as.Date("2015-10-01"), by="year")  
years.breaks = as.numeric(format(breaks,"%Y"))
labels.wy = years.breaks[2:length(breaks)]  # Take all but the first year in the breaks.  Why?
daymet$WaterYear <- cut(daymet$Date, breaks,labels=labels.wy)

Look at your data. Is WaterYear correct?

# First look at the first 600 lines
daymet[1:600,]

# Now look at the last 600 lines.  Is WaterYear correct?
nrecords = length(daymet$Date)
daymet[(nrecords-600):nrecords,]
  1. Aggregate the data to annual total precipitation (last chunk of code in Ex 2 IVC), rename the columns of the data.frame, and look at the first 6 lines.
##   WaterYear Precip.mm
## 1      1981  285.5686
## 2      1982  458.7020
## 3      1983  733.0353
## 4      1984  197.8824
## 5      1985  336.6706
## 6      1986  517.6235

Look at the class of the y

class(daymet.annual.wy$WaterYear)
## [1] "factor"

It’s a factor, which is a difficult class to use. For example, it results in crummy plots:

plot(daymet.annual.wy$WaterYear,daymet.annual.wy$Precip.mm,ylab="Annual Precipitation, mm",xlab="",type="l")

What if you try to convert this factor to numeric?

as.numeric(daymet.annual.wy$WaterYear)
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
## [24] 24 25 26 27 28 29 30 31 32 33 34 35

So unfortunately, we need to do “as.character” first to change the factor into a character, then use as.numeric to convert the character to a number. You can do this in one nested command:

as.numeric(as.character(daymet.annual.wy$WaterYear))
##  [1] 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994
## [15] 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008
## [29] 2009 2010 2011 2012 2013 2014 2015
daymet.annual.wy$WaterYear = as.numeric(as.character(daymet.annual.wy$WaterYear))

And see the improvement in the plot.

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 = "11022480"  #  The USGS streamgage code for San Diego River at Mast Blvd
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 = ""
  #  Use your site code in the line below:
x = readNWISdv("11022480",parameter.code,start.date,end.date)
head(x)
##   agency_cd  site_no       Date X_00060_00003 X_00060_00003_cd
## 1      USGS 11022480 1912-05-01            33                A
## 2      USGS 11022480 1912-05-02            33                A
## 3      USGS 11022480 1912-05-03            30                A
## 4      USGS 11022480 1912-05-04            25                A
## 5      USGS 11022480 1912-05-05            23                A
## 6      USGS 11022480 1912-05-06            23                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 11022480 1912-05-01     33  A
## 2      USGS 11022480 1912-05-02     33  A
## 3      USGS 11022480 1912-05-03     30  A
## 4      USGS 11022480 1912-05-04     25  A
## 5      USGS 11022480 1912-05-05     23  A
## 6      USGS 11022480 1912-05-06     23  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 for the DAYMET data, but with mean.

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("1969-10-01"), to=as.Date("2015-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("Year","MeanQcfs")
annual.mean.cfs$Year = as.numeric(as.character(annual.mean.cfs$Year))
annual.mean.cfs$MeanQcfs = as.numeric(as.character(annual.mean.cfs$MeanQcfs))
head(annual.mean.cfs)
##   Year  MeanQcfs
## 1 1970  3.663178
## 2 1971  5.548438
## 3 1972  5.624809
## 4 1973 13.597260
## 5 1974  9.771233
## 6 1975 18.263014
2. Calculate annual runoff in mm/year.
  1. 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] 368
  1. Now calculate annual discharge from your watershed in mm (write out the units conversion first on paper!).

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,daymet.annual.wy,by.x="WaterYear",by.y="WaterYear")

#  Print Q.P.merge to make sure it worked.
plot(Q.P.merge$Precip.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) and plot its timeseries.

Q.P = Q.P.merge$Q.ann.mm/Q.P.merge$Precip.mm
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$Precip.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$Precip.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("1980-1990","1991-2000","2001-2015"),pch=22,pt.bg=c("black","grey","white")) 

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