Areal reduction factors - some edge cases https://tonyladson.wordpress.com/2020/04/14/arr2019-areal-reduction-factors-some-edge-cases/

ARF_edge_cases.R

library(tidyverse)

devtools::source_url("https://gist.githubusercontent.com/TonyLadson/fc870cf7ebfe39ea3d1a812bcc53c8fb/raw/d8112631a92a32be749cabe334a22931c035711e/ARF2019.R?raw=TRUE")

#source(file.path('ARR2019_ARF', "ARF_2019.R"))
#source(file.path('ARR2019_ARF', "ARF_tests.R"))


# Functions and data ------------------------------------------------------



log_breaks = function(maj, radix=10) {
  function(x) {
    minx         = floor(min(logb(x,radix), na.rm=T)) - 1
    maxx         = ceiling(max(logb(x,radix), na.rm=T)) + 1
    n_major      = maxx - minx + 1
    major_breaks = seq(minx, maxx, by=1)
    if (maj) {
      breaks = major_breaks
    } else {
      steps = logb(1:(radix-1),radix)
      breaks = rep(steps, times=n_major) +
        rep(major_breaks, each=radix-1)
    }
    radix^breaks
  }
  
  
  
  
}
scale_x_log_eng = function(..., radix=10) {
  scale_x_continuous(...,
                     trans=scales::log_trans(radix),
                     breaks=log_breaks(TRUE, radix),
                     minor_breaks=log_breaks(FALSE, radix))
}
scale_y_log_eng = function(..., radix=10) {
  scale_y_continuous(...,
                     trans=scales::log_trans(radix),
                     breaks=log_breaks(TRUE, radix),
                     minor_breaks=log_breaks(FALSE, radix))
}

# panel labels for ggplot facet_wrap

aep.labs <- str_c('AEP ', c(0.5, 1, 2, 5, 10, 20, 50),'%')
names(aep.labs) <- c(0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5)


params <-
  structure(list(`East Coast North` =     c(0.327,  0.241, 0.448, 0.36,  0.00096,   0.48,  -0.21,   0.012,   -0.0013), 
                 `Semi-arid Inland QLD` = c(0.159,  0.283, 0.25,  0.308, 7.3e-07,   1,      0.039,  0,        0), 
                 Tasmania =               c(0.0605, 0.347, 0.2,   0.283, 0.00076,   0.347,  0.0877, 0.012,   -0.00033), 
                 `SW WA` =                c(0.183,  0.259, 0.271, 0.33,  3.845e-06, 0.41,   0.55,   0.00817, -0.00045), 
                 `Central NSW` =          c(0.265,  0.241, 0.505, 0.321, 0.00056,   0.414, -0.021,  0.015,   -0.00033), 
                 `SE Coast` =             c(0.06,   0.361, 0,     0.317, 8.11e-05,  0.651,  0,      0,        0), 
                 `Southern Semi-arid` =   c(0.254,  0.247, 0.403, 0.351, 0.0013,    0.302,  0.058,  0,        0), 
                 `Southern Temperate` =   c(0.158,  0.276, 0.372, 0.315, 0.000141,  0.41,   0.15,   0.01,    -0.0027), 
                 `Northern Coastal` =     c(0.326,  0.223, 0.442, 0.323, 0.0013,    0.58,  -0.374,  0.013,   -0.0015), 
                 `Inland Arid` =          c(0.297,  0.234, 0.449, 0.344, 0.00142,    0.216, 0.129,  0,        0)), 
            class = "data.frame", row.names = c("a", "b", "c", "d", "e", "f", "g", "h", "i"))




# Plots -------------------------------------------------------------------




## Creating a plot of ARF agaist area

area <- 10^seq(0,5,length.out = 500)
duration <- c(30, 60, 120, 180, 360, 540, 720, 1080, 1440, 2880, 4320, 5760, 7200, 8640, 10080)
AEP <- c(0.5)



xdf <- expand_grid(area = area, duration = duration) %>% arrange(duration)


p <- xdf %>% 
  rowwise() %>% 
  mutate(ARF_calc = ARF(area, duration, aep = AEP, region = 'Tasmania', params = params))  %>% 
  ggplot(aes(area, ARF_calc, colour = factor(duration))) +
  geom_line(size = 0.2) + 
  scale_x_log_eng(name = bquote('Area'~km^2), labels = scales::label_comma(accuracy = 1)) +
  scale_y_continuous(name = 'Areal Reduction Factor') +
  scale_colour_discrete(name = 'Duration\n(min)') +
  guides(colour = guide_legend(reverse=T)) +
  labs(title = 'Tasmania', subtitle = '50% AEP') +
  theme_grey(base_size = 5) +
  theme(legend.key.size = unit(2, 'mm'))

p

# We can ignore the warnings because they arise from attempting calculate ARFs in areas were equations aren't appropriate

ggsave(file.path('figure','ARF-area.png'), plot = p, width = 4, height = 3)


area <- 10^seq(0,5,length.out = 500)
duration <- c(30, 60, 120, 180, 360, 540, 720, 1080, 1440, 2880, 4320, 5760, 7200, 8640, 10080)
AEP <- c(0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5)
xdf <- expand_grid(area = area, duration = duration, AEP)

aep.labs <- str_c('AEP ', 100*AEP, '%')
names(aep.labs) <- AEP




p <- xdf %>% 
  rowwise() %>% 
  mutate(ARF_calc = ARF(area, duration, aep = AEP, region = 'Tasmania', params = params))  %>% 
  ggplot(aes(area, ARF_calc, colour = rev(factor(duration)))) +
  geom_line(size = 0.2) + 
  scale_x_log_eng(name = bquote('Area'~km^2), labels = scales::label_comma(accuracy = 1)) +
  scale_y_continuous(name = 'Areal Reduction Factor') +
  scale_colour_discrete(name = 'Duration\n(min)') +
  labs(title = 'Tasmania') +
  theme_grey(base_size = 5) +
  theme(legend.key.size = unit(2, 'mm')) +
  facet_wrap(~AEP, labeller = labeller(AEP = aep.labs))

p

ggsave(file.path('figure','ARF-area-facet_AEP.png'), plot = p, width = 4, height = 3)


## Plots against duration


AEP <- c(0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5)
area <- c(5, 10, 100, 1000, 10000, 20000)
duration <- 10^seq(0,4,length.out = 500)

xdf <- expand_grid(area = area, AEP = AEP, duration = duration)

p <- xdf %>% 
  rowwise() %>% 
  mutate(ARF_calc = ARF(area, duration, aep = AEP, region = 'Tasmania', params = params))  %>% 
  ggplot(aes(x = duration, y = ARF_calc, colour = factor(area))) +
  geom_line(size = 0.2) +
  geom_vline(xintercept = 720, linetype = 'dotted', size = 0.2, alpha = 0.5) +
  geom_vline(xintercept = 1440, linetype = 'dotted', size = 0.2, alpha = 0.5) +
  scale_x_log_eng(name = 'Duration (min)', labels = scales::label_comma(accuracy = 1)) +
  scale_y_continuous(name = 'Areal Reduction Factor') +
  scale_colour_hue(name = bquote('Area ('~km^2*')'))+
  labs(title = 'Tasmania') +
  theme_grey(base_size = 5) +
  theme(legend.position = c(0.6, 0.15)) + 
  theme(legend.key.size = unit(2, 'mm')) +
  facet_wrap(~AEP, labeller = labeller(AEP = aep.labs))

p

ggsave(file.path('figure','ARF-duration-facet_AEP.png'), plot = p, width = 4, height = 3)



##########################

# Downhill segment in Figure 3

# For some scenarios, the 1440 min long duration ARF is smaller than the 720 min short duration ARF
# This means the ARF v duration graph is not monotonic and instead slopes downhill for a short period

ARF_short(26, 720, aep = 0.0005)
#[1] 0.9377527

ARF_long(area = 26, duration = 1440, aep = 0.0005, region = 'Tasmania', params = params)
#[1] 0.9322746



########################

# Plot short and long duration ARF curves for Tasmania for a catchment area of 26 km2

duration <- 10^seq(0,4,length.out = 500)
area <- c(26)

xdf <- expand_grid(duration, area) %>% 
  rowwise() %>% 
  mutate(ARF_s = ARF_short(area, duration, aep = 0.005)) %>% 
  mutate(ARF_l = ARF_long(area, duration, aep = 0.005, region = 'Tasmania', params = params)) 
  



xdf %>% 
ggplot(aes(x = duration, y = ARF_s)) +
  geom_line() +
  geom_line(data = xdf, aes(x = duration, y = ARF_l), colour = 'orange') +
  scale_x_continuous(trans = 'log10') +
  geom_vline(xintercept = 720) +
  geom_vline(xintercept = 1440) +
  scale_y_continuous(limits = c(0.8, NA))
           
# the short duration curve exceed the long duration curve beyond 720 min