system.file(“F:/ENG330/Flow_Duration_curves.png”,package = “imager”)
The following code produces a flow duration curve for all daily data for Mooloolah river at Mooloolah. This data was accessed from the Queensland Government Water Monitoring Portal.
tfe=0.2 la<-c(0.1,1,10,100) ylim = range(pretty(c(0.1, zd[,3]))) fdc(zd[,3],ylab=varn[1],thr.shw=FALSE,col=“blue”,lty=1,pch=0,lwd=1,ylim=ylim,yat=la,new=TRUE) grid(nx= NULL, ny = NULL, lty = 1, lwd = 1, col=“gray60”,equilogs = TRUE)
qv<-quantile(zd[,3],1-tfe,na.rm=T) arrows(tfe,0.1,tfe,qv,col=“blue”) arrows(tfe,qv,0,qv,col=“blue”) points(tfe,qv,col=“blue”,cex=2) text(tfe,qv,paste(“A flow of”,ceiling(qv)," ML/day is exceeded“, tfe*100, “% of the time”),pos=4,offset=2,col=“blue”)
The following code produces a Flow Duration Curve with quantile distributions. #FLOW DURATION CURVE (QUANTILE DISTRIBUTION) tfe=0.2
stm1<-1 finm1<-2 if(stm1>finm1) { mon_ind<-c(seq(from=stm1,to=12),seq(from=1,to=finm1)) }else {mon_ind<-seq(from=stm1,to=finm1)}
maxy<-signif(max(zd[which.max(zd[,3]),]),2)
la<-c(0.1,1,10,100) fdc(zd[,3],ylab=varn[1],thr.shw=FALSE,col=“grey”,lty=3,lwd=0.2,pch=0,cex=0,ylim=ylim,yat=la,new=TRUE) grid(nx= NULL, ny = NULL, lty = 1, lwd = 1, col=“gray60”,equilogs = TRUE) zd1<-zd[month(index(zd)) %in% mon_ind,] fdc(zd1[,3],thr.shw=FALSE,lty=1,new=FALSE,col=“red”) qvw<-quantile(zd1[,3],1-tfe,na.rm=T) stm<-8 finm<-9 if(stm>finm) { mon_ind<-c(seq(from=stm,to=12),seq(from=1,to=finm)) }else {mon_ind<-seq(from=stm,to=finm)}
zd2<-zd[month(index(zd)) %in% mon_ind,] fdc(zd2[,3],thr.shw=FALSE,lty=1,new=FALSE,col=“brown”) qv<-quantile(zd[,3],1-tfe,na.rm=T)
qvd<-quantile(zd2[,3],1-tfe,na.rm=T) arrows(tfe,qvd,0,qvd,col=“brown”) text(0.01,(qvd-100),paste(“A flow of”,ceiling(qvd)," ML/day exceeded“, tfe*100, “% of the time in the months”, stm,“-”,finm),pos=4,col=“brown”)
arrows(tfe,0.1,tfe,qvw,col=“red”) arrows(tfe,qvw,0,qvw,col=“red”) text(tfe,qvw,paste(“A flow of”,ceiling(qvw)," ML/day is exceeded“, tfe*100, “% of the time in the months”, stm1,“-”,finm1),pos=4,col=“red”,offset=2)
The following is the code used to produce the Annual Maximum Series from the daily data for Mooloolah River at Mooloolah. #Annual Series Code library(hydroTSM) stn<- “141006A” fn<-paste(“~/ENG Assignment/cf141006A/141006A.csv”) ncol <- max(count.fields(fn, sep = “,”))
nrows<-length(count.fields(fn, sep = “,”)) nv<-as.integer((ncol-2)/6) hd<-read.csv(fn,header=FALSE,nrows=4,sep=“,”) varn<-character(0) cn<-2 for(i in 1:nv) {
varn[i]<-as.character(hd[3,cn]) cn=cn+6
}
nrws<-length(count.fields(fn, sep = “,”)) ddtemp<-read.csv(fn,header=FALSE,skip=4,sep=“,”,nrows=nrws-6) dd<-ddtemp[1:(nrow(ddtemp)),1:(ncol-1)] cnames<-NULL for(c in 1:nv){cnames<-c(cnames,paste(c(“Mean”,“Qual”,“Min”,“Qual”,“Max”,“Qual”),rep(c,6)))} names(dd)<-c(“Date”,cnames) stdate<- as.POSIXct(dd\(Date[1],format="%H:%M:%S %d/%m/%Y")-100*365.25 findate<-as.POSIXct(dd\)Date[nrow(dd)],format=“%H:%M:%S %d/%m/%Y”) date_ind<-seq(stdate,findate,by=“day”) zd<-zoo(dd[,c(2,4,6)],date_ind) nma<-4 zd\(yr<-as.numeric(as.yearmon(time(zd)) + nma/12) %/% 1 minyr<-as.integer(zd\)yr[1]) maxyr<-as.integer(zd$yr[nrow(zd)])
nyrs<-as.integer(maxyr-(minyr-1)) maxv<-numeric(nyrs) indd<-as.Date(nyrs) for(i in 1:nyrs){ temp<-subset(zd,zd$yr==(i+minyr-1)) maxv[i]<-max(temp[,1:3],na.rm=T) indd[i]<-as.Date(index(temp[which.max(temp[,3])])) }
indd2 <- as.POSIXct(paste(c(toString(indd[1]), “13:51:15”),collapse=" “), format=”%Y-%m-%d %H:%M:%S“)
for (i in 2:length(indd)){ indd2[i] <- as.POSIXct(paste(c(toString(indd[i]), “13:51:15”),collapse=" “), format=”%Y-%m-%d %H:%M:%S“) }
peaks <- zoo(maxv,indd2)
To construct an empirical distribution function using the Annual Maximum Series, firstly the Annual Maximum peaks must be found. The code for extracting these is shown below. #Plot of Peaks Code plot(zd[,3],col=“blue”,xlab=“”,ylab=varn[1]) points(peaks,col=“red”,pch=20)
abline(v=seq(index(zd[1]),index(zd[nrow(zd)]),by=“year”),col=“grey”)
A Density vs Discharge Histogram was produced showing PDF, normal and log-normal distibutions. The code used to construct the histogram is shown below: #Histogram Code library(lattice) wyv<-aggregate(zd, zd$yr, max,na.rm=TRUE) hv<-hist(coredata(wyv[,3]),breaks=15,plot=FALSE) hist(coredata(wyv[,3]),breaks=15,xlab=varn[1],freq=F,col=“blue”,main=“”)
lines(density(coredata(wyv[,3])),col=“red”) mv<-mean(coredata(wyv[,3])) sdv<-sd(coredata(wyv[,3])) curve(dnorm(x,mean=mv,sd=sdv),add=TRUE,col=“green”) curve(dlnorm(x,mean=log(mv),sd=log(sdv)),add=TRUE,col=“darkorange4”) yl<-max(hv\(density) xl<-hv\)mids[length(hv$mids)-5] legend(xl,yl,legend=c(“histogram”,“pdf”,“normal”,“log-normal”),col=c(“blue”,“red”,“green”,“darkorange4”),pch=c(15,NA,NA,NA),lty=c(0,1,1,1),cex=0.8)
rind<-rank(maxv) n<-length(maxv) pp<-(rind-0.4)/(n+0.2) plot(qnorm(pp),maxv,log=“y”,col=“blue”,pch=20,xlab=“z”,ylab=varn[1])
x <- seq(0, 100, 1) plot(x,dbinom(x,100,1/50),col=“red”,lwd=3,type=“h”, las=1, xlab=“Number of Occurrences of 50 Year Floods in 100 Years”, ylab=“Probability of Occurrence”, main=“Probability Distribution for number of Occurrences of 50 Year Floods in 100 Years”)
points(2,dbinom(2,100,1/50), pch=19, col =“blue”)
Shown below, is the code for the probability plot using the annual maximum series for estimating the AEP of an observed flood. A Distribution has been added to the plot with confidence limits to the distribution. The largest flood on record is a 1 in 25 year event with an annual peak of approximately 160m^3/s. #FFA AEP library(zoo) library(EnvStats) stn<- “141006A” fn<-paste(dir<-“~/ENG Assignment/cf141006A/141006A.csv”,stn,“.csv”,sep=“”) ncol <- max(count.fields(fn, sep = “,”)) nrows<-length(count.fields(fn, sep = “,”)) hd<-read.csv(fn,header=FALSE,nrows=4,sep=“,”) vn<-unique(unlist(hd[3,])) vn[vn==“”]<-NA vn<-na.omit(vn) nv<-length(vn)-1 cnames<-NULL varn<-character(0) cn<-2 for(i in 1:nv) {
varn[i]<-as.character(hd[3,cn]) if(hd[4,cn]==“Total”) { for(vni in 0:1){cnames<-c(cnames,paste(as.character(hd[4,(cn+vni)]),i,sep=“”))} cn<-cn+2} else { for(vni in 0:5){cnames<-c(cnames,paste(as.character(hd[4,(cn+vni)]),i,sep=“”))} cn<-cn+6}
}
nrws<-length(count.fields(fn, sep = “,”)) ddtemp<-read.csv(fn,header=FALSE,skip=4,sep=“,”,nrows=nrws-6) dd<-ddtemp[1:(nrow(ddtemp)),1:(ncol-1)] dd<-dd[complete.cases(dd),] colnames(dd)<-c(“Datev”,cnames) dd\(Datev<-strptime(dd\)Datev,format=“%H:%M:%S %d/%m/%Y”) dd\(Datev<-as.POSIXct(dd\)Datev)
zd<-zoo(dd[,c(2,4,6)],dd\(Datev) nma<-3 zd\)yr<-as.numeric(as.yearmon(time(zd)) + nma/12) %/% 1 minyr<-as.integer(zd\(yr[1]) maxyr<-as.integer(zd\)yr[nrow(zd)]) nyrs<-as.integer(maxyr-(minyr-1)) maxv<-numeric(nyrs) indd<-as.Date(nyrs) for(i in 1:nyrs) {temp<-subset(zd,zd$yr==(i+minyr-1)) maxv[i]<-max(temp[,3]) indd[i]<-as.Date(index(temp[which.max(temp[,3])])) } maxv<-na.omit(maxv) maxv<-maxv[1:length(maxv)] a<-0.4 n=length(maxv) pp<-(rank(maxv)-a)/(n+1-2*a)
T<-1/(1-pp)
Ttick<-c(1.001,1.01,1.1,1.5,2:20,seq(25,50,5),seq(60,100,10)) xtlab<-c(1.001,1.01,1.1,1.5,2,rep(NA,2),5, rep(NA,4),10,rep(NA,4),15,rep(NA,4),20,NA,30,rep(NA,3),50,rep(NA,4),100)
x<- -log(-log(1-1/T)) xtick<- -log(-log(1-1/Ttick)) xmin<-min(min(x),min(xtick)) xmax<-max(xtick) KTtick<–(sqrt(6)/pi)(0.5772+log(log(Ttick/(Ttick-1)))) QTtick<-mean(maxv)+KTticksd(maxv) se<-(sd(maxv)sqrt((1+1.14KTtick+1.1KTtick^2)))/sqrt(n) LB<-QTtick-qt(0.975,n-1)se UB<-QTtick+qt(0.975,n-1)*se max<-max(UB) Qmax<-max(QTtick) plot(x,maxv, ylab=expression(“Annual Peak (m^3/s)”), xaxt=“n”,xlab=“Annual Exceedance Probability (AEP) [1 in x years]”, ylim=c(0,Qmax), xlim=c(xmin,xmax), pch=21,bg=“red”, main=“Flood Frequency Analysis Mooloolah River at Mooloolah” ) par(cex=0.65) axis(1,at=xtick,labels=as.character(xtlab))
lines(xtick,QTtick,col=“red”,lty=1,lwd=2) lines(xtick,LB,col=“blue”,lty=1,lwd=1.5) lines(xtick,UB,col=“blue”,lty=1,lwd=1.5) abline(v=xtick,lty=3,col=“light gray”) abline(h=seq(50,floor(Qmax),50),lty=3,col=“light gray”)
ms<-sort(maxv,decreasing=TRUE) m<-order(ms,decreasing=TRUE) n<- length(maxv) a<-0.4 pp<-(m-a)/(n+1-2*a) z<-qnorm(1-pp)
qqv<-qqnorm(ms,datax=FALSE,plot=FALSE)
plot(qqv\(x,qqv\)y,log=“y”,pch=21,ylim=c(10,pretty(max(qqv$y))[2]*10),bg=“red”,axes=FALSE,ylab=varn[1],xlab=“Annual Exceedance Probability”)
Exceedances <- c(0.99, 0.98, 0.95, 0.90, 0.80, 0.50, 0.2, 0.05, 0.02, 0.01) Locations.x <- qnorm(Exceedances) axis(side=1, at=Locations.x, labels=(1-Exceedances), las=2)
box() abline(v=Locations.x, col=“gray”, lty=2)
digits = floor( log10( pretty(max(qqv$y))[2] ) ) + 1; Powers10 <- 0:digits axis(side=2, at=10^Powers10, labels=10^Powers10, las=1)
Steps <- 0:9 for(i in 1:digits) { abline(h=Steps,col=“Gray”,lty=2) Steps=Steps*10 }
nr<-1000000 #random samples from a log normal distribution to fit a line - not an ideal method dist_data<-rlnorm(nr,meanlog=mean(log(maxv),na.rm=T),sdlog=sd(log(maxv),na.rm=T)) qqd<-qqnorm(sort(dist_data,decreasing=TRUE),datax=FALSE,plot=FALSE)
xp<-c(qqd\(x[1],qqd\)x[nr]) yp<-c(qqd\(y[1],qqd\)y[nr]) points(xp,yp,type=“l”,col=“red”,cex=1)
nex<-length(Exceedances) civll_ul<-array(,dim=c(nex,2)) colnames(civll_ul)<-c(“LL”,“UL”) for(e in 1:nex) { civ<-eqlnorm(maxv, p =Exceedances[e] , method = “qmle”, ci = T, ci.method = “exact”, ci.type = “two-sided”, conf.level = 0.95, digits = 0) civll_ul[e,]<-civ\(interval\)limits
} points(Locations.x ,civll_ul[,1],type=“l”,col=“blue”) points(Locations.x ,civll_ul[,2],type=“l”,col=“blue”)
Below is the code to produce RFFE outputs by scraping the RFFE tool. This code can be easily manipulated to change return periods of AEP events. #RFFE library(RCurl) library(jsonlite) library(stringr) library(dplyr) library(ggplot2)
rffe.data <- postForm(“https://rffe.arr-software.org/”, catchment_name = “test1”, lato = “-26.686”, lono = “153.13”, latc = “-26.762”, lonc = “152.982”, area = “39” )
rffe.data <- as.character(rffe.data)
x <- stringr::str_match_all(rffe.data, ‘\[\{.*\}\]’ ) gauges.ffa.JSON <- x[[1]][1,] # information from surrounding gauges RFFE.res <-x[[1]][2,]
gauges.ffa.df <- jsonlite::fromJSON(gauges.ffa.JSON, flatten = TRUE)
gauges.ffa <- gauges.ffa.df %>% select(station_id, area, record.length = sflength, latc, lonc, lato, lono, bcf, i2_6h, i50_6h, mar, shape.factor = sf, flow_1pc = q1, flow_2pc = q2, flow_5pc = q5, flow_10pc = q10, flow_20pc = q20, flow_50pc = q50, flow_1pc_LCL = lower1, flow_2pc_LCL = lower2, flow_5pc_LCL = lower5, flow_10pc_LCL = lower10, flow_20pc_LCL = lower20, flow_50pc_LCL = lower50, flow_1pc_UCL = upper1, flow_2pc_UCL = upper2, flow_5pc_UCL = upper5, flow_10pc_UCL = upper10, flow_20pc_UCL = upper20, flow_50pc_UCL = upper50, region.name, region.number, region.version, statistics.mean = statistics.mean, statistics.mean_se, statistics.skew, statistics.skew_se, statistics.stdev, statistics.stdev_se)
stat.correlation <-t(sapply(gauges.ffa.df$statistics.correlations, unlist)) stat.correlation <- as.data.frame(stat.correlation) stat.correlation <- stat.correlation[ ,c(2,4,5)] names(stat.correlation) = c(‘cor_mean_sd’, ‘cor_mean_skew’, ‘cor_sd_skew’)
gauges.ffa <- cbind(gauges.ffa, stat.correlation)
RFFE.res RFFE.res <- str_replace_all(RFFE.res, “[{\[\]}]”, “”) RFFE.res <- str_replace_all(RFFE.res, c(“:” = “”, “‘“=”“,”aep" = “”, “lower_limit”= “”, “upper_limit” = “”, “flow” = “”, “\s+” = “”)) RFFE.res <- unlist(str_split(RFFE.res,’,‘)) RFFE.res <- as.data.frame(matrix(as.numeric(RFFE.res), ncol = 4, byrow = TRUE)) names(RFFE.res) <- c(’ARI’, ‘upper_limit’, ‘lower_limit’, ‘flow’)
MyTheme = theme_bw() + theme( panel.background = element_rect(fill=“gray98”), axis.title.x = element_text(colour=“grey20”, size=14, margin=margin(20,0,0,0)), axis.text.x = element_text(colour=“grey20”,size=12), axis.title.y = element_text(colour=“grey20”,size=14, margin = margin(0,20,0,0)), axis.text.y = element_text(colour=“grey20”,size=12), legend.title = element_text(colour=“grey20”,size=12), plot.margin = unit(c(0.5, 0.5, 1, 0.5), “cm”), # top, right, bottom, left panel.grid.minor = element_line(colour=“grey90”, size=0.2), panel.grid.major = element_line(colour=“grey90”, size=0.4))
my.x.labels = c(2, 20, 50,100) my.x.breaks <- qnorm(1/my.x.labels, lower.tail = FALSE)
my.y.breaks <- c(10, 100, 1000) my.y.breaks my.y.minor_breaks <- c(1:10, seq(10,100,10), seq(100, 1000, 100))
RFFE.res %>% mutate(AEP = 1/ARI) %>% mutate(z.score = qnorm(AEP, lower.tail = FALSE)) %>% ggplot() + geom_line(aes(x = z.score, y = flow), linetype = 1, colour = ‘blue’) + geom_point(aes(x = z.score, y = flow), colour = ‘black’, size = 3) + geom_line(aes(x = z.score, y = lower_limit), linetype = 2, colour = ‘blue’) + geom_line(aes(x = z.score, y = upper_limit), linetype = 2, colour = ‘blue’) + scale_x_continuous(name = ‘AEP 1 in Y years’, breaks = my.x.breaks, labels = my.x.labels ) + scale_y_log10(name = bquote(‘Peak discharge (m’ ^3 ‘s’ ^-1 ‘)’), limits = c(1,1000), breaks = my.y.breaks, labels = my.y.breaks, minor_breaks = my.y.minor_breaks) + MyTheme