################################################################################
##  Import & Save USGS Flow Data - Daily or Instantaneous (15 min)
##  Joel S Sholtes - PhD Candidate
##  jsholtes@gmail.com
##  Colorado State University, Dept. of Civl and Env. Engr.
##  January 2014
##
##  This script reads a table within USGS station and date information and uses
##  the importDV function (or importIV) to query NWIS, download the data into a
##  data frame, and save it in a specified folder. This function works for any
##  number of sites as input into the table. It requires the XML package along
##  with the ImportDV function, or the waterData package and the importDVs func.
##
##  Note 1: for gages with shorter records to be extended with regression with 
##  nearby gage, maybe need to go back and download original dataset because
##  this script finds first instance of YYYY-10-01 to trim the flow file. 
##  In some cases, seasonal flow records have mnay summers of data before a 
##  YYYY-10-01 records appears. This is flow data that would be lost otherwise
##  for flow extension regression analysis.
##  
##  Note 2: **This is not the most robust code** 
##  Known Issues: 1) Sometimes NWIS website doesn't respond and no flow file
##  is downloaded. Just start the loop back at the i value it stopped and try 
##  again. 
##  2) Sometimes the start date while loop and if/else statements become
##  stuck on one start date value. Can't figure this out. It usually works if 
##  you set the i value to the site no you want to download and run it manually
##  inside the loop (select code inside of for loop and CTRL+R). 
##  3) If no complete years exist on the flow record or NWIS doesn't produce a 
##  flow file, then there will be a gap in the site record causing problems in 
##  the loops further down. Just be aware, can set i to start after that problem
##  record and re-start loop.
###############################################################################

# Clear data
rm(list = ls())
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis")
# load libraries
library(waterData) # fill missing year function... doesn't work very well!
# library(moments) 
library(chron) # leap.year()

# load custom water download function modified from importDVs{waterData}
# can also use the importDVs function from waterData instead.
source("ImportDV.R")

################################################################################
##  Initial import and cleaning of flow data from NWIS
##  After this, run subsequent code (below) to remove partial years
################################################################################

## Set working dir for bedload or suspended load

# Suspended Load Data
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/SuspLoad")
path <- "susp_sites_v10.csv"
sites <- read.table(path, header = TRUE, sep = ",", stringsAsFactors = F)

# # Bed Load Data
# setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/Bedload")
# path <- "additional_sites.csv" 
# sites <- read.table(path, header = TRUE, sep = ",", 
#                     colClasses = c("character", "character"))

# create table to summary info for each site
cnames <- c("site.no", "site.nm", "start.date", "end.date", "yrs", "yrs.miss")
flow.summary <- data.frame(matrix(data = 0, nrow = nrow(sites), ncol = length(cnames)),
                      row.names = NULL, stringsAsFactors = F)
names(flow.summary) <- cnames  # Assigns col names to data frame

# # Import existing partial flow.summary table?
# flow.summary <- read.table("flow.summary2.csv",sep = ",", stringsAsFactors = F, header = T)

# 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

for ( i in 1:nrow(sites) ) {
  staid <- sites[i,1] # Station ID colums (character)
  #sdate <- sites[i,6] # Start date column (character)
  #edate <- sites[i,7] # End date column (character)
  sdate <- "1890-10-10"
  edate <- "2014-09-30"
  if( nchar(staid) < 8) {
    staid <- paste("0", staid, sep = "")
  } else{
    staid <- as.character(staid)
  }
  
  flow.file <- importDV(staid, code = "00060", stat = "00003", sdate, edate) 
  
#   ## Use this if you just want to import and edit a flow file
#   staid <- "13011000"
#   path <- paste("DV_Flows/", staid, ".csv", sep = "")
#   flow.file <- read.table(path, header = T, stringsAsFactors = F, sep = ",")

  # re-format "dates" col as.Date()
  flow.file$dates <- as.Date(flow.file$dates, format = "%Y-%m-%d")
  
  # Remove negative entries and replace with NA
  flow.file$val <- with(flow.file, replace(flow.file$val, flow.file$val < 0, NA))
  
  # Remove leap year entries
  flow.file <- subset(flow.file, !(format(flow.file$date, "%m") == "02" & 
                                     format(flow.file$date, "%d") == "29"))
  
  # start flow.file at YYYY-10-01
  # if more than 274 days in start date (10-02 is 275th day) then start on *following* year
  # 30.3333 days per month Jan through Sept 30th
  start.days <- (as.numeric(format(flow.file$dates[1], "%m")) - 1)*30.3333 + 
    as.numeric(format(flow.file$dates[1], "%d"))
  
  if ( start.days > 274) {
    start.yr <- as.numeric(format(flow.file$dates[1], "%Y")) + 1
    start.date <- as.Date(paste(start.yr,"-10-01", sep = ""))
    start.idx <- which(start.date == flow.file$dates)
  } else {
    # if false -> start *current* year at 10-01
    start.yr <- as.numeric(format(flow.file$dates[1], "%Y"))
    start.date <- as.Date(paste(start.yr,"-10-01", sep = ""))
    start.idx <- which(start.date == flow.file$dates)
    }
  
  # Sometimes the date YYY-10-01 does not exist in the first year of record,
  # so we have to cycle through the years untile we find one and start.idx contains a value
  if (length(start.idx) == 0) {
    k = 1
    while ( length(start.idx) == 0 ) {
      if ( start.days > 274) {
        start.yr <- as.numeric(format(flow.file$dates[k], "%Y")) + 1
        start.date <- as.Date(paste(start.yr,"-10-01", sep = ""))
        start.idx <- which(start.date == flow.file$dates)
      } else {
        # if false -> start *current* year at 10-01
        start.yr <- as.numeric(format(flow.file$dates[k], "%Y"))
        start.date <- as.Date(paste(start.yr,"-10-01", sep = ""))
        start.idx <- which(start.date == flow.file$dates)
      }
      
      k = k + 365 # jump to next flow year (assuming the next year has a whole year of data...)
      
      # break while loop if we go beyond teh flow record (i.e., no YYYY-10-01 dates on record)
      if ( is.na(start.yr) ) {
        # this just sets start.date to first entry and I'll have to go back and manually check flow record
         start.yr <- as.numeric(format(flow.file$dates[1], "%Y")) + 1
         start.date <- as.Date(paste(start.yr,"-10-01", sep = ""))
         start.idx <- 1 
        break
      }
    }
  }
  
  # Extract end date from flow.file
  end.date <- flow.file$dates[nrow(flow.file)]
  end.days <- (as.numeric(format(end.date, "%m")) - 1)*30.375 + as.numeric(format(end.date, "%d"))
  
  # End flow.file at 09-30
  # If end.date is >= 09-30 (273 days), then end *that* year at 09-30, else end in *previous* year
  if ( end.days >= 273 ) {
    end.yr <- as.numeric(format(end.date, "%Y"))
    end.date <- as.Date(paste(end.yr,"-09-30", sep = ""))
    end.idx <- which(end.date == flow.file$dates)
  } else {
    end.yr <- as.numeric(format(end.date, "%Y")) - 1
    end.date <- as.Date(paste(end.yr,"-09-30", sep = ""))
    end.idx <- which(end.date == flow.file$dates)
    }
 
  # Subset flow.file using start and end date indices
  flow.file <- flow.file[start.idx:end.idx,]
  
  # Fill in any NA values with waterData interpolation tool <- doesn't work well, so no fill
  fill = F
  
  # Save flow data file
  file.name <- staid
  write.table(flow.file, file = paste("DV_Flows/", file.name, ".csv", sep = ""), 
              row.names = FALSE, sep = ",")
  
  # Write row of flow.summary results by site
  years <- end.yr - start.yr
  years.diverge <- nrow(flow.file)/365 - years
  days.diverge <- nrow(flow.file) - (end.yr - start.yr)*365
  
  # write results as a row in the flow.summary table
  flow.summary[i,] <- c(sites$site.no[i], sites$site.nm[i], as.character(start.date), 
                        as.character(end.date), years, years.diverge)
}  

write.table(sites, file = "Susp_Sed_Sites_v8.csv", row.names = F, sep = ",")

write.table(flow.summary, file = "Flow.Summary_v7.csv", 
            row.names = FALSE, sep = ",")

########################################################################
## Script to Remove years with partial data
## Doing this separately after downloading all flow data
## Need to run this again after incorporating extended flow records
########################################################################

setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/SuspLoad")
path <- "Susp_Sed_Sites_v8.csv"
sites <- read.table(path, header = TRUE, sep = ",", stringsAsFactors = F)
site.no <- sites$site.no
site.no <- ifelse( nchar(site.no) < 8, paste("0", site.no, sep = ""), site.no) 
site.nm <- sites$site.nm

# setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/Bedload")
# path <- "bedload_sites.csv"
# sites <- read.table(path, header = TRUE, sep = ",", 
#                     colClasses = c("character", "character"))
# site.no <- sites$site.no
# site.nm <- sites$site.nm

## Create table that shows that alll flow data are divisible by 365 to 
# verify partial year culling from above & create table to summary info for each site
cnames <- c("site.no", "site.nm", "start.date", "end.date", "yrs", "yrs.miss")

flow.summary2 <- data.frame(matrix(data = 0, nrow = nrow(sites), ncol = length(cnames)),
                            row.names = NULL, stringsAsFactors = F)
names(flow.summary2) <- cnames  # Assigns col names to data frame

for (i in 1:nrow(sites)) {
  setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/SuspLoad/DV_Flows2")
  path <- paste(site.no[i], ".csv", sep = "")
  data <- read.table(path, header = T, stringsAsFactors = F, sep = ",")
#   setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/SuspLoad/DV_Flows2")
  
  # Remove leap year entries (only in USFS files with 2 cols)
  if (ncol(data) < 4){
  data$date <- as.Date(data$date, format = "%Y-%m-%d")
  data <- subset(data, !(format(data$date, "%m") == "02" & 
                           format(data$date, "%d") == "29"))
  }
  if (ncol(data) >= 4){
    data$dates <- as.Date(data$dates, format = "%Y-%m-%d")
    data <- subset(data, !(format(data$dates, "%m") == "02" & 
                             format(data$dates, "%d") == "29"))
  }
  
  # this checks if there are whole years with NAs and takes them out
  if ( nrow(data)/365 == floor(nrow(data)/365) &
      sum(is.na(data$val))/365 == floor(sum(is.na(data$val))/365) ){
      data <- data[!is.na(data$val),]
      write.table(data, file = paste(site.no[i], ".csv", sep = ""),
                  row.names = FALSE, sep = ",")
      if( ncol(data) > 3 ) {
        data$dates <- as.Date(data$dates, format = "%Y-%m-%d")
        data <- data[!is.na(data$val),]
        WY <- numeric(nrow(data))
        data <- data.frame(data, WY)
        start.yr <- as.numeric(format(data$dates[1], "%Y"))
        start.date <- as.Date(paste(start.yr,"-10-01", sep = ""))
        end.date <- data$dates[nrow(data)]
        end.yr <- as.numeric(format(end.date, "%Y"))
        days <- (end.yr - start.yr)*365
        yrs <- nrow(data)/365
        years.miss <- (days - nrow(data))/365
      }else{ # for flow files with only 2 columns
        data$date <- as.Date(data$date, format = "%Y-%m-%d")
        data <- data[!is.na(data$val),]
        WY <- numeric(nrow(data))
        data <- data.frame(data, WY)
        start.yr <- as.numeric(format(data$date[1], "%Y"))
        start.date <- as.Date(data$date[1])
        end.date <- data$date[nrow(data)]
        end.yr <- as.numeric(format(end.date, "%Y"))
        days <- (end.yr - start.yr)*365
        years.miss <- (days - nrow(data))/365
        yrs <- nrow(data)/365
        years.miss <- (days - nrow(data))/365
      }
    }else{
    if( ncol(data) > 3 ) {
      data$dates <- as.Date(data$dates, format = "%Y-%m-%d")
      data <- data[!is.na(data$val),]
      WY <- numeric(nrow(data))
      data <- data.frame(data, WY)
      start.yr <- as.numeric(format(data$dates[1], "%Y"))
      start.date <- as.Date(paste(start.yr,"-10-01", sep = ""))
      end.date <- data$dates[nrow(data)]
      end.yr <- as.numeric(format(end.date, "%Y"))
      days <- (end.yr - start.yr)*365
      yrs <- nrow(data)/365
      years.miss <- (days - nrow(data))/365
      
      # this loop first checks to see if nrows = 365 * # of years
      # is not (i.e., if days missing in a year), then:
      if( nrow(data)/365 != floor(nrow(data)/365) ) {
          # 1) we assign a water year to each day accounting for leap years (which have been removed)
          for (j in 1:nrow(data)) {
            if( leap.year(as.numeric(format(data$dates[j], "%Y"))) ) {
              if(strptime(data$dates[j], "%Y-%m-%d")$yday+1 >= 275){
                data$WY[j] = as.numeric(format(data$dates[j], "%Y")) + 1
              } else {
                data$WY[j] = as.numeric(format(data$dates[j], "%Y"))
                }
            } else{
                if(strptime(data$dates[j], "%Y-%m-%d")$yday+1 >= 274){
                  data$WY[j] = as.numeric(format(data$dates[j], "%Y")) + 1
                } else {
                  data$WY[j] = as.numeric(format(data$dates[j], "%Y"))
                  }
              }
          }
          
          # 2) creates data frame with count of rows per water year -> table() does this
          WY.count <- as.data.frame(table(data[,5]))
          bad.yrs <- subset(WY.count, WY.count$Freq < 365)
          
          # 3) remove any years where 365 entries do not exist
          if (nrow(bad.yrs) > 0) {
            data <- subset(data, !data$WY %in% bad.yrs$Var1)
          }
          # 4) re-write flow file with partial years removed
          write.table(data, file = paste(site.no[i], ".csv", sep = ""),
                      row.names = FALSE, sep = ",")
          # 5) If there were any bad years, then we write a table with those years in it.
          if (nrow(bad.yrs) > 0) {
          write.table(bad.yrs, file = paste(site.no[i], "_Bad_Yrs_NA.csv", sep = ""),
                      row.names = FALSE, sep = ",")
          }
          # Re-Calc these variables to fill in table
          start.date <- as.Date(data$dates[1])
          end.date <- data$dates[nrow(data)]
          days <- (end.yr - start.yr)*365
          yrs <- nrow(data)/365
          years.miss <- (days - nrow(data))/365
        }
      }else{ # for flow files with only 2 columns
          data$date <- as.Date(data$date, format = "%Y-%m-%d")
          data <- data[!is.na(data$val),]
          WY <- numeric(nrow(data))
          data <- data.frame(data, WY)
          start.yr <- as.numeric(format(data$date[1], "%Y"))
          start.date <- as.Date(data$date[1])
          end.date <- data$date[nrow(data)]
          end.yr <- as.numeric(format(end.date, "%Y"))
          days <- (end.yr - start.yr)*365
          years.miss <- (days - nrow(data))/365
          yrs <- nrow(data)/365
          years.miss <- (days - nrow(data))/365
          
          # this loop first checks to see if nrows = 365 * # of years
          # is not (i.e., if days missing in a year), then:
          
          if( nrow(data)/365 != floor(nrow(data)/365) ) {
            # 1) we assign a water year to each day accounting for leap years (which have been removed)
            for (j in 1:nrow(data)) {
              if( leap.year(as.numeric(format(data$date[j], "%Y"))) ) {
                if(strptime(data$date[j], "%Y-%m-%d")$yday+1 >= 275){
                  data$WY[j] = as.numeric(format(data$date[j], "%Y")) + 1
                } else {
                  data$WY[j] = as.numeric(format(data$date[j], "%Y"))
                }
              } else{
                if(strptime(data$date[j], "%Y-%m-%d")$yday+1 >= 274){
                  data$WY[j] = as.numeric(format(data$date[j], "%Y")) + 1
                } else {
                  data$WY[j] = as.numeric(format(data$date[j], "%Y"))
                }
              }
            }
            
            # 2) creates data frame with count of rows per water year -> table() does this
            WY.count <- as.data.frame(table(data[,3]))
            bad.yrs <- subset(WY.count, WY.count$Freq < 365)
            
            # 3) remove any years where 365 entries do not exist
            if (nrow(bad.yrs) > 0) {
              data <- subset(data, !data$WY %in% bad.yrs$Var1)
            }
            # 4) re-write flow file with partial years removed
            write.table(data, file = paste(site.no[i], ".csv", sep = ""),
                        row.names = FALSE, sep = ",")
            # 5) If there were any bad years, then we write a table with those years in it.
            if (nrow(bad.yrs) > 0) {
              write.table(bad.yrs, file = paste(site.no[i], "_Bad_Yrs_NA.csv", sep = ""),
                          row.names = FALSE, sep = ",")
            }
          }
          # Re-Calc these variables to fill in table
          start.date <- as.Date(data$date[1])
          end.date <- data$date[nrow(data)]
          days <- (end.yr - start.yr)*365
          yrs <- nrow(data)/365
          years.miss <- (days - nrow(data))/365
      }
    }
  
    flow.summary2[i,] <- c(site.no[i], site.nm[i], as.character(start.date), 
                             as.character(end.date), yrs, years.miss)
}  

# Note: Years.miss should be an integer meaning complete years missing
write.table(flow.summary2, file = "flow.summary2.csv", row.names = FALSE, sep = ",")
