Non-Normal Flood Frequency Distributions

library(copula)
library(nsRFA)
library(lmomco)

Challenge 1: Generate Flood Discharge Data and Calculate the Return Interval

# read input file & assign data pairs matrix
load("D:/Users/Michelle Newcomer/Documents/Berkeley/Wohler Project/Data/USGS Gauging Stations/Guerneville Station 11467000/Guerne_PeakVolDur.RData")
#n<-1000
#df <- data.frame(year=sort(runif(n,1900,2014)), volume = 1:n + rnorm(n,mean=1000, sd=200),peak = 1:n + rnorm(n,mean=250,sd=100),duration = (1:n+rnorm(n,mean=25,sd=25))*0.1)

peak = duration[1:273,5] #m3/s
vol = duration[1:273,4]*1/1000000 #convert m3 to million cubic m
dur = duration[1:273,2] #days

Plot the peak discharge versus the volume. Plot the volume versus the duration

Challenge 2:

Is there a method we can use to model the volume based on the duration? Write the equation we can use to predict the volume from the duration. Write the equation we can use to predict the volume from the peak discharge. Write your equation here in the R markdown document using standard Latex equation input. Use this website to help you build the equation code: http://mathurl.com/

Challenge 3:

Plot a histogram of the volume, the duration, and the discharge. This will help us determine the normality of the data.

Challenge 4: Estimate mean volume, mean duration, and mean peak discharge.

#=============== fitting the gumbel distribution ========================
ML.parGUMB.v = ML_estimation(vol,dist="GUMBEL")
ML.parGUMB.d = ML_estimation(dur,dist="GUMBEL")
ML.parGUMB.p = ML_estimation(peak,dist="GUMBEL")
mean_v = ML.parGUMB.v[1]; alpha_v =ML.parGUMB.v[2];
mean_d = ML.parGUMB.d[1]; alpha_d =ML.parGUMB.d[2];
mean_p = ML.parGUMB.p[1]; alpha_p =ML.parGUMB.p[2];

Challenge 5: Calculate the exceedence probability, non-exceedence, and the return interval

#=============== determine the probability distribution =================
# probability of non-exceedance
F.vol = F.gumb (vol, mean_v, alpha_v)
F.dur = F.gumb (dur, mean_d, alpha_d)
F.peak = F.gumb (peak, mean_p, alpha_p)
# probability of exceedance & Return Period
P.vol = 1-F.vol ; T.vol = 1/P.vol;
P.dur = 1-F.dur ; T.dur = 1/P.dur;
P.peak = 1-F.peak ; T.peak = 1/P.peak; 

#=============== fitting the Kendall'tau ; theta =========================
F.vd = as.matrix(cbind(F.vol,F.dur))
F.pv = as.matrix(cbind(F.peak,F.vol))

gumbel.cop = gumbelCopula(dim = 2)
fit.tau <- fitCopula(gumbel.cop, F.vd, method="itau") #Find the theta parameter for each copula family directly
theta.vd = coef(fit.tau)
fit.tau <- fitCopula(gumbel.cop, F.pv, method="itau") #Find the theta parameter for each copula family directly
theta.pv = coef(fit.tau)

#============== compute the Copula distribution ===========================
gumbel.cop = gumbelCopula(theta.vd, dim = 2)
F.Cop.vd = pCopula(F.vd, gumbel.cop) # Probability of Non-exceedance
P.Cop.vd = 1-F.Cop.vd # Probability of exceedance
T.Cop.vd = 1/P.Cop.vd # Return Period
gumbel.cop = gumbelCopula(theta.pv, dim = 2)
F.Cop.pv = pCopula(F.pv, gumbel.cop) # Probability of Non-exceedance
P.Cop.pv = 1-F.Cop.pv # Probability of exceedance
T.Cop.pv = 1/P.Cop.pv # Return Period

Cop.vd = as.matrix(cbind(vol,dur,F.Cop.vd,P.Cop.vd,T.Cop.vd)) # result for observe data
Cop.pv = as.matrix(cbind(peak,vol,F.Cop.pv,P.Cop.pv,T.Cop.pv)) # result for observe 


#============== Create The Contour Plot of Bivariate Distribution
# Define the copula distribution
theta.vd

## param 
##   7.2

theta.pv

## param 
##  4.72

# Define the range of proability non-exceedance of each variables (vol,dur,peak)
n.sim =1000
fVol = seq(0.05,0.99998,length.out=n.sim) # Sim. probability of flood volume from their distribution
fDur = seq(0.05,0.99998,length.out=n.sim) # Sim. probability of flood durationfrom their distribution
fPeak = seq(0.05,0.99998,length.out=n.sim) # Sim. probability of flood peak fromtheir distribution

Vol = invF.gumb(fVol,mean_v,alpha_v) # Cal. Flood Volume
Dur = invF.gumb(fDur,mean_d,alpha_d) # Cal. Flood Volume
Peak = invF.gumb(fPeak,mean_p,alpha_p) # Cal. Flood Volume

gumbel.cop.vd = gumbelCopula(theta.vd, dim = 2)
gumbel.cop.pv = gumbelCopula(theta.pv, dim = 2)

pair.PV = cbind(expand.grid(fPeak,fVol)$Var1,expand.grid(fPeak,fVol)$Var2)
sim.PV = pCopula(pair.PV, gumbel.cop.pv)
sim.Cop.PV = matrix(sim.PV,n.sim,n.sim)
rnd.PV = rCopula(5000,gumbel.cop.pv)
rnd.PV.Vol = invF.gumb(rnd.PV[,2],mean_v,alpha_v)
rnd.PV.Peak = invF.gumb(rnd.PV[,1],mean_p,alpha_p)


pair.VD = cbind(expand.grid(fVol,fDur)$Var1,expand.grid(fPeak,fVol)$Var2)
sim.VD = pCopula(pair.VD, gumbel.cop.vd)
sim.Cop.VD = matrix(sim.VD,n.sim,n.sim)
rnd.VD = rCopula(5000,gumbel.cop.vd)
rnd.VD.Vol = invF.gumb(rnd.VD[,1],mean_v,alpha_v)
rnd.VD.Dur = invF.gumb(rnd.VD[,2],mean_d,alpha_d)