# 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)