#########################################################################################
## 1: import list of sites to download field discharge meas. data from NWIS
## 2: Plot various channel geometry and discharge relationships to visually ID breaks
##    in relationships indicating changes in channel geometry, hence bankfull
## 3: Manually record this relationship in a table then run FindQbf.R script
##    to determine other bankfull channel geometry values using linear interpolation
##########################################################################################

# Clear data
rm(list = ls())
# set working dir
library(Hmisc) # for minor tick marks
source("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/minor.ticks.axis.R")

setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/Bedload/NWIS XS Analysis")
# import table of USGS site numbers
site.table <- read.csv("site_summary.csv", header = T, stringsAsFactors = F)

sites <- site.table[c("site.no", "site.nm")]

# append bankfull variables as columns to site table
sites[c("n", "Qbf.m", "Wbf.m", "Dbf.m")] <- NA

pdf( "XS_Plots_Qbf.pdf", width = 11, height = 8.5, paper = "special") # turn on PDF creator
# set plot to be square and 4x4 plots
par(mfrow = c(2,2), pty = "s", oma = c(0,0,0,0), 
    cex = 1, cex.lab = 0.8, cex.axis = 0.8, cex.main = 0.8)

for( i in 1:nrow(sites)) {
  site.no <- as.numeric(sites$site.no[i])
  # attach leading zero to site.no if necessary
  if(!is.na(site.no)) {

    if(nchar(site.no) < 8){
      site.no <- paste("0", site.no, sep = "")
    } else {site.no <- as.character(site.no)}
    
    path <- paste("http://waterdata.usgs.gov/nwis/measurements?site_no=",
                  site.no, "&agency_cd=USGS&format=rdb_expanded", sep = "")
    # Note that comment.char = "#" tells R to skip lines beginning with "#" as in
    # USGS RDB files.
    xs <- read.table(path, comment.char="#", header = T, sep = "\t", 
                      stringsAsFactors = F, na.strings = "", fill = T)
    
    if(!is.null(dim(xs)) & ncol(xs) > 2) { 
    
      xs <- xs[-1,] # removes first row, which does not contain data
      
      # Extract numeric values and asign to variables for plotting
      h <- as.numeric(xs$gage_height_va)
      Q <- as.numeric(xs$discharge_va)
      w <- as.numeric(xs$chan_width)
      A <- as.numeric(xs$chan_area)
      wd <- w[which(w != 0 | A != 0)]^2 / A[which(w != 0 | A != 0)]
      wd <- wd[which(!is.infinite(wd))]
      Q1 <- Q[which(w != 0 | A != 0)]
      Q1 <- Q1[which(!is.infinite(wd))]
      d <- A / w
      
      w <- w[which(Q > 0)]
      A <- A[which(Q > 0)]
      wd <- wd[which(Q > 0)]
      h <- h[which(Q > 0)]
      Q <- Q[which(Q > 0)]
      Q1 <- Q1[which(Q > 0)]
      d <- d[which(Q > 0)]
      
#       hg.table <- data.frame(w = w, A = A, h = h, Q = Q, Q1 = Q1, d = d)
      
      # sort data to find bounds for plotting ONLY, NOT interpolation nor analysis
      # need to maintain Q, A, W measurements
      sort.Q <- sort(Q[which(Q != 0)])
      sort.Q1 <- sort(Q1[which(Q1 != 0)])
      sort.wd <- sort(wd[which(wd != 0)])
      sort.A <- sort(A[which(A != 0)])
      sort.h <- sort(h[which(h != 0)])
      sort.w <- sort(w[which(w != 0)])
      
      # Calc Wbf and Dbf based on Qbf estimate and linear interpolation
      Qbf <- site.table$Qbf.est.cfs[i]
      
      if(!is.na(Qbf)){
        xs[,c(9,10,25,26)] <- as.numeric(unlist(xs[,c(9,10,25,26)]))
        xs <- xs[order(xs$discharge_va),]
        n <- sum(!is.na(xs$chan_width))
        
        if( n > 3){
          wbf <- round(approx(xs$discharge_va, xs$chan_width, xout = Qbf)$y, 0)
          dbf <- round(approx(xs$discharge_va, xs$chan_area, xout = Qbf)$y / wbf, 0)
          Qbf.m <- round(Qbf * 0.028316847, 2)
          wbf.m <- round(wbf * 0.3048, 2)
          dbf.m <- round(dbf * 0.3048, 2)
        } else {
          Qbf.m <- NA
          wbf.m <- NA
          dbf.m <- NA
        }
      } else {
        Qbf.m <- NA
        wbf.m <- NA
        dbf.m <- NA
      }
    
    if(!is.na(Qbf)){
    ## Write results to table
    sites[i,3:6] <- c(n, Qbf.m, wbf.m, dbf.m)
    } else {
      sites[i,3:6] <- NA
    }
    ## Plot Q vs. W/D
    plot(log10(Q1), log10(wd), xlab = "Q [cfs]", ylab = "Width:Depth", cex = .5,
         cex.main = 0.7,
         axes = F)
    
    x.min.idx <- min(which(sort.Q1 != 0 & sort.Q1 > 0)) # find first nonzero value
    x.min <- floor(log10(sort.Q1[x.min.idx]))
    x.max <- ceiling(log10(max(sort.Q1)))
    
    y.min.idx <- min(which(sort.wd != 0)) # find first nonzero value
    y.min <- floor(log10(sort.wd[y.min.idx]))
    y.max <- ceiling(log10(max(sort.wd)))
    
    minor.ticks.axis(1,9,mn = x.min, mx = x.max, cex = 0.8)
    minor.ticks.axis(2,9,mn = y.min, mx = y.max, cex = 0.8)
    abline( v = log10(Qbf), lty = 2 )
    box()
    
    ## Plot Q vs stage
    plot(log10(Q), log10(h), xlab = "Q [cfs]", ylab = "stage [ft]", cex = .5,
         axes = F)
    
    x.min.idx <- min(which(sort.Q != 0 & sort.Q > 0)) # find first nonzero value
    x.min <- floor(log10(sort.Q[x.min.idx]))
    x.max <- ceiling(log10(max(sort.Q)))
    
    y.min.idx <- min(which(sort.A != 0)) # find first nonzero value
    y.min <- floor(log10(sort.A[y.min.idx]))
    y.max <- ceiling(log10(max(sort.A)))
    
    minor.ticks.axis(1,9,mn = x.min, mx = x.max, cex = 0.8)
    minor.ticks.axis(2,9,mn = y.min, mx = y.max, cex = 0.8)
    abline( v = log10(Qbf), lty = 2 )
    box()
    
    ## Plot width vs. Q
    plot(log10(Q), log10(w), xlab = "Q [cfs]", ylab = "width [ft]", cex = .5,
         axes = F)
    
    x.min.idx <- min(which(sort.Q != 0)) # find first nonzero value
    x.min <- floor(log10(sort.Q[x.min.idx]))
    x.max <- ceiling(log10(max(sort.Q)))
    
    y.min.idx <- min(which(sort.w != 0)) # find first nonzero value
    y.min <- floor(log10(sort.w[y.min.idx]))
    y.max <- ceiling(log10(max(sort.w)))
    
    minor.ticks.axis(1,9,mn = x.min, mx = x.max, cex = 0.8)
    minor.ticks.axis(2,9,mn = y.min, mx = y.max, cex = 0.8)
    abline( v = log10(Qbf), lty = 2 )
    box()
    
    ## Plot width vs. area
    plot(log10(w), log10(A), xlab = "Width [ft]", ylab = "Area [ft^2]", cex = .5,
         axes = F)
    
    x.min.idx <- min(which(sort.w != 0)) # find first nonzero value
    x.min <- floor(log10(sort.w[x.min.idx]))
    x.max <- ceiling(log10(max(sort.w)))
    
    y.min.idx <- min(which(sort.A != 0)) # find first nonzero value
    y.min <- floor(log10(sort.A[y.min.idx]))
    y.max <- ceiling(log10(max(sort.A)))
    
    minor.ticks.axis(1,9,mn = x.min, mx = x.max, cex = 0.8)
    minor.ticks.axis(2,9,mn = y.min, mx = y.max, cex = 0.8)
    box()
      
    mtext(paste(site.no, "\n", sites$site.nm[i]), side = 3, line = -4, outer = T)
  }
  
  # create data frame for XS data and write table to "XS Analysis" folder
  chan.table <- data.frame(discharge = Q, gage.height = h, width = w, area = A, 
                          depth = d)
  write.csv(chan.table, paste("chan_meas/" , site.no, ".csv", sep = ""), 
           row.names = F) # Need to make this folder first
  
 } 
  
}
dev.off()

write.csv(sites, "results_bankfull.csv", row.names = F)

## Compare EStimated Vs Measured Values

Qbf.meas <- c(425, 197.2, 313, 362.1, 281.9, 30, 220.1, 135.4, 75.4, 
              690.5, 1351.5, 3404, 1306.6, 463.6, 2, 651.3, 24.1, 453.1, 446)
Qbf.est <- c(22000, 6000, 14000, 12000, 18000, 1100, 9000, 11000, 3500, 
             22000, 120000, 13000, 22000, 20000, 80, 4000, 500, 3000, 7000)*0.0283

p.diff <- (Qbf.est - Qbf.meas)/Qbf.meas
sum(p.diff)
mean.pdiff<- round(mean(p.diff),2)
mape <- round(sum(abs(p.diff))/length(p.diff)*100,0)

setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/Figures")
pdf("Qbf_Est_Meas.pdf", height = 6, width = 4)
  par(mfrow = c(2,1), pty = "s", cex.axis = 0.75,
      mar = c(0,4,0,0), oma = c(4,0.5,0.5,0.5))
  plot(log10(Qbf.meas), log10(Qbf.est), 
       xlim = c(0,log10(4000)), 
       ylim = c(0,log10(4000)),
       axes = F,
       xlab = NA, 
       ylab = expression(italic(Q[BF])~"Estimated"~~"("*m^3/s*")"))
  abline(0,1)
  minor.ticks.axis(2,9,mn = 0, mx = 4, cex = 0.8)
  box()
  
  
  plot(log10(Qbf.meas), p.diff,
       xlim = c(0,log10(4000)), 
       ylab = "Percent Error",
       axes = F)
  minor.ticks.axis(1,9,mn = 0, mx = 4, cex = 0.8)
  axis(side = 2, at = c(-0.5, 0, 0.5, 1, 1.5), labels = c(-0.5, 0, 0.5, 1, NA))
  box()
  abline(h = 0)
  abline(h = mean.pdiff, lty = 2)
  text(x = 0.1, y = max(p.diff)*0.9, 
       bquote("MPE = "~.(mean.pdiff*100)~"%"), 
       cex = 0.8,
       adj = c(0,0))
text(x = 0.1, y = max(p.diff)*0.75, 
     bquote("MAPE = "~.(mape)~"%"), 
     cex = 0.8,
     adj = c(0,0))
  mtext(expression(italic(Q[BF])~"Measured"~~"("*m^3/s*")"), side = 1, line = 3)
dev.off()
