Distribution of pre-burst rainfall. See the blog https://tonyladson.wordpress.com/2019/08/06/the-distribution-of-pre-burst-rainfall/

Pre-burst.R

# Data from the ARR data hub for Toomuc Creek

pre_burst <- tribble(~perc, ~depth, ~ ratio,
                     10,0,0,
                     25,0,0,
                     50,1.5, 0.029,
                     75,13.3, 0.256,
                     90,38.4, 0.739)


# Extrapolate to 100

(depth_100 <- Hmisc::approxExtrap(x = pre_burst$perc, y = pre_burst$depth, xout = 100)$y)
#55.13
(ratio_100 <- Hmisc::approxExtrap(x = pre_burst$perc, y = pre_burst$ratio, xout = 100)$y)
#1.061

pre_burst <- pre_burst %>% 
  mutate(percentile = perc/100)

pre_burst

pre_burst_ext <- pre_burst %>% 
  add_row(perc = 100, depth = depth_100, ratio = ratio_100, percentile = 1) %>% 
  add_row(perc = 0, depth = 0, ratio = 0, percentile = 0, .before = 1)


pre_burst
pre_burst_ext

plot(ratio ~ percentile, data= pre_burst, type = 'b')
plot(depth ~ percentile, data= pre_burst, type = 'b')



# plot the distribution of pre-burst values
# As a cumulative distribution function
p <- pre_burst_ext %>% 
  ggplot(aes(depth, percentile )) +
  geom_point() +
  geom_line(linetype = 'dashed', colour = 'blue') +
  geom_line(data = pre_burst, colour = 'blue') +
  scale_x_continuous(name = 'Preburst depth (mm)', breaks = seq(0,60,10)) +
  scale_y_continuous(name = 'Percentile') +
  geom_hline(yintercept = 1.00) +
  theme_grey(base_size = 7)

p
ggsave(filename = file.path('../figures', 'depth_cdf.png'), plot = p, width = 4, height = 3)



# -------------------------------------------------------------------------

# Randomly generate pre-burst depths


PB_f <- approxfun(x = pre_burst_ext$percentile, y= pre_burst_ext$depth)

set.seed(2001)

# Random initial loss
PB_f(runif(10))

pbs <- tibble(pb = PB_f(runif(10000)))

median(pbs$pb)
# 1.5

mean(pbs$pb)
# 10.66

mean(pbs$pb > 20)
# 0.2139

# plot
hist(pbs$pb, breaks = 'FD')
# Check the FD breaks

hist(pbs$pb, breaks = 'FD', plot = 'FALSE')



# Create graph for blog
p <- pbs %>% 
  ggplot(aes(pb)) +
  geom_histogram(binwidth = 1, fill = 'blue', colour = 'black') +
  scale_x_continuous(name = 'Pre-burst rainfall depth (mm)', breaks = seq(0,60,5)) +
  scale_y_continuous(name = 'Count') +
  theme_grey(base_size = 7)

p  
ggsave(file.path('../figures', 'Pre_burst_hist.png'), p, width = 4, height = 3)


# -------------------------------------------------------------------------

# Randomly generate pre-burst ratios

PB_ratio_f <- approxfun(x = pre_burst_ext$percentile, y= pre_burst_ext$ratio)

set.seed(2001)

# Random initial loss
PB_ratio_f(runif(10))

pb_ratios <- tibble(ratio = PB_ratio_f(runif(10000)))

median(pb_ratios$ratio)
#[1] 0.029

mean(pb_ratios$ratio)
# 1] 0.2051929

# plot
hist(pb_ratios$ratio, breaks = 'FD')
# Check the FD breaks

hist(pb_ratios$ratio, breaks = 'FD', plot = 'FALSE')



# Create graph for blog
p <- pb_ratios %>% 
  ggplot(aes(ratio)) +
  geom_histogram(binwidth = 0.02, fill = 'blue', colour = 'black') +
  scale_x_continuous(name = 'Pre-burst to burst ratio', breaks = seq(0, 2, 0.1)) +
  scale_y_continuous(name = 'Count') +
  theme_grey(base_size = 7)

p  
ggsave(file.path('../figures', 'Pre_burst_ratio_hist.png'), p, width = 4, height = 3)