Recipe for creating a flood frequency plot see https://tonyladson.wordpress.com/2017/06/09/flood-frequency-plots-using-ggplot/

FFA_plot.R

library(scales)
library(tidyverse)  
library(stringr)  

# The tribble function provides a good way of entering small data sets
# flows are in ML/d

ff_comparison <- tribble(
  ~aep, ~broken_casey, ~seven_kialla,
  0.5, 6900, NA,
  0.2, 23450, 21400,
  0.1, 37800, 33400,
  0.05, 53600, 46300,
  0.02, 76000, 64100,
  0.01, 95000, 77700
)

# Add a column of normal quantiles that correspond to 
# AEP values and a column of 1 in X AEP values

ff_comparison <- ff_comparison %>% 
  mutate(z = qnorm(1 - aep),
         aep_1inX = 1/aep) 


# Copy this dataframe so we can use it for labelling
my_labels <-  ff_comparison
my_labels


# Change the format of the dataframe to make plotting easier
ff_comparison  <- ff_comparison %>% 
  gather(key = 'source', value = peak_flow, -aep, -z, -aep_1inX) # one observation per row

# Plot
ff_comparison %>% 
  mutate(peak_flow = peak_flow/86.4) %>% # convert to cumec
  ggplot(aes(x = z, y = peak_flow, colour = source)) + 
  geom_point() +
  geom_line() +
  scale_y_continuous(name = expression(Peak~flow~(m^{3}~s^{-1})), labels = scales::comma,
                     breaks = seq(0, 1000, 200)) +
  #    scale_y_log10(name = expression(Peak~flow~(m^{3}~s^{-1})), labels = comma,  # log scale
  #                  breaks = c(seq(10000, 100000, 10000))) +
  scale_x_continuous(name = 'AEP', breaks = my_labels$z, labels = str_c(100* my_labels$aep,'%'),
                     sec.axis = dup_axis(name = 'AEP (1 in X years)', labels = my_labels$aep_1inX)) +
  scale_colour_manual(name = 'Waterway', labels = c('Broken River at Casey Weir', 'Sevens Cks at Kialla West'),
                     values = c('light blue', 'dark blue')) +
  theme_grey(base_size = 12) +
  theme(legend.position = c(0.8, 0.2))

#  Can also use the Probit transformation to achieve much the same result
# Note trans = 'probit' has been added to scale_x_continuous and the x aesthetic is (1 - aep)
# We need the 1 - aep for the exceedance probability

ff_comparison %>% 
  mutate(peak_flow = peak_flow/86.4) %>% # convert to cumec
  ggplot(aes(x = 1-aep, y = peak_flow, colour = source)) + 
  geom_point() +
  geom_line() +
  scale_y_continuous(name = expression(Peak~flow~(m^{3}~s^{-1})), labels = scales::comma,
                     breaks = seq(0, 1000, 200)) +
  #    scale_y_log10(name = expression(Peak~flow~(m^{3}~s^{-1})), labels = comma,  # log scale
  #                  breaks = c(seq(10000, 100000, 10000))) +
  scale_x_continuous(name = 'AEP', 
                     trans = 'probit',
                     breaks = 1-my_labels$aep, 
                     labels = str_c(100* my_labels$aep,'%'),
                     sec.axis = dup_axis(name = 'AEP (1 in X years)', 
                                         labels = my_labels$aep_1inX)) +
  scale_colour_manual(name = 'Waterway', labels = c('Broken River at Casey Weir', 'Sevens Cks at Kialla West'),
                      values = c('light blue', 'dark blue')) +
  theme_grey(base_size = 12) +
  theme(legend.position = c(0.8, 0.2))