Take results from Flike and produce a nice looking frequency plot (see https://tonyladson.wordpress.com/2015/10/20/better-frequency-plots-from-web-based-flike/)

FlikePlot.R

##########################################################
#
#  Plotting Flike Output
#  see https://tonyladson.wordpress.com/2015/10/20/better-frequency-plots-from-web-based-flike/
#
############################################################



library(stringr)
library(R.utils)

# my.path <- c('As required')
# see example file at https://goo.gl/mDeHJp
my.file <- c('Tyers-test-graph.csv')

fname <- str_c(my.path, my.file, sep = '/')

# find the line numbers of those lines that contain the word 'Number'
lineNum.number <- grep('Number', readLines(fname)) # line numbers of heading

lineNum.eof <- as.vector(countLines(fname) ) # end of file


# Read in the gauged data and deviates

flike.data <- read.csv(fname, header = TRUE, nrows = lineNum.number[2] - 2 )


# # It looks like flike sets any zero flow values  to 0.001
# # we can delete those if necessary.
# 
# flike.data <- flike.data[-which(flike.data$Gauged_value == -3), ]


my.ari <- c(1.5, 2, 5, 10, 20, 50, 100)
my.aep <- 1/my.ari
my.z <- qnorm(1 - my.aep)

# Change suggested by Phil Pedruco, thank you
#set y limits to Q100 upper CI and min gauged data

# Read in the confidence limites and expected parameter quantile

flike.cl <- read.csv(fname, header = TRUE, skip = lineNum.number[2] - 1, nrows = lineNum.number[3] -  lineNum.number[2] - 1)

ylim_min <- signif(min(10^flike.data$Gauged_value), 3)
ylim_max <- max(flike.cl$Upper_90._probability_limit[11])
ylims <- c(ylim_min, 10^ylim_max)

par(oma = c(5,3,0,0))
plot(10^Gauged_value ~ Deviate, 
     data = flike.data,
     xaxt = 'n',
     yaxt = 'n',
     log = 'y',
     ylab = '',
     xlab = '',
     pch = 21,
     bg = 'grey',
     main = '',
     ylim = ylims
       )

axis(side = 1, at = my.z, my.ari)
mtext(text = 'ARI (years)', side = 1, line = 2)

my.label = str_c(round(100*my.aep), '%')

axis(side = 1, at = my.z, my.label, outer = TRUE)
mtext(text = 'AEP', side = 1, line = 7)

my.label = str_c(round(100*my.aep), '%')

axis(side = 1, at = my.z, my.label, outer = TRUE)
mtext(text = 'AEP', side = 1, line = 7)

my.labels <- prettyNum(axTicks(2), 
          scientific = FALSE, 
          big.mark = ',')

axis(side=2, at=axTicks(2), 
     labels=my.labels,
     las = 2)


mtext(side = 2, line = 4.5, text ='Instantaneous maximum flow (ML/d)')

abline(h = seq(0.1,1,0.1), lty = 3, col = 'grey' )
abline(h = seq(1,10,1), lty = 3, col = 'grey' )
abline(h = seq(10,100,10), lty = 3, col = 'grey' )
abline(h = seq(100,1000,10), lty = 3, col = 'grey' )
abline(h = seq(1000,10000,1000), lty = 3, col = 'grey' )
abline(h = seq(10000,100000,10000), lty = 3, col = 'grey' )

abline(v =my.z, lty = 3, col = 'grey' )



# Read in the confidence limites and expected parameter quantile
# plot on graph


flike.cl <- read.csv(fname, header = TRUE, skip = lineNum.number[2] - 1, nrows = lineNum.number[3] -  lineNum.number[2] - 1)

lines(10^Expected_par_quantile ~ Deviate, data = flike.cl) # Expected parameter quantile
lines(10^Lower_90._probabilty_limit ~ Deviate, data = flike.cl, lty = 2) # 
lines(10^Upper_90._probability_limit ~ Deviate, data = flike.cl, lty = 2) #

# Readin and plot the expected probability quantile

flike.exp.prob <- read.csv(fname, header = TRUE, skip = lineNum.number[3] - 1, nrows = lineNum.eof -  lineNum.number[3])

flike.exp.prob

# remove all rows where the Expected probably quantile value is zero

flike.exp.prob <- flike.exp.prob[flike.exp.prob$Expected_probability_quantitle > 0, ]
lines(10^Expected_probability_quantitle ~ Deviate, data = flike.exp.prob, col = 'red') # Expected parameter quantile

legend('bottomright', 
       legend = c('90% CL', 'Gauged', 'Expected param', 'Expected prob'), 
       lty = c(2, -1, 1, 1), 
       pch = c(-1, 21, -1, -1), 
       pt.bg = 'grey', 
       bty = 'n',
       col = c(1, 1, 1, 2),
       inset = 0.01,
       cex=0.9)