# GumbelPlot.R # # Code for plotting annual peak flow series on # extreme-value (Gumbel) paper. # # This code illustrates how to customize graph axes, and # also how to use superscripts in axis labels. # # RDM 2014 Feb 6 ############################################################### rm(list = ls()) # Specify Greata Creek peak flows for 1971 - 1980 Q = c(1.23,2.37,0.085,1.69,1.2,0.898,0.176,0.96,0.212,0.266) graphlab = "Greata Creek 1971-1980" # Generate plotting positions n = length(Q) r = n + 1 - rank(Q) # highest Q has rank r = 1 T = (n + 1)/r # Set up x axis tick positions and labels Ttick = c(1.001,1.01,1.1,1.5,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,25,30,35,40,45,50,60,70,80,90,100) xtlab = c(1.001,1.01,1.1,1.5,2,NA,NA,5,NA,NA,NA,NA,10,NA,NA,NA,NA,15,NA,NA,NA,NA,20,NA,30,NA,NA,NA,50,NA,NA,NA,NA,100) y = -log(-log(1 - 1/T)) ytick = -log(-log(1 - 1/Ttick)) xmin = min(min(y),min(ytick)) xmax = max(ytick) # Fit a line by method of moments, along with 95% confidence intervals KTtick = -(sqrt(6)/pi)*(0.5772 + log(log(Ttick/(Ttick-1)))) QTtick = mean(Q) + KTtick*sd(Q) nQ = length(Q) se = (sd(Q)*sqrt((1+1.14*KTtick + 1.1*KTtick^2)))/sqrt(nQ) LB = QTtick - qt(0.975, nQ - 1)*se UB = QTtick + qt(0.975, nQ - 1)*se max = max(UB) Qmax = max(QTtick) # Plot peak flow series with Gumbel axis plot(y, Q, ylab = expression( "Peak Flow ("*m^3*s^{-1}*")" ) , xaxt = "n", xlab = "T (yr)", ylim = c(0, Qmax), xlim = c(xmin, xmax), pch = 21, bg = "red", main = paste( "Annual Peak Flows,",graphlab ) ) par(cex = 0.65) axis(1, at = ytick, labels = as.character(xtlab)) # Add fitted line and confidence limits lines(ytick, QTtick, col = "black") lines(ytick, LB, col = "black", lty = 3) lines(ytick, UB, col = "black", lty = 3) # Draw grid lines abline(v = ytick, lty = 3) abline(h = seq(0.5, floor(Qmax), 0.5), lty = 3) par(cex = 1)