floodfreq=function(filed,monthmean=8,meanpref=c(0,0.5,1,1.5)) { #Be careful: the time series must include complete calendar years of data (start date maust be Jan 01 - End date must be Dec 31st) #Create the vector of dates startdate=scan(file=filed, what="character", n=1, skip=4) enddate=scan(file=filed, what="character", n=1, skip=6) delbis=scan(file=filed, what="character", n=1, skip=8) numonth=scan(file=filed, n=1, skip=10) floodseas=scan(file=filed, n=numonth, skip=12) flows=scan(file=filed, skip=14) pr1=seq(as.Date(startdate),as.Date(enddate),1) if(delbis==F) { #Remove 29s Februaries if delbis=T pr3=flows[strftime(pr1,"%m-%d")!="02-29"] pr1=pr1[strftime(pr1,"%m-%d")!="02-29"] } else { pr1=pr1[strftime(pr1,"%m-%d")!="02-29"] } #Compute number of years of the time series iiii=length(pr1[strftime(pr1,"%m-%d")=="01-01"]) #Start the loop to extract the annual maxima in the flood season and the mean flow in each year in the pre-flood season pr10=0 pr11=0 #Computation of startyear startyear=strftime(pr1[1],"%Y") #Extraction of the annual maxima in the flood season for (i in as.numeric(startyear):(as.numeric(startyear)+iiii-1)) { pr10=c(pr10,max(flows[as.numeric(strftime(pr1,"%Y"))==i & (as.numeric(strftime(pr1,"%m")) %in% floodseas)==T])) } pr10=pr10[2:length(pr10)] #Extraction of the mean flow in the pre-flood season for (i in as.numeric(startyear):(as.numeric(startyear)+iiii-1)) { pr11=c(pr11,mean(flows[as.numeric(strftime(pr1,"%Y"))==i & (as.numeric(strftime(pr1,"%m")) %in% monthmean)==T])) } pr11=pr11[2:length(pr11)] #Transformation of the two vectors to the normal space pr20=qqnorm(pr10,plot=F)$x pr21=qqnorm(pr11,plot=F)$x #Computation of the correlation between annual maxima and mean flow corcoef=cor(pr20,pr21) #Computation conditional mean and variance meancond=corcoef*meanpref stdcond=(1-corcoef^2)^0.5 #Plot annual maxima distribution alfa=1.283/sd(pr10) u=mean(pr10)-0.45*sd(pr10) quant=seq(min(flows)/2, max(flows)*1.5,length.out=1000) proba=exp(-exp(-alfa*(quant-u))) #Computation return period trit=1/(1-proba) #plot(quant,trit,log="y",type="l",lwd=2,col="red") pdf("output.pdf") plot(quant,-log(-log(proba)),type="l",lwd=2,col="red",axes=F,ylim=c(-log(-log(0.2)),-log(-log(0.995))),ylab="Return Period (years)", xlab=expression('Peak river flow (m'^3*'/s)')) axis(1) axis(2,at=-log(-log(c(0.2,0.5,0.8,0.9,0.95,0.98,0.99,0.995))),labels=c("","2","5","10","20","50","100","200")) segments(y0=c(-log(-log(c(0.2,0.5,0.8,0.9,0.95,0.98,0.99,0.995)))),x0=c(0,0,0,0,0,0,0,0),y1=c(-log(-log(c(0.2,0.5,0.8,0.9,0.95,0.98,0.99,0.995)))),x1=c(15000,15000,15000,15000,15000,15000,15000,15000),lty="dotted") segments(y0=rep(-log(-log(0.2)),15),x0=seq(0,15000,1000),y1=rep(-log(-log(0.995)),15),x1=seq(0,15000,1000),lty="dotted") #Generation conditional distributions qp=0 for (i in 1:length(meanpref)) { pr100=rnorm(10000,meancond[i],stdcond) pr200=approx(pr20,pr10,xout=pr100)$y pr200=pr200[is.na(pr200)==F] #Estimation of the Gumbel distribution parameters alfa=1.283/sd(pr200) u=mean(pr200)-0.45*sd(pr200) proba=exp(-exp(-alfa*(quant-u))) quant1=u-1/alfa*log(-log(0.995)) qp[i]=quant1 lines(quant,-log(-log(proba)),col="black") #lines(quant,trit,log="y",type="l",lwd=2,col="black") } dev.off() #return(corcoef) #return(qp) }