#######################################################################################
# importXS.R <- A script to import USGS "Field Measurements" taken when discharge is
# measured at a gage site: top width, flow area, stage
#
# 1) Creates a series of figures showing scatter plot relationships between flow and
# various estimates of channel geometry from these data.
# 2) Following determination of bankfull discharge from at-a-station hydraulic 
# geometry relations(Wiliams 1978), one must manually fill in estimates of Qbf into 
# a table. 
# 3) From there, this script is run again and estimates of channel geometry at bankfull
# discharge are made.
#
# NOTE: script must be altered a bit to run twice: First time to create plots to estimate
# bankfull discharge. Second time to estimate channel geometry at bankfull discharge.
#
# Joel Sholtes, PhD Candidate, Colorado State University, August 2014
#######################################################################################

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

# fucntion to return a vector for placing minor tick locations.
# int: # of ticks total (major + minor) from beginning major tick to final minor tick
source("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/minor.ticks.axis.R")

# import table of USGS site numbers, this table should have a column labeled "Qbf.est".
# after the first run, estimates of Qbf are entered here and the script is run again
# to calculate the width and average depth at Qbf
site.table <- read.table("susp_sites_v10.csv", header = TRUE, sep = ",", 
                         stringsAsFactors = F)
sites <- site.table[c("site.no", "site.nm")]
# append GSD variables as columns to site table
sites[c("n", "Qbf.m", "Wbf.m", "Dbf.m")] <- NA

pdf( "XS_Plots_Qbf.pdf", width = 8.5, height = 11, paper = "special") # turn on PDF creator
  # set plot to be square and 4x4 plots
  par(mfrow = c(2,2), 
      pty = "s", 
      mar = c(5.1,5.1,1,1),
      oma = c(2,2,2,2),
      cex = 1, cex.lab = 1.1, cex.axis = 1, cex.main = 1)
  
  for( i in 1:nrow(sites)) {
    site.no <- sites$site.no[i]
    # attach leading zero to site.no if necessary
    if(nchar(site.no) < 8){
      site.no <- paste("0", site.no, sep = "")
    }
    # Create URL to download site-specific field measurements
    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)
    # Check to see if any data came back, then proceed
    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
      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))]
      h <- as.numeric(xs$gage_height_va)
      Q <- as.numeric(xs$discharge_va)
      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)]
      
      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)])
      
      # NOTE: This comes in step 2 after initial Qbf estimates have been made (Qbf.est)
      # Calc Wbf and Dbf based on Qbf estimate and linear interpolation
      Qbf <- site.table$Qbf.est[i]
      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
      }
        
      ## Write results to table
      sites[i,3:6] <- c(n, Qbf.m, wbf.m, dbf.m)
      
      ## 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()
      
      mtext(paste(sites$site.nm[i], "\n", site.no), outer = T, side = 3, line = -2, font = 2)
      
      ## 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 = expression("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()
    }
  }

dev.off()

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