################################################################################
# Script to download instantaneous maximum peak discharge data from USGS and 
# calculate the Q1.5, Q2, Qavg, and Qmed floods from this series. Q1.5 and Q2
# are calculated via linear interpolation using the Weibull Plotting Method.
#
# Joel S Sholtes - PhD Candidate
# Colorado State University - Dept. Civil & Env. Engr
# May, 2014
# jsholtes@gmail.com
################################################################################

# Clear data
rm(list = ls())

# set working dir
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/SuspLoad")

# import table of USGS site numbers for Suspended Sites
site.table <- read.table("susp_sites_v10.csv", header = TRUE, sep = ",", 
                         stringsAsFactors = F)
sites <- site.table[c("site.no", "site.nm")]

# import table of USGS site numbers for Suspended Sites
# sites <- read.table("bedload_sites.csv", header = TRUE, sep = ",", 
#                          stringsAsFactors = F)

# Append new peak flow metric columns to site table
sites[c("n.yrs", "Qavg", "Qmed", "Q1.5", "Q2")] <- NA

# Loop to read site ID, download peak discharge data, and calc peak flow values
for (i in 1:nrow(sites)){
  
  site.no <- sites$site.no[i]
  
  # some site IDs are made up (extended flow records) therefore no
  # no USGS peak data will be available
  if( !is.na( as.numeric(site.no) ) ) { 
  
  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 = "")
  
  
  # Web query URL
  path <- paste("http://nwis.waterdata.usgs.gov/nwis/peak?site_no=",
                site.no, "&agency_cd=USGS&format=rdb", sep = "")
  
  # reads in tab-delimited RDB file as dataframe and skips lines beginning with "#"
  pk.table <- read.table(path, comment.char = "#", header = T, sep = "\t", 
                     stringsAsFactors = F)
  
    # Check to see if any data from NWIS came back
    if(!is.null(dim(pk.table)) & ncol(pk.table) > 2) { 
      # Remove first row, which does not contain data
      pk.table <- pk.table[-1,]
      # convert these columns from characters to text
      pk.table$peak_va <- as.numeric(pk.table$peak_va)
      pk.table$gage_ht <- as.numeric(pk.table$gage_ht)
      
      # Weibull Plotting Method
      peak <- sort(pk.table$peak_va, decreasing = T)
      rank <- seq(1,length(peak), by = 1)
      prob <- rank / (length(peak) + 1)
      RI <- 1 / prob
      
      Qavg <- round( mean(peak), 0)
      Qmed <- round( median(peak), 0)
      
      # Calc Q1.5 and Q2 using linear interpolation from Weibull Plotting Series
      Q1.5 <- round( approx(RI, peak, xout = 1.5)$y, 0)
      Q2 <- round( approx(RI, peak, xout = 2)$y, 0)
      n.yrs <- nrow(pk.table)
      
      sites[i,3:7] <- c(n.yrs, Qavg, Qmed, Q1.5, Q2)
      
      # Write table of peak discharge values and dates in new "peaks/" folder
      peak.table <- data.frame(pk.table$peak_dt, pk.table$peak_va)
      names(peak.table) <- c("date", "peak.cfs")
      write.csv(peak.table, paste("peaks/" , site.no, ".csv", sep = ""), 
                row.names = F) # Need to make this folder first
    }
  }
}

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