Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment