Skip to content

Instantly share code, notes, and snippets.

@TonyLadson

TonyLadson/ARF_interp.R

Created May 12, 2020 00:03
Show Gist options
  • Save TonyLadson/b9cc317273d6b99cb6dc0f9b9505246c to your computer and use it in GitHub Desktop.
Save TonyLadson/b9cc317273d6b99cb6dc0f9b9505246c to your computer and use it in GitHub Desktop.
Download ZIP
Code to reproduce the figures in the blog: ARR2019 – Areal Reduction Factors: interpolating between short and long duration ARFs
Raw
ARF_interp.R
# Load the functions we need
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")) # Check that we pass tests
# Functions and constants ------------------------------------------------------
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))
}
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
# Figure 1 from blog
xdf <- tribble(~area_min, ~area_max, ~dur_min, ~dur_max, ~type,
0, 1000, 0, 12, 'short',
0, 30000, 24, 168, 'long',
0, 30000, 12, 24, 'interpolation')
p <- xdf %>%
ggplot(aes(xmin = area_min,
xmax = area_max,
ymin = dur_min,
ymax = dur_max)) +
geom_rect(aes(fill = type)) +
scale_x_continuous(name = bquote('Catchment Area'~km^2), breaks = c(0, 1000, 30000), labels = c('', '1000', '30,000')) +
scale_y_continuous(name = 'Duration (hours)', breaks = c(0, 12, 24, 168)) +
scale_fill_brewer(palette = "Blues") +
theme_classic(base_size = 7) +
annotate(geom = 'text', x = 1200, y = 6, label = 'Short duration', hjust = 'left', size = 2) +
annotate(geom = 'text', x = 15000, y = 18, label = 'Interpoation', size = 2) +
annotate(geom = 'text', x= 15000, y = 96, label = 'Long duration', size = 2)+
guides(fill = 'none')
p
ggsave(filename = file.path('figure','interp.png'), plot = p, width = 4, height = 3)
# Figure 2 from blog
AEP <- c(0.005)
area <- c(30000)
duration <- 10^seq(0,4,length.out = 500)
duration <- append(duration, c(719, 720, 722, 725, 727, 1440)) %>% sort()
ARF_df <- expand_grid(area = area, AEP = AEP, duration = duration)
ARF_long_df <-
ARF_df %>%
rowwise() %>%
mutate(ARF_calc = ARF_long(area, duration, aep = AEP, region = 'Tasmania', params = params)) %>%
mutate(ARF_calc = if_else(duration < 720, NA_real_, ARF_calc))
p <- ARF_df %>%
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(data = ARF_long_df, aes(x = duration, y = ARF_calc), linetype = 'dashed', colour = scales::hue_pal()(3)[3], size = 0.2) +
geom_line(size = 0.2, colour = scales::hue_pal()(3)[3]) + # size = 0.2
annotate(geom = 'text', x = 1420, y = 0, angle = 90, label = '24 hours', size = 2, hjust = 0) +
annotate(geom = 'text', x = 700, y = 0, angle = 90, label = '12 hours', size = 2, hjust = 0) +
annotate(geom = 'text', x = 5000, y = 0.93, angle = 0, label = bquote('1000'~km^2), colour = scales::hue_pal()(3)[1]) +
annotate(geom = 'text', x = 4000, y = 0.89, angle = 0, label = bquote('1100'~km^2), colour = scales::hue_pal()(3)[2]) +
annotate(geom = 'text', x = 4000, y = 0.65, angle = 0, label = bquote('30,000'~km^2), colour = scales::hue_pal()(3)[3]) +
annotate(geom = 'point', x = 720, y = ARF_short(area = 30000, duration = 720, aep = 0.005), colour = 'black', size = 1) + # colour = scales::hue_pal()(3)[3]
annotate(geom = 'point', x = 720, y = ARF_long(area = 30000, duration = 720, aep = 0.005, params = params,region = 'Tasmania'), colour = 'black', size = 1) +
annotate(geom = 'point', x = 1440, y = ARF_long(area = 30000, duration = 1440, aep = 0.005, params = params,region = 'Tasmania'), colour = 'black', size = 1) +
#annotate(geom = "segment", x = 1440, y=0, xend = 1440, yend = 0.01) +
annotate(geom = 'text', x = 1475, y = 0.55, label = 'Long duration, 24 hours', hjust = 'left', size = 2) +
annotate(geom = 'segment', x = 1470, y =0.55, xend = 1440, yend = 0.625, size = 0.2 ) +
annotate(geom = 'text', x = 805, y = 0.45, label = 'Short duration, 12 hours', hjust = 'left', size = 2) +
annotate(geom = 'segment', x = 800, y =0.45, xend = 720, yend = ARF_short(area = 30000, duration = 720, aep = 0.005), size = 0.2 ) +
annotate(geom = 'text', x = 815, y = 0.640, label = 'Long duration, 12 hours', hjust = 'left', size = 2) +
annotate(geom = 'segment', x = 810, y =0.640, xend = 720, yend = ARF_long(area = 30000, duration = 720, aep = 0.005, params = params, region = 'Tasmania'), 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)') +
scale_x_continuous(name = 'Duration (min)') +
scale_y_continuous(name = 'Areal Reduction Factor') +
scale_colour_hue(name = bquote('Area ('~km^2*')'))+
labs(title = 'Tasmania',
subtitle = bquote("AEP = 0.5%, Area = 30,000"~km^2)) +
theme_grey(base_size = 6) + #5
theme(legend.position = c(0.6, 0.3)) +
coord_cartesian(xlim = c(500, 2000), ylim = c(0,1)) +
guides(colour = 'none')
p
ggsave(filename = file.path('figure','ARF_interp.png'), plot = p, width = 4, height = 3)
# Figure 3 from blog
AEP <- c(0.005)
area <- c(1000, 1100, 30000)
duration <- 10^seq(0,4,length.out = 500)
duration <- append(duration, c(719, 720, 722, 725, 727, 1440)) %>% sort()
ARF_df <- expand_grid(area = area, AEP = AEP, duration = duration)
area <- c(1100, 30000)
duration <- 10^seq(0,4,length.out = 500)
duration <- append(duration, c(719, 720, 722, 725, 727, 1440)) %>% sort() # add a few extra points so there is a clean plot near 720 min
ARF_long_df <-
expand_grid(area = area, AEP = AEP, duration = duration) %>%
rowwise() %>%
mutate(ARF_calc = ARF_long(area, duration, aep = AEP, region = 'Tasmania', params = params)) %>%
mutate(ARF_calc = if_else(duration < 720, NA_real_, ARF_calc))
ARF_short(area = 30000, duration = 720, aep = 0.005) #[1] 0.5151772
my_line_size <- 0.2
my_text_size <- 1.8
my_point_size <- 1
# log scale
p <- ARF_df %>%
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(data = ARF_long_df, aes(x = duration, y = ARF_calc), linetype = 'dashed', size = my_line_size ) +
geom_line(size = my_line_size) + # size = 0.2
annotate(geom = 'text', x = 1480, y = 0, angle = 90, label = '24 hours', size = my_text_size, hjust = 'left') +
annotate(geom = 'text', x = 700, y = 0, angle = 90, label = '12 hours', size = my_text_size, hjust = 'left') +
annotate(geom = 'text', x = 5000, y = 0.93, angle = 0, label = bquote('1000'~km^2), colour = scales::hue_pal()(3)[1], size = my_text_size) +
annotate(geom = 'text', x = 4000, y = 0.88, angle = 0, label = bquote('1100'~km^2), colour = scales::hue_pal()(3)[2], size = my_text_size) +
annotate(geom = 'text', x = 4000, y = 0.65, angle = 0, label = bquote('30,000'~km^2), colour = scales::hue_pal()(3)[3], size = my_text_size) +
annotate(geom = 'point', x = 720, y = ARF_short(area = 30000, duration = 720, aep = 0.005), colour = scales::hue_pal()(3)[3], size = my_point_size) +
annotate(geom = 'point', x = 720, y = ARF_long(area = 30000, duration = 720, aep = 0.005, params = params,region = 'Tasmania'), colour = scales::hue_pal()(3)[3], size = my_point_size) +
annotate(geom = 'point', x = 1440, y = ARF_long(area = 30000, duration = 1440, aep = 0.005, params = params,region = 'Tasmania'), colour = scales::hue_pal()(3)[3], size = my_point_size) +
annotate(geom = 'point', x = 720, y = ARF_short(area = 1100, duration = 720, aep = 0.005), colour = scales::hue_pal()(3)[2], size = my_point_size) +
annotate(geom = 'point', x = 720, y = ARF_long(area = 1100, duration = 720, aep = 0.005, region = 'Tasmania', params = params), colour = scales::hue_pal()(3)[2], size = my_point_size) +
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)') +
#scale_x_continuous(name = 'Duration (min)') +
scale_y_continuous(name = 'Areal Reduction Factor') +
scale_colour_hue(name = bquote('Area ('~km^2*')'))+
labs(title = 'Tasmania',
subtitle = "AEP = 0.5%") +
theme_grey(base_size = 6) + #5
theme(legend.position = c(0.6, 0.3)) +
coord_cartesian(xlim = c(500, 10000), ylim = c(0,1)) +
guides(colour = 'none') +
NULL
ggsave(filename = file.path('figure','ARF_interp_3.png'), plot = p, width = 4, height = 3)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment