0. Functions in R

Functions are scripts that you want to run over and over again, but not have to type the details every time. For example, if you want to calculate the area of a circle (A = pi*r2) but don’t want to have to write the equation, you could make a function:

circle.area <- function(r){   #  Input value is "r", the radius.
  A = pi*r^2
  return(A)
} 

A.mycircle = circle.area(r=10)  # Try for a circle of radius 10 (arbitrary units)
print(A.mycircle)
## [1] 314.1593

I. C-Q plots for hysteresis

A. Define arrowplot function

You can write your own functions in R. Here’s one that plots arrows between consecutive points. Copy, paste and execute this code to define the function (nothing happens…it just puts the function in memory):

arrowLine <- function(x, y, N, ...){
  lengths <- c(0, sqrt(diff(x)^2 + diff(y)^2))
  l <- cumsum(lengths)
  tl <- l[length(l)]
  el <- seq(0, to=tl, length=N+1)[-1]
  plot(x, y, t="l", ...)
  for(ii in el){
    int <- findInterval(ii, l)
    xx <- x[int:(int+1)]
    yy <- y[int:(int+1)]
    ## points(xx,yy, col="grey", cex=0.5)
    dx <- diff(xx)
    dy <- diff(yy)
    new.length <- ii - l[int]
    segment.length <- lengths[int+1]
    ratio <- new.length / segment.length
    xend <- x[int] + ratio * dx
    yend <- y[int] + ratio * dy
    #points(xend,yend, col="black", pch=19)
    arrows(x[int], y[int], xend, yend, length=0.15)
  }
  points(x,y,pch=20)
}

B. Now use it to plot Q-C for a fake storm:

par(mfrow=c(1,1)) # get back to a single-panel plot

# Make up a fake storm with Q and C
Q = c(1,3,10,15,12,8,3,2,1)
C = c(1,10,50,40,30,10,5,4,1)

# Plot it.
arrowLine(Q,C,N=10,xlab="Q",ylab="C")

C. Now do for “real” data at Los Coches Creek:

1. Load the water quality data

Url of the CSV file is “https://docs.google.com/spreadsheets/d/e/2PACX-1vRSE-bOj9P0AbsZNblF4rg18p9trJK7MmtlFXg8Bqt62eNeX22NTVTTw1B6UkDJZpdaCi3b_J8Xp8ib/pub?output=csv”

Load the data into a data.frame called “xwq”. I’m not showing the code for this…

2. Plot points, no arrows
plot(xwq$Q.cfs,xwq$TN.mgL,xlab="Q cfs",ylab="TN mg/L")

3. Now plot with arrows using the arrowline function from Part A. I’m not showing the code.

4. Calculate the event-mean concentration (EMC):
EMC.TN = sum(xwq$Q.cfs * xwq$TN.mgL)/sum(xwq$Q.cfs)
print(EMC.TN)
## [1] 3.03619
5. Write a function to calculate the EMC. Call it “EMC.calc”.

I’m not showing the code for this.

6. Now, check to make sure your function works by comparing the results of the function with the results from Step 4.
EMC.TN.using.function = EMC.calc(Q=xwq$Q.cfs,C=xwq$TN.mgL)
print(EMC.TN.using.function)
## [1] 3.03619