#######################################################################################
## Urbanization influence on flow record: Nonstationarity
## 
## Joel Sholtes, PhD Candidate
## Colorado State University
## Dept. of Civil and Env. Engr
## March, 2015
##
## Imports flow records from USGS NWIS and divides flow up into decades
## to study how urbanization influence the FDC and MFA analysis.
## I don't have sediment load data for these sites and therefore
## explore sensitivity to changing means and variance over time
## on MFA as a function of beta, the slope of the sediment rating curve.
#########################################################################################

# Clear data
rm(list = ls())
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis")

# load custom water download function modified from importDVs{waterData}
# can also use the importDVs function from waterData instead.
source("ImportDV.R")
library(XML) # To run ImportDV function
library(Hmisc)
library(RColorBrewer)
library(MASS) # fitdistr

# Urban sites from Tyler Rosberg, MS student 2015
k = 2
sites <- as.character(c(12120500, 12120000, 12127100, 12091200, 12073500, 12069550, 12135000, 12108500))
names <- c("Juanita Ck nr Kirkland, WA", # urban
           "Mercer Ck nr Bellevue, WA", # urban
           "Swamp Ck at Kenmore, WA", # urban
           "Leach Ck nr Fircrest, WA", # urban
           "Huge Ck nr Wauna, WA",  # rural
           "Big Beef Ck nr Seabeck, WA",
           "Wallace River at Gold Bar ,WA",
           "Newakum Ck nr Black Diamond, WA") # rural
short.names <- c("Juanita", "Mercer", "Swamp", "Leach", "Huge", "Beef", "Wallace", "Newakum")

# This loop calls on either the importDV or importIV functions to download and write the flow series
# directly from NWIS. Note that importIV is not currently working. See email from USGS R lady
# it also replaces any negative flow values with NAs and removes leap year days
# finally, it starts and stops the flow record on 10-01 and 09-30, respectively

## **ONLY DO THIS ONCE**  ##
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Uncertainty/Urban")
for ( i in 1:length(sites) ) {
    staid <- sites[i] # Station ID colums (character)
    sdate <- "1900-09-30"
    edate <- "2014-10-01"
    
    flow.file <- importDV(staid, code = "00060", stat = "00003", 
                            sdate, edate) 

    file.name <- staid
    write.table(flow.file, file = paste("UrbanFlow/", file.name, ".csv", sep = ""), 
                row.names = FALSE, sep = ",")
}

###################################################################################
# Now that flows have been downloaded and saved, read in flow file and get to work!
###################################################################################

setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Uncertainty/Urban")
staid <- sites[k]
path <- paste("UrbanFlow/", staid, ".csv", sep = "")
flow.file <- read.csv(path, stringsAsFactors = F)
flow.file$dates <- as.Date(flow.file$dates, format = "%Y-%m-%d")

# Consult range of flow data
flow.file$dates[1]
flow.file$dates[nrow(flow.file)]

# Plot Daily flow data for POR
  par(mfrow = c(1,1), mar = c(4,5,1,1), oma = c(0,0,0,0), pty = "m")
  plot(flow.file$dates, flow.file$val/35.3146667, type = "l", 
       xlab = "", ylab = expression("Discharge ("*m^3/s*")"),
       xaxt="n")
  axis(1,at = as.Date(paste(decades[2:length(decades)], "-01-01", sep = "")), labels = decades[2:length(decades)])
  label = paste(staid, " - ", names[k], sep = "")
  text(x=grconvertX(.05, from='npc'), y = grconvertY(.9, from='npc'), labels = label, cex = 1, adj = c(0,0))

# Manually set start dates based on full or nearly full decades
start.date <- as.Date("1960-01-01",format = "%Y-%m-%d")
end.date <- as.Date("2009-12-31",format = "%Y-%m-%d")

decades <- seq(as.numeric(format(start.date, "%Y")), as.numeric(format(end.date, "%Y")), by = 10)
#decades <- decades[c(1,2,4,5,6)] # remove 1970 from flow record

for (i in 1:length(decades)){
  assign(paste("flow.",decades[i], sep = ""), 
         subset(flow.file, 
                as.numeric(substr(as.character(format(flow.file$dates, "%Y")),1,3)) == 
                as.numeric(substr(as.character(decades[i]),1,3))
                )
        )
}

decades <- c(1900, decades)
flow.1900 <- flow.file

# Calc mean and sd by decade
mean <- numeric(length(decades))
sd <- numeric(length(decades))

for (i in 1:length(decades)){
  txt <- paste("x <- flow.",decades[i],"$val/35.3146667", sep = "")
  eval(parse(text=txt))
  mean[i] <- mean(x, na.rm = T)
  sd[i] <- sd(x, na.rm = T)
}

# Plot Daily flow data for POR
name <- paste("HG_",short.names[k], ".pdf",sep="")
pdf(name, width = 8, height = 3)
  par(mfrow = c(1,1), mar = c(4,5,1,4), oma = c(0,0,0,0), pty = "m")
  begin <- as.Date(paste(decades[2:length(decades)], "-01-01", sep = ""))
  end <- as.Date(paste(decades[2:length(decades)]+9, "-12-31", sep = ""))
  plot(flow.file$dates, flow.file$val/35.3146667, type = "l", 
       xlab = "", ylab = expression("Discharge ("*m^3/s*")"))
  par(new=TRUE)
  plot(flow.file$dates, flow.file$val/35.3146667, type = "n", ylim = c(0,3), axes = F, xlab = "", ylab = "")
  #abline(h = mean[i], lty = 5, col = "gray50", lwd = 2) 
  segments(x0=begin, y0=mean[2:length(decades)], x1=end, y1=mean[2:length(decades)], 
           lty = 6, col = "royalblue", lwd = 2)
  #abline(h = sd[1], lwd = 2, lty = 4, col = "firebrick4")
  segments(x0=begin, y0=sd[2:length(decades)], x1=end, y1=sd[2:length(decades)], 
         lty = 2, col = "firebrick", lwd = 2)
  axis(4, col = "royalblue", col.axis = "royalblue")
  box()
  label = paste(staid, " - ", names[k], sep = "")
  text(x=grconvertX(.05, from='npc'), y = grconvertY(.9, from='npc'), labels = label, cex = 1, adj = c(0,0))
  legend("topright", c("Mean", "SD"), lty = c(6,2), lwd = 2, bty = "n", col = c("royalblue", "firebrick"), cex = 0.8)
dev.off()

###############################################################################################
## Empirical MFA - *** NOTE does not work well with short term (10 year) daily flow records ***
###############################################################################################

## Create empirical PDF of flows...
# number of points to differentiate
npts <- 50 
pdf.dec <- matrix(nrow = npts, ncol = length(decades))
xpts.dec <- matrix(nrow = npts, ncol = length(decades))

for (i in 1:length(decades)){
  eq <- paste("flow <- flow.", decades[i], "$val/35.3146667", sep = "")
  eval(parse(text=eq))
  # create logarithmic sequence of Q's at which to differentiate
  xpts <- exp( seq( log(min(flow) + 0.001), 
                    log(max(flow) - 0.001), length = npts) )
  
  # create empirical CDF using flow record
  ecdf.prob <- ( 1:length(flow) ) / length(flow)
  flow.sort <- sort(flow)
  
  # create sequence of interpolated eCDF values evaluated at 'xpts1'
  ypts <- approx(flow.sort, ecdf.prob, xout = xpts)$y
  
  # add extreme from sorted flow and eCDF to vectors for differentiation
  xpts <- c(min(flow), xpts, max(flow))
  ypts <- c(min(ecdf.prob), ypts, max(ecdf.prob))
  
  # Numerical differentiation of eCDF.
  # 1) Eval eCDF at logarithmically-spaced intervals using central difference (j+1) - (j-1)
  # 2) Assign central diff to PDF at mid point flow value (j)
  pdf.e <- numeric(length(xpts) - 2)
  n <- length(xpts) - 1
  q.e <- xpts[2:n]
  for (j in 2:n) { # j starts at 2, but assigns central diff to PDF at j-1
    pdf.e[j-1] <- ( ypts[j+1] - ypts[j-1] ) / ( xpts[j+1] - xpts[j-1] )
  }
  pdf.dec[,i] <- pdf.e
  xpts.dec[,i] <- q.e
}

# Calc Eff Curve and Sediment Yield Record by Decade as function of Beta
a <- 0.1
b <-  c(1.5, 2, 3, 4)

## Qeff Calc ##
qe.1900 <- matrix(ncol = length(b), nrow = npts)
qe.1950 <- matrix(ncol = length(b), nrow = npts)
qe.1960 <- matrix(ncol = length(b), nrow = npts)
qe.1970 <- matrix(ncol = length(b), nrow = npts)
qe.1980 <- matrix(ncol = length(b), nrow = npts)
qe.1990 <- matrix(ncol = length(b), nrow = npts)
qe.2000 <- matrix(ncol = length(b), nrow = npts)

for (i in 1:length(b)){
  qe.1900[,i] <- a*(xpts.dec[,1])^b[i]*pdf.dec[,1]
  qe.1950[,i] <- a*(xpts.dec[,2])^b[i]*pdf.dec[,2]
  qe.1960[,i] <- a*(xpts.dec[,3])^b[i]*pdf.dec[,3]
  #qe.1970[,i] <- a*(xpts.dec[,3])^b[i]*pdf.dec[,3]
  qe.1980[,i] <- a*(xpts.dec[,4])^b[i]*pdf.dec[,4]
  qe.1990[,i] <- a*(xpts.dec[,5])^b[i]*pdf.dec[,5]
  qe.2000[,i] <- a*(xpts.dec[,6])^b[i]*pdf.dec[,6]
}

## Plot Eff Curves ##
pal <- colorRampPalette(c("gray", "cadetblue", "firebrick"))
cols <- pal(length(decades)) #brewer.pal(length(decades)+2,"YlOrRd")
#cols <- c("gray", cols[2:length(decades)])
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Uncertainty/Urban")

pdf("Eff_Huge.pdf", width = 6, height = 5)
  par(mfrow=c(length(b),1), mar = c(0,5,0,0), oma = c(5,0,1,1), pty = "m")
  
  for(i in 1:length(b)){
    plot(xpts.dec[,1], qe.1960[,i], type = "n", 
         xlim=range(xpts.dec),
         ylim = c(0, max(qe.1960[,i], qe.1960[,i], qe.1980[,i], qe.1990[,i], qe.2000[,i])),
         xaxt = ifelse(i < length(b), "n", "s"),
         yaxt = "n",
         las = 1,
         ylab = bquote(beta == .(b[i])),
         cex.lab = 1.2,
         cex.axis = 1.2
         )
  
    if( i == 1){
      legend("topright",
             legend = c("POR",decades[2:length(decades)]),
             col = cols,
             lty = 1:length(decades),
             bty = "n",
             cex = 1,
             y.intersp = 0.8,
             lwd = c(2.5, rep(1.5, length(decades)-1)))
      label = paste(staid, " - ", names[k], sep = "")
      text(x=grconvertX(.5, from='npc'), y = grconvertY(.9, from='npc'), labels = label, cex = 1)
    }
    
    for(j in 1:length(decades)){
      txt1 <- paste("lines(xpts.dec[,j], qe.",decades[j],"[,i], col = cols[j],
                    lwd = ifelse(j==1,2.5,1.5), lty = j)", sep = "")
      eval(parse(text=txt1))
    }
    for(j in 1:length(decades)){
    txt2 <- paste("points(xpts.dec[which.max(qe.", decades[j],"[,i]),j], max(qe.", decades[j],"[,i]),",
                  "pch = 20, cex = 2, col = cols[j])", sep = "")
    eval(parse(text=txt2))
    }
  }
  mtext(expression("Discharge ("*m^3/s*")"), side = 1, line = 3.5)
  mtext("Sediment Yield Density", side = 2, line = -4, outer = T)

dev.off()

### Qs50 Calc ###
nrow = c(nrow(flow.1900),nrow(flow.1960),nrow(flow.1970),nrow(flow.1980), nrow(flow.1990), nrow(flow.2000))

cum.1900 <- matrix(ncol = length(b), nrow = nrow[1])
cum.1960 <- matrix(ncol = length(b), nrow = nrow[2])
cum.1970 <- matrix(ncol = length(b), nrow = nrow[3])
cum.1980 <- matrix(ncol = length(b), nrow = nrow[4])
cum.1990 <- matrix(ncol = length(b), nrow = nrow[5])
cum.2000 <- matrix(ncol = length(b), nrow = nrow[6])

for (i in 1:length(b)){
  for (j in 1:length(decades)) {
  eq <- paste("cum.", decades[j], "[,i] <- cumsum(a*(sort(flow.", decades[j], "$val/35.3146667))^b[i])/",
              "sum(a*(sort(flow.", decades[j], "$val/35.3146667))^b[i])", sep = "")
  eval(parse(text=eq))
  }
}

## Plot Cum Yield Curves ##
pal <- colorRampPalette(c("gray", "purple4"))
cols <- pal(length(decades))

setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Uncertainty/Urban")

pdf("Cum_Leach.pdf", width = 10, height = 3.5)
  par(mfrow=c(1,length(b)), mar = c(5,0,2,0), oma = c(0,5,1,1), pty = "s")
  
  for(i in 1:length(b)){
    plot(sort(flow.1900$val), cum.1900[,i], type = "n", 
         xlim=range(flow.file$val),
         ylim = c(0,1),
         yaxt = ifelse(i == 1, "s","n"),
         xlab = "",
         las = 1,
         cex.lab = 1.2,
         cex.axis = 1.2,
         main = bquote(beta == .(b[i]))
    )
    abline(h = 0.5, lty = 2, col = "gray")
    
    if( i == 1){
      legend("bottomright",
             legend = c("POR",decades[2:length(decades)]),
             col = cols,
             lty = 1:6,
             bty = "n",
             cex = 1.2,
             y.intersp = 0.8,
             lwd = c(2.5, rep(1.5, length(decades)-1)))
    }
    
    for(j in 1:length(decades)){
      txt1 <- paste("lines(sort(flow.", decades[j],"$val), cum.",decades[j],"[,i], col = cols[j],
                    lwd = ifelse(j==1,2.5,1.5), lty = j)", sep = "")
      eval(parse(text=txt1))
    }
  }
  mtext(expression("Discharge ("*m^3/s*")"), side = 1, line = -1.5, outer = T)
  mtext("Cum. Sediment Yield", side = 2, line = 3, outer = T)
  mtext(paste(staid, " - ", names[k], sep = ""), side = 3, line = -.5, outer = T, cex = 0.8)
  
dev.off()

####################################
### Calc and plot FDCs by decade ###
####################################
p = seq(.0001,.9999, length.out = 500)

fdc.dec <- matrix(ncol = length(decades), nrow = length(p))

n <- length(decades)
rhs <- paste("hdquantile(flow.", decades[1:n], "$val/35.3146667, probs = p)", sep = "")
lhs <- paste("fdc.dec[,", 1:n, "]", sep = "")
eq <- paste(paste(lhs, rhs, sep="<-"), collapse=";")
eval(parse(text=eq))


setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Uncertainty/Urban")

### Plot FDCs ###
name <- paste("FDC_",short.names[k], ".pdf",sep="")
pdf(name, width = 8, height = 5)
  par(mar = c(5,5,1,1), mfrow = c(1,1), oma = c(0,0,0,0), pty = "m")
  x.norm <-  qnorm(1 - p)
  plot(x.norm,fdc.dec[,1], type = "l", 
       ylim = range(fdc.dec[,1:length(decades)]), 
       log = "y",
       xaxt = "n",
       xlab = "Probability of Exceedence",
       ylab = expression("Discharge ("*m^3/s*")"),
       lwd = 2,
       cex.lab = 1.2,
       col = "grey50")
  for(j in 2:length(decades)){
    lines(x.norm, fdc.dec[,j], lty = j, lwd = 1.5)
  }

  xtic <- c(qnorm(c(0.1, 1, 2, 5, 10, 25, 50, 75, 90, 95, 99, 99.9)/100))
  xlab <- c(0.1, 1, 2, 5, 10, 25, 50, 75, 90, 95, 99, 99.9) 
  axis(1, at = xtic, labels = xlab) 
  label = paste(staid, " - ", names[k], sep = "")
  text(x=grconvertX(.7, from='npc'), y = grconvertY(.95, from='npc'), labels = label, cex = 0.9)
  legend("bottomleft",
         c("POR", decades[2:length(decades)]),
         lty = 1:6,
         bty = "n",
         lwd = c(2,1.5,1.5,1.5, 1.5, 1.5, 1.5, 1.5),
         col = c("gray40","black","black","black", "black", "black", "black"),
         cex = 0.8
         )

dev.off()

#################################################################################################
## Theoretical Fitting to decadal flow records..
## They are all too short! This means that the flow in the tails override the effectiveness calc
#################################################################################################

# Fit Lognormal Distribution using MLE estimators for parameters
# NOTE: MLE estimators do not reflect trends in mean or variance well
# therefore, I did not use MLE
mean.log <- numeric(length(decades))
sd.log <- numeric(length(decades))

# for (i in 1:length(decades)){
#   txt <- paste("x <- flow.",decades[i],"$val/35.3146667", sep = "")
#   eval(parse(text=txt))
#   fit<-fitdistr(x,"lognormal")
#   mean.log[i] <- as.numeric(fit$estimate[1])
#   sd.log[i] <- as.numeric(fit$estimate[2])
# }

# Method of Moments estimators <- This reflects trends better
for (i in 1:length(decades)){
  txt <- paste("x <- flow.",decades[i],"$val/35.3146667", sep = "")
  eval(parse(text=txt))
  mu <- mean(x)
  sig <- sd(x)
  cv <- sig/mu
  mean.log[i] <- 1/2*log(mu^2/(1+cv^2))
  sd.log[i] <- log(1+cv^2)
}

flow <- seq(min(flow.file$val/35.3146667), max(flow.file$val/35.3146667*4), length.out = 1000)

# Create Lognormal Effectiveness curves
eff.dec <- array(0,dim = c(length(decades), length(flow), length(b)))

for (i in 1:length(decades)){
  for (j in 1:length(b)){
    eff.dec[i,,j] <- a*flow^b[j]*dlnorm(flow, meanlog = mean.log[i], sdlog = sd.log[i])
  }
}

## Create Lognormal Cumulative Sed Yield Curves
cum.dec <- array(0,dim = c(length(decades), length(flow), length(b)))

for (i in 1:length(decades)){
  for (j in 1:length(b)){
    cum.dec[i,,j] <- pnorm((log(flow) - mean.log[i] - b[j]*sd.log[i]^2) / sd.log[i])
  }
}

## Plot Eff Curves ##
pal <- colorRampPalette(c("gray", "cadetblue", "firebrick"))
cols <- pal(length(decades)) #brewer.pal(length(decades)+2,"YlOrRd")
#cols <- c("gray", cols[2:length(decades)])
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Uncertainty/Urban")
name <- paste("Eff_",short.names[k], "L.pdf",sep="")
pdf(name, width = 6, height = 5)
par(mfrow=c(length(b),1), mar = c(0,5,0,0), oma = c(5,0,1,1), pty = "m")

for(i in 1:length(b)){
  plot(flow, eff.dec[1,,i], type = "n", 
       xlim=range(flow/2),
       ylim = c(0, max(eff.dec[,,i])),
       xaxt = ifelse(i < length(b), "n", "s"),
       yaxt = "n",
       las = 1,
       ylab = bquote(beta == .(b[i])),
       cex.lab = 1.2,
       cex.axis = 1.2
  )
  
  if( i == 1){
    legend("topright",
           legend = c("POR",decades[2:length(decades)]),
           col = cols,
           lty = 1:length(decades),
           bty = "n",
           cex = 1,
           y.intersp = 0.8,
           lwd = c(2.5, rep(1.5, length(decades)-1)))
    label = paste(staid, " - ", names[k], sep = "")
    text(x=grconvertX(.5, from='npc'), y = grconvertY(.9, from='npc'), labels = label, cex = 1)
  }
    for(j in 1:length(decades)){
      lines(flow, eff.dec[j,,i], col = cols[j], lwd = ifelse(j==1,2.5,1.5), lty = j)
    }
    for(j in 1:length(decades)){
      points(flow[which.max(eff.dec[j,,i])], max(eff.dec[j,,i]), pch = 20, cex = 2, col = cols[j])
    }
  }
  mtext(expression("Discharge ("*m^3/s*")"), side = 1, line = 3.5)
  mtext("Sediment Yield Density", side = 2, line = -4, outer = T)

dev.off()

## Plot Cum Curves ##
pal <- colorRampPalette(c("gray", "purple4"))
cols <- pal(length(decades))

setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Uncertainty/Urban")
name <- paste("Cum_",short.names[k], "L.pdf",sep="")
pdf(name, width = 10, height = 3.5)
par(mfrow=c(1,length(b)), mar = c(5,0,2,0), oma = c(0,5,1,1), pty = "s")

for(i in 1:length(b)){
  plot(flow, cum.dec[1,,i], type = "n", 
       xlim=range(flow),
       ylim = c(0,1),
       yaxt = ifelse(i == 1, "s","n"),
       xlab = "",
       las = 1,
       cex.lab = 1.2,
       cex.axis = 1.2,
       main = bquote(beta == .(b[i]))
  )
  abline(h = 0.5, lty = 2, col = "gray")
  
  if( i == 1){
    legend("bottomright",
           legend = c("POR",decades[2:length(decades)]),
           col = cols,
           lty = 1:6,
           bty = "n",
           cex = 1.2,
           y.intersp = 0.8,
           lwd = c(2.5, rep(1.5, length(decades)-1)))
  }
  
  for(j in 1:length(decades)){
    lines(flow, cum.dec[j,,i], col = cols[j], lwd = ifelse(j==1,2.5,1.5), lty = j)
  }
}
  mtext(expression("Discharge ("*m^3/s*")"), side = 1, line = -1.5, outer = T)
  mtext("Cum. Sediment Yield", side = 2, line = 3, outer = T)
  mtext(paste(staid, " - ", names[k], sep = ""), side = 3, line = -.5, outer = T, cex = 0.8)

dev.off()


###################################################################
## Theoretical using Mercer as info and smooth increase in sd ##
###################################################################

# NOTE: Ended up not using this b/c MoM parameter estimation shows trends well

# Assuming that only flow variability increases as with Mercer
mu.x <- 0.6339469 #mean[1], POR
sd.x <- 0.7740903 # sd[1], POR
min.sd.x <- 0.602178 #sd[2] (lowest val)
meanlog <- -0.8271773 #mean.log[1] from MLE fit
sdlog <- 0.7966851 # sd.log[2]

# MLE parameter fits from Mercer Ck data for lognormal 1960:2000 by decade
sdlog.trend <- c(0.7966851, 0.7351751, 0.7016280, 0.7961895, 0.8065082)
sd.trend <- c(0.6021780, 0.7435384, 0.7705797, 0.9265376, 0.8562760)
time <- 1:5
sd.lm <- lm(sd.trend ~ time)
sd.coef <- coef(sd.lm)
sd.x <- c(sd.x, min.sd.x, sd.coef[1] + sd.coef[2]*2:5)
cv <- sd.x/mu.x

mu.y <- 1/2*log(mu.x^2/(1+cv^2))
sigma.y <- log(1+cv^2)

# Create Lognormal Effectiveness curves
eff.dec <- array(0,dim = c(length(decades), length(flow), length(b)))

n <- length(decades)

for (i in 1:n){
  for (j in 1:length(b)){
    eff.dec[i,,j] <- a*flow^b[j]*dlnorm(flow, meanlog = mu.y[i], sdlog = sigma.y[i])
  }
}

## Create Lognormal Cumulative Sed Yield Curves
cum.dec <- array(0,dim = c(length(decades), length(flow), length(b)))

for (i in 1:n){
  for (j in 1:length(b)){
    cum.dec[i,,j] <- pnorm((log(flow) - mu.y[i] - b[j]*sigma.y[i]^2) / sigma.y[i])
  }
}

## Plot Eff Curves ##
pal <- colorRampPalette(c("gray", "cadetblue", "firebrick"))
cols <- pal(length(decades)) #brewer.pal(length(decades)+2,"YlOrRd")
#cols <- c("gray", cols[2:length(decades)])
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Uncertainty/Urban")

pdf("Eff_MercerLT.pdf", width = 6, height = 5)
par(mfrow=c(length(b),1), mar = c(0,5,0,0), oma = c(5,0,1,1), pty = "m")

for(i in 1:length(b)){
  plot(flow, eff.dec[1,,i], type = "n", 
       xlim=range(flow/2),
       ylim = c(0, max(eff.dec[,,i])),
       xaxt = ifelse(i < length(b), "n", "s"),
       yaxt = "n",
       las = 1,
       ylab = bquote(beta == .(b[i])),
       cex.lab = 1.2,
       cex.axis = 1.2
  )
  
  if( i == 1){
    legend("topright",
           legend = c("POR",decades[2:length(decades)]),
           col = cols,
           lty = 1:length(decades),
           bty = "n",
           cex = 1,
           y.intersp = 0.8,
           lwd = c(2.5, rep(1.5, length(decades)-1)))
    label = paste(staid, " - ", names[k], sep = "")
    text(x=grconvertX(.5, from='npc'), y = grconvertY(.9, from='npc'), labels = label, cex = 1)
  }
  
  for(j in 1:length(decades)){
    lines(flow, eff.dec[j,,i], col = cols[j], lwd = ifelse(j==1,2.5,1.5), lty = j)
  }
  for(j in 1:length(decades)){
    points(flow[which.max(eff.dec[j,,i])], max(eff.dec[j,,i]), pch = 20, cex = 2, col = cols[j])
  }
}
mtext(expression("Discharge ("*m^3/s*")"), side = 1, line = 3.5)
mtext("Sediment Yield Density", side = 2, line = -4, outer = T)

dev.off()

## Plot Cum Curves ##
pal <- colorRampPalette(c("gray", "purple4"))
cols <- pal(length(decades))

setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Uncertainty/Urban")

pdf("Cum_MercerTL.pdf", width = 10, height = 3.5)
par(mfrow=c(1,length(b)), mar = c(5,0,2,0), oma = c(0,5,1,1), pty = "s")

for(i in 1:length(b)){
  plot(flow, cum.dec[1,,i], type = "n", 
       xlim=range(flow),
       ylim = c(0,1),
       yaxt = ifelse(i == 1, "s","n"),
       xlab = "",
       las = 1,
       cex.lab = 1.2,
       cex.axis = 1.2,
       main = bquote(beta == .(b[i]))
  )
  abline(h = 0.5, lty = 2, col = "gray")
  
  if( i == 1){
    legend("bottomright",
           legend = c("POR",decades[2:length(decades)]),
           col = cols,
           lty = 1:6,
           bty = "n",
           cex = 1.2,
           y.intersp = 0.8,
           lwd = c(2.5, rep(1.5, length(decades)-1)))
  }
  
  for(j in 1:length(decades)){
    lines(flow, cum.dec[j,,i], col = cols[j], lwd = ifelse(j==1,2.5,1.5), lty = j)
  }
}
mtext(expression("Discharge ("*m^3/s*")"), side = 1, line = -1.5, outer = T)
mtext("Cum. Sediment Yield", side = 2, line = 3, outer = T)
mtext(paste(staid, " - ", names[k], sep = ""), side = 3, line = -.5, outer = T, cex = 0.8)

dev.off()
