Distribution of burst initial loss https://tonyladson.wordpress.com/2019/08/13/the-distribution-of-burst-initial-loss/

Burst_IL.R

# Burst initial loss = storm initial loss - preburst rainfall


# Storm initial loss ------------------------------------------------------


# emprical distribution of initial loss
# From ARR Table Table 5.3.13
# http://book.arr.org.au.s3-website-ap-southeast-2.amazonaws.com

loss_std <- structure(list(Percentile = c(100, 90, 80, 70, 60, 50, 40, 30, 
                                          20, 10, 0), IL = c(0.14, 0.39, 0.53, 0.68, 0.85, 1, 1.2, 1.4, 
                                                             1.71, 2.26, 3.19), CL = c(0.15, 0.35, 0.48, 0.61, 0.79, 1, 1.24, 
                                                                                       1.5, 1.88, 2.48, 3.85), prob = c(0, 0.1, 0.2, 0.3, 0.4, 0.5, 
                                                                                                                        0.6, 0.7, 0.8, 0.9, 1)), class = c("tbl_df", "tbl", "data.frame"
# re-formatting and relabelling                                                                                                                        ), row.names = c(NA, -11L))

loss_std <- loss_std %>% 
  mutate(prob = (100-Percentile)/100) %>% 
  arrange(prob)

# Function to interpolate the table to get a loss for a given percentile
# (percentile need to be defined the 'correct' way i.e. high percentile = high loss
# This is consistent with way percentiles are defined for pre-burst

IL_f <- approxfun(x = 100-loss_std$Percentile, y= loss_std$IL)

# The median ILs value for Toomuc Creek 
# From the data hub 
# -38.064520
# 145.463277
# 25 mm


# Generate 10000 random values of initial loss (given a median initial loss of 25 mm)
ILs <-  25*IL_f(100*runif(10000))



# Pre-burst rainfall ------------------------------------------------------

# Pre-burst rainfall for Toomuc Creek
# 2 hour 1% AEP event
# from the data hub https://data.arr-software.org
# Lat = -38.064520 Lon = 145.463277


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
# The data hub preburst values only have the following percentiles 10, 25, 50, 75, 90
# To randomly generate values accross the full range of percentiles, we need to 
# Extrapolate to 100 percent, and include a value for the zeroth percentile

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

# function to generate a preburst value for a given percentile (by interpolating the data hub values)
PB_f <- approxfun(x = pre_burst_ext$percentile, y= pre_burst_ext$depth)

set.seed(2001)

# Random initial loss
PB <- PB_f(runif(10000))


# Burst inital loss -------------------------------------------------------

# Burst initial loss = Storm initial loss - Pre-burst rainfall

hist(ILs)
median(ILs)

hist(PB)
median(PB)

ILb = ILs - PB

median(ILb)
#17.6

hist(ILb, breaks = 'FD')

hist(ILb, plot = 'FALSE', breaks = 'FD' )



ILb <- tibble(ILb = ILb)

p <- ILb %>% 
  ggplot(aes(x = ILb)) + 
  geom_histogram(binwidth = 2, fill = 'blue', colour = 'black') +
  scale_x_continuous(name = 'Burst initial loss (mm)', breaks = seq(-60,80,10)) +
  scale_y_continuous(name = 'Count') +
  theme_grey(base_size = 7)

p

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

ILb %>% 
  mutate(ILb = if_else(ILb < 0, 0, ILb)) %>% View
  summarise(median(ILb))
  
  # 17.6


# with vales less than zero set to zero
p <- ILb %>% 
  mutate(ILb = if_else(ILb < 0, 0, ILb)) %>% 
  ggplot(aes(x = ILb)) + 
  geom_histogram(binwidth = 2, fill = 'blue', colour = 'black') +
  scale_x_continuous(name = 'Burst initial loss (mm)', breaks = seq(-60,80,10)) +
  scale_y_continuous(name = 'Count') +
  theme_grey(base_size = 7)

p

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