# Clear everything. Clean start -------------------------------------------
rm(list = ls())

# Hydrological data retrieval ------------------------------------------------------
library(EGRET)

# Flow history analysis
# Gather discharge data: Unit in SI
siteNumber <- "12422500" # SPOKANE RIVER AT SPOKANE, WA
startDate <- "1990-10-01" 
endDate <- "2005-09-30"
Q_daily <- readNWISDaily(siteNumber, "00060", startDate, endDate)

# Check data
str(Q_daily)
head(Q_daily)
summary(Q_daily)

# Create a year column
library(lubridate)
Q_daily$Year <- year(Q_daily$Date)

# Gather site and parameter information:
# Here user must input some values for
# the default (interactive=TRUE)
INFO <- readNWISInfo(siteNumber, "00060")

# Data Visualization ------------------------------------------------------
library(xts)
Q_daily_xts <- xts(x = (Q_daily[, 2]), order.by = (Q_daily$Date))
str(Q_daily_xts)
plot(Q_daily_xts)

library(hydroTSM)
hydroplot(Q_daily_xts, var.type = "Flow", main = "HydroTSM",
          pfreq = "dma", var.unit = "ft3", from = startDate)

hydroplot(Q_daily_xts, var.type = "Flow", main = "HydroTSM",
          pfreq = "seasonal", var.unit = "ft3", from = startDate)

library(dygraphs)
dygraph(Q_daily_xts, main = "Q: SPOKANE RIVER AT SPOKANE") %>%
  dyOptions(fillGraph = TRUE, fillAlpha = 0.4)

dygraph(Q_daily_xts, main = "Q: SPOKANE RIVER AT SPOKANE") %>%
  dyOptions(fillGraph = TRUE, fillAlpha = 0.4) %>% 
  dyRangeSelector()

# Data Analysis -----------------------------------------------------------
# Check flow history data:
eList <- as.egret(INFO, Q_daily, NA, NA)
eList
str(eList)
# A part of the flowHistory system. The index of the flow statistics is istat. These statistics are: 
#   (1) 1-day minimum, (2) 7-day minimum, (3) 30-day minimum, (4) median 
#   (5) mean, (6) 30-day maximum, (7) 7-day maximum, and (8) 1-day maximum
plotFlowSingle(eList, istat = 3, qUnit = "cfs")
plotFourStats(eList, qUnit=3)

# Seasonal average
create_season_df <- function(data_frame_object) {
  data_frame_object$Season <- with(data_frame_object,
                                    ifelse((Month == 10 | Month == 11 | Month ==
                                    12), "Fall", ifelse((Month == 1 | Month ==
                                    2 | Month == 3), "Winter", ifelse((Month ==
                                    4 | Month == 5 | Month == 6), "Spring",
                                    "Summer"))))
  return(data_frame_object)
}

Q_daily <- create_season_df(Q_daily)
head(Q_daily)

library(dplyr)
Q_daily_seasonal_ave <- Q_daily %>% 
  tbl_df() %>% 
  select(Q, Season, Year) %>% 
  group_by(Year, Season) %>% 
  summarise_each(funs(mean = mean(., na.rm = TRUE)))
Q_daily_seasonal_ave

library(ggplot2)
q_plot <- ggplot(data = Q_daily_seasonal_ave, aes(x = Year, y = Q))
q_plot + facet_grid(~ Season) + geom_line() + theme_grey()
last_plot() + theme(axis.text.x = element_text(angle = 90))

################################################################################
#### Model result evaluation in R ####
################################################################################
# Load package EcoHydRology for Baseflow Separation
library(hydroGOF)

obs <- data.frame(Q_daily$Date, Q_daily$Q)
sim <- data.frame(Q_daily$Date, 0.9*Q_daily$Q)

library(zoo)
Observed <- read.zoo(obs[, 1:2], drop = FALSE)
Simulated <- read.zoo(sim[, 1:2], drop = FALSE)

#### Getting the new numeric goodness-of-fit measures
gof(sim = Simulated, obs = Observed)

# Plot 'obs' vs 'sim' for the daily, monthly and annual time series 
# graphics.off()
ggof(sim = Simulated, obs = Observed, ylab = "Q (ft3/day)",
     ftype = "dma", FUN = mean)