Created
May 12, 2020 00:03
-
-
Save TonyLadson/b9cc317273d6b99cb6dc0f9b9505246c to your computer and use it in GitHub Desktop.
Code to reproduce the figures in the blog: ARR2019 – Areal Reduction Factors: interpolating between short and long duration ARFs
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # 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