##################################################################################
## ImportDV.R Script authored by Ryberg and Vecchia, USGS 2012
## This function imports daily flow data from USGS NWIS website
## Required input for this function are "staid" (station id as character), 
## and start and end dates. If you input a very early date ("1880-10-01"), 
## NWIS will return the earliest recorded flow. 
##################################################################################

importDV <- function (staid, code = "00060", stat = "00003", sdate, 
                      edate) 
{
  if (is.character(staid) == FALSE) 
    stop("staid needs to have quotes around it")
  if (nchar(staid) < 8) 
    stop("staid must be at least 8 characters")
  base_url <- "http://waterservices.usgs.gov/nwis/dv?"
  url <- paste(base_url, "site=", staid, "&parameterCd=", code, 
               "&statCd=", stat, sep = "")
  url <- paste(url, "&startDt=", sdate, "&endDt=", edate, sep = "")
  doc <- xmlTreeParse(url, getDTD = FALSE, useInternalNodes = TRUE)
  r <- xmlRoot(doc)
  i <- 1
  val <- vector(mode = "numeric", length = 1)
  while (xmlName(r[[2]][[3]][[i]]) == "value") {
    val[i] <- as.numeric(xmlValue(r[[2]][[3]][[i]]))
    i <- i + 1
  }
  Attribute <- xmlApply(r[[2]][[3]], xmlAttrs)
  N <- length(val)
  NoDataValue <- xmlValue(r[["timeSeries"]][["variable"]][["NoDataValue"]])
  NoDataValue <- as.integer(NoDataValue)
  dates <- vector(mode = "character", length = 1)
  qualcode <- vector(mode = "character", length = 1)
  if (N > 1) {
    for (z in 1:N) {
      dates[z] <- as.character(strsplit(Attribute[z][[1]][[2]], 
                                        "T")[[1]][1])
      qualcode[z] <- Attribute[z][[1]][[1]]
    }
    dates <- as.Date(dates, "%Y-%m-%d")
    df <- data.frame(staid, val, dates, qualcode)
    beginDate <- df$dates[1]
    endDate <- df$dates[dim(df)[[1]]]
    myDates <- as.data.frame(seq.Date(beginDate, endDate, 
                                      by = 1))
    dimnames(myDates)[[2]][1] <- "dates"
    ndays <- dim(myDates)[1]
    nobs <- dim(df)[1]
  }
  attributes(df)$code <- code
  attributes(df)$stat <- stat
  df
}