# This script creates time slices for a list of USGS gauge

library(rwrfhydro)
library(dataRetrieval)
library(foreach)
outPath <- "./TimeSlices"
#***********************************************************************************************************************
# get the instantaeous data from the website 
#***********************************************************************************************************************
# Provide the list of the gages 
assimilate_gages <- read.csv("./list_of_assimilation_gages.csv", header = TRUE, colClasse = "character")

print(assimilate_gages)
               
discharge <- foreach (siteNumber = assimilate_gages, .combine = rbind.data.frame) %do% {
  parameterCd <- "00060"  # Discharge
  startDate <- "2019-12-31"
  endDate <- "2020-12-31"
  
  discharge <- readNWISuv(siteNumber, parameterCd, startDate, endDate)
  return(discharge)
}

#************************************************************************************************************************
#   Prepare the format for the USCG time slices write
#************************************************************************************************************************

# add the discharge values in cms
discharge$discharge.cms <- discharge$X_00060_00000 * 0.028316847

# add the quality 
Qc <- function(dataDf) {
  ## assume the worst and then recover what fits out mental model of this chaos.
  ## JLM: a different assumption off the bat.
  dataDf$quality <- dataDf$discharge.cms * 0
  
  ## Recover quality for valid/sane flow range
  ## JLM this deviates from their script
  ## in their script they consider 0 a valid flow, but I'd argue otherwise...
  ## why would any streams with no flow be gaged?? Then again, drought. IT's a really dicey value,
  ## JLM also limit maximum flows to approx twice what I believe is the largest gaged flow
  ## on MS river *2
  ## http://nwis.waterdata.usgs.gov/nwis/peak?site_no=07374000&agency_cd=USGS&format=html
  ## baton rouge 1945: 1,473,000cfs=41,711cms
  ## multiply it roughly by 2
  
  isValidFlow <- dataDf$discharge.cms > 0 & dataDf$discharge.cms < 90000
  wh100 <- which(isValidFlow)
  if(length(wh100)) dataDf$quality[wh100] <- 100
}
discharge$quality <- Qc(discharge)

# add query time to the data
queryTime <- Sys.time()
attr(queryTime, "tzone") <- "UTC"
discharge$queryTime <- queryTime   # add the query time

# reformat the name of th gauge to 15 character
discharge$site_no <- formatC(discharge$site_no, width=15)


# add the Round minutes to the discharge data frame 
RoundMinutes <- function (POSIXct, nearest = 5) {
  if ((60%%nearest) != 0)
    warning(paste0("The nearest argument (passed: ", nearestMin,
                   ") is mean to divide 60 with no remainder."), immediate. = TRUE)
  nearestInv <- 1./nearest
  theMin <- as.numeric(format(POSIXct, "%M")) + as.numeric(format(POSIXct, "%S"))/60
  floorDiff <- (theMin - nearest * (floor(theMin/nearest))) / nearest # added by Arezoo
  whFloor <- which(floorDiff < 0.5)
  roundMin <- (ceiling(theMin * nearestInv)/nearestInv)
  roundMin[whFloor] <- (floor(theMin * nearestInv)/nearestInv)[whFloor]
  diffMin <- roundMin - theMin
  lubridate::floor_date(POSIXct, "hour") + lubridate::minutes(floor(roundMin))
}

discharge$dateTimeRound <- RoundMinutes(discharge$dateTime,nearest=15)


discharge <- subset(discharge, lubridate::hour(dateTimeRound) == 23 & lubridate::minute(dateTimeRound) == 0)
print(head(discharge))
discharge$dateTimeRound <- discharge$dateTimeRound - 23*3600
discharge$dateTime <- discharge$dateTime - 23*3600
print(head(discharge))
#*********************************************************************************************************************************
# Now it is time to write it into time slices
#*********************************************************************************************************************************

#Loop through the times and and write then into files
for (i in 1:length(rev(sort(unique(discharge$dateTimeRound))))) {
  
  dfByPosix <- subset(discharge,dateTimeRound == rev(sort(unique(discharge$dateTimeRound)))[i])
  dfByPosix$discharge.quality <- dfByPosix$quality
  rwrfhydro::WriteNcTimeSlice(dfByPosix,
                   outPath=outPath,
                   sliceResolution = 15)
}