############################################################################################
# FlowMetrics.R <- A script that analyses flow records and calculates various statistical
# properties of them. This script also calculates the Median Annual Flow Duration Curve
# following Vogel and Fennessey (1994)
#
# Joel Sholtes, PhD Candidate, Colorado State University, March 2014
############################################################################################

rm(list = ls())
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/SuspLoad")
path <- "susp_sites_v10.csv"
library(Hmisc) # for hdquantile estimator
library(moments) # skewness

## set length of percentile vector for AFDC
p.length <- 1000
# Assigns lognormal (?) distributed percentile probability of exceedence
# based on n = 365 days
# p <- pnorm(seq(qnorm(1/(365+1)), qnorm(.9999), length.out = 100))
p <- pnorm(seq(qnorm(0.0001), qnorm(.9999), length.out = p.length))

sites <- read.table(path, header = TRUE, sep = ",", stringsAsFactors = F)
site.no <- ifelse(nchar(sites$site.no) < 8, paste("0", sites$site.no, sep = ""),
                  sites$site.no)
site.nm <- sites$site.nm

# Create empty dataframe to put flow metric results in
cnames <- c("site.no", "site.nm", "yrs", "min", "max", "mean", "med", "var", "cv", 
            "spread", "skewness", 
            "Q1", "Q1.5", "Q2", 
            "Q1:Qavg", "Q1:mean", "Q1.5:Qavg", "Q1.5:mean", "Q2:Qavg", "Q2:mean",
            "flash.RB", "Flash.Tavg")
results <- data.frame(matrix(data = 0, nrow = length(site.no), ncol = length(cnames)),
                      row.names = NULL, stringsAsFactors = F)
names(results) <- cnames  # Assigns col names to data frame

# Create empty data frame to place median AFDC and 95 p-tile AFDCs
afdc <- data.frame(matrix(data = 0, nrow = length(site.no), ncol = p.length),
                   row.names = site.no, stringsAsFactors = F)

afdc.975 <- data.frame(matrix(data = 0, nrow = length(site.no), ncol = p.length),
                       row.names = site.no, stringsAsFactors = F)

afdc.025 <- data.frame(matrix(data = 0, nrow = length(site.no), ncol = p.length),
                       row.names = site.no, stringsAsFactors = F)

for (i in 1:length(site.no)) {
  path <- paste("DV_Flows2/", site.no[i], ".csv", sep = "")
  data <- read.table(path, header = T, stringsAsFactors = F, sep = ",")
  flow <- data$val*0.028316847
  nyears <- (length(flow)/365)
  
  # Calc daily stats 
  # *** consider going back and making these "L-moments" 
  min <- round(min(flow),3)
  max <- round(max(flow),3)
  mean <- round(mean(flow),3)
  med <- round(median(flow),3)
  var <- round(var(flow),3)
  cv <- round(sqrt(var)/mean,3)
  spread <- round( ( as.numeric(quantile(flow, 3/4)) - as.numeric(quantile(flow, 1/4)) ) / med, 3 )
  skewness <- round(skewness(flow),3)
  flash.RB <- round(sum(abs(diff(flow)))/sum(flow) , 3)  # Richards- Baker (2004) Flashiness index
  flash.Tavg <- round( sum(flow > mean)/length(flow) , 3) # Double check with Sarah
  
  # calc annual max flood stats based on daily flow record
  # Create matrix with years as columns and 365 rows
  q.mat <- matrix(flow, ncol = nyears, nrow = 365, byrow = FALSE)
  # Sort matrix so that flows are descending by row
  #   sort.q <- apply(q.mat, 2, sort, decreasing = T)
  q.pk <- apply(q.mat, 2, max) # creates vector max annual flood series. apply(x, by, fun): by-> 1 row, 2 col
  q.pk <- sort(q.pk, decreasing = T) # ranks with largest flood first
  rank <- seq(1, nyears, 1)
  prob <- rank / ( nyears + 1 ) # creates weibull probs vector
  RI <- 1/prob
  
  # now find the index of the flows associated with RI = 1, 1.5, 2
  Q1 <- round(q.pk[which(abs(RI-1) == min(abs(RI-1)))],3)
  Q1.5 <- round(q.pk[which(abs(RI-1.5) == min(abs(RI-1.5)))],3)
  Q2 <- round(q.pk[which(abs(RI-2) == min(abs(RI-2)))],3)
  Qavg <- round(mean(q.pk),3)
  Q1.Qavg <- round(Q1.5/Qavg,3)
  Q1.mean <- round(Q1.5/mean,3)
  Q1.5.Qavg <- round(Q1.5/Qavg,3)
  Q1.5.mean <- round(Q1.5/mean,3)
  Q2.Qavg <- round(Q1.5/Qavg,3)
  Q2.mean <- round(Q1.5/mean,3)
  
  results[i,] <- c(site.no[i], site.nm[i], nyears, min, max, mean, med, var, cv, 
                   spread, skewness, 
                   Q1, Q1.5, Q2, Q1.Qavg, Q1.mean, Q1.5.Qavg, Q1.5.mean,Q2.Qavg, Q2.mean,
                   flash.RB, flash.Tavg)

    ## AFDC Analysis
    # Harrell Davis quantile estimator (smoothest) -> hdquantile is the function that does this
    # hmisc package
    # calcs value of p-tiles for each year (individual AFDCs) by column
    # 'qp3' refers to Vogel and Fennessey (1993)'s notation for the HD quantile estimator (percentile?)
    qp3 <- matrix(data = 0, nrow = length(p), ncol = nyears)
      for (j in 1:nyears){
        qp3[,j] <- hdquantile(q.mat[,j], probs = p,
                              se = FALSE, na.rm = FALSE, names = FALSE, weights=FALSE)
      }
    
    # Calc median, 0.025 and 0.0975 CI bands based on nyears of Qp3 AFDCs
    # the 97.5%, 0.25% and median (50%) quantile AFDCs are in this matrix
    # the median AFDC is the q3_tile[,2] -> the middle columns of this matrix
    q3_tile <- matrix(data = 0, ncol = 3, nrow = length(p))
      for (j in 1:length(p)){
        q3_tile[j,] <- hdquantile(qp3[j,], probs = c( 0.025, 0.5, 0.975), 
                                  se = FALSE, na.rm = FALSE, names = FALSE, weights=FALSE)
      }
    
    # Write AFDCs into tables by row
    afdc.025[i,] <- q3_tile[1,]  
    afdc[i,] <- q3_tile[2,]
    afdc.975[i,] <- q3_tile[3,]
}

# Write results
write.table(results, file = "results_flow_metrics.csv", row.names = FALSE, sep = ",")

write.table(afdc.025, file = "DV_Flows/afdc_025.csv", row.names = FALSE, sep = ",")
write.table(afdc, file = "DV_Flows/afdc.csv", row.names = FALSE, sep = ",")
write.table(afdc.975, file = "DV_Flows/afdc_975.csv", row.names = FALSE, sep = ",")
