Skip to content

Instantly share code, notes, and snippets.

@TonyLadson

TonyLadson/ARF_sd.R

Created Apr 19, 2020
Embed
What would you like to do?
library(tidyverse)
library(optimx)
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 -------------------------------------------------------------------
## Plots against duration
AEP <- c(0.005)
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',
subtitle = "AEP = 0.5%") +
theme_grey(base_size = 5) +
theme(legend.position = c(0.6, 0.3)) +
theme(legend.key.size = unit(2, 'mm')) +
annotate(geom = 'text',x = 400, y = 0.02, label = 'Short duration', hjust = 'right', size = 1.5) +
annotate(geom = 'text',x = 2640, y = 0.02, label = 'Long duration', hjust = 'left', size = 1.5) +
annotate(geom = "segment", x = 720, y = 0.02, xend = 420, yend = 0.02,
arrow = arrow(angle = 20, length = unit(1, 'mm'), ends = "last", type = "closed" ),
size = 0.2) +
annotate(geom = "segment", x = 1440, y = 0.02, xend = 2440 , yend = 0.02,
arrow = arrow(angle = 20, length = unit(1, 'mm'), ends = "last", type = "closed" ),
size = 0.2)
# These plots look odd on-screen but are fine when saved and loaded into a blog
# The warning messages relate to attempting to calculate short duration ARFs for areas
# larger than 1000 km-squared.
p
ggsave(file.path('figure','ARF-duration-smallAEP.png'), plot = p, width = 4, height = 3)
# Derivative of the Short Duration ARF equation
# First Derivative with respect to duration
DARF_short <- function(A, D, P){
(0.10332 * (A^0.265 - 0.190655 * log(D)))/(D^1.36) +
(0.0002825 *A^0.226 * (log(P)/log(10) + 0.3))/D^0.875 -
9.46938e-7*10^(-0.0000145833 * (D - 180)^2) * A^0.213 * (D - 180) * (log(P)/log(10) + 0.3) +
0.0547181/D^1.36
}
# Second derivative with respect to duration
D2ARF_short <- function(A, D, P){
#Second derivative
-(0.140515 * (A^0.265 - 0.190655 * log(D)))/D^2.36 -
(0.000247188 * A^0.226 * (log(P)/log(10) + 0.3))/D^1.875 +
0.0141 * A^0.213 * (4.5103e-9 * 10^(-0.0000145833 * (D - 180)^2) * (D - 180)^2 - 0.0000671587 * 10^(-0.0000145833 * (D - 180)^2)) * (log(P)/log(10) + 0.3) -
0.0941151/D^2.36
}
# Create plot that shows the short duration ARF and its derivative
duration <- 10^seq(0,4,length.out = 500)
area <- c(100)
aep <- 0.005
p1 <- expand_grid(duration, area, aep) %>%
rowwise() %>%
mutate(ARF_s = ARF_short(area = area, duration = duration, aep = aep) ) %>%
ggplot(aes(x = duration, y = ARF_s)) +
geom_line(size = 0.3) +
geom_vline(xintercept = 180, linetype = 'dotted', size = 0.2) +
scale_x_continuous(name = 'Duration (min)', trans = 'log10', limits = c(NA, 720), breaks = c(1, 10, 100, 180, 720)) +
scale_y_continuous(name = 'Areal Reduction Factor') +
labs(title = 'Tasmania',
subtitle=expression(paste("0.5% AEP, Area = 100 ", km^2))) +
theme_grey(base_size = 5)
p1
xdf_ARFs <- expand_grid(duration, area, aep) %>%
rowwise() %>%
mutate(ARF_s = ARF_short(area = area, duration = duration, aep = aep))
p2 <- expand_grid(duration, area, aep) %>%
rowwise() %>%
mutate(DARF_s = DARF_short(A = area, D = duration, P = aep)) %>%
ggplot(aes(duration, DARF_s)) +
geom_line(colour = 'red', size = 0.2) +
geom_vline(xintercept = 180, linetype = 'dotted', size = 0.2) +
#geom_line(data = xdf_ARFs, aes(duration, ARF_s)) +
scale_x_continuous(name = 'Duration (min)', trans = 'log10', limits = c(NA, 720), breaks = c(1, 10, 100, 180, 720)) +
scale_y_continuous(name = 'Derivative(Areal Reduction Factor)', trans = 'log10') +
theme_grey(base_size = 5)
p2
gridExtra::grid.arrange(p1, p2)
p <- gridExtra::arrangeGrob(p1, p2)
ggsave(file.path('figure','ARF-derivative.png'), plot = p, width = 4, height = 3)
# Breaking the short duration equation up into three parts
# Constants for short duration ARF
a <- 0.287
b <- 0.265
c <- 0.439
d <- 0.36
e <- 0.00226
f <- 0.226
g <- 0.125
h <- 0.0141
i <- -0.021
j <- 0.213
ARF_short <- function(area, duration, aep) {
min(1, (1 - a*(area^b - c*log10(duration)) * duration^(-d) +
e*area^f*duration^g * (0.3 + log10(aep)) +
h * area^j * 10^(i* (1/1440) * (duration - 180)^2) * (0.3 + log10(aep))))
}
ARF_short1 <- function(area, duration, aep){
-a*(area^b - c*log10(duration)) * duration^(-d)
}
ARF_short2 <- function(area, duration, aep){
e*area^f*duration^g * (0.3 + log10(aep))
}
ARF_short3 <- function(area, duration, aep){
h * area^j * 10^(i* (1/1440) * (duration - 180)^2) * (0.3 + log10(aep))
}
# Test the equations
ARF_short1(area = 1000, duration = 200, aep = 0.01)
ARF_short2(area = 1000, duration = 200, aep = 0.01)
ARF_short3(area = 1000, duration = 200, aep = 0.01)
# Plot the three parts
duration <- 10^seq(0,4,length.out = 500)
p <- expand_grid(area = 1000, duration = duration, aep = 0.005) %>%
rowwise() %>%
mutate(s1 = ARF_short1(area, duration, aep)) %>%
mutate(s2 = ARF_short2(area, duration, aep)) %>%
mutate(s3 = ARF_short3(area, duration, aep)) %>%
mutate(s4 = 1+s1+s2+s3) %>%
pivot_longer(cols = c(s1, s2, s3, s4), names_to = 'component') %>%
ggplot(aes(x = duration, y = value, colour = component)) +
geom_line(size = 0.2) +
geom_vline(xintercept = 180, linetype = 'dotted', size = 0.2) +
scale_colour_hue(name = 'Eqn parts', labels = c('Part 1', 'Part 2', 'Part 3', 'Sum')) +
scale_x_continuous(name = 'Duration (min)',
trans = "log10",
breaks = c(1, 10, 100, 180, 720),
limits = c(NA, 720)) +
#theme(legend.position = c(0.65, 0.2)) +
theme_gray(base_size = 5) +
theme(legend.key.size = unit(2, 'mm'), legend.position = c(0.65, 0.2))
p
ggsave(file.path('figure','ARF-short_parts.png'), plot = p, width = 4, height = 3)
# test of the (D-180)^2 part
T1 <- 10^(-0.021* (1/1440)*(D-180)^2)
D <- 100:300
plot(T1)
## Plots against aep
# Make a sequence of AEP values evenly spaced in as quantiles of the normal distribution
z_seq <- seq(qnorm(0.5), qnorm(0.005), length.out = 200)
AEP <- pnorm(z_seq)
area <- c(5, 10, 100, 1000, 10000, 20000)
duration <- c(1, 5, 10, 100, 180, 720, 1000, 1440, 10080)
aep.labs <- 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)
aep_breaks = c(0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5)
duration.labs <- str_c(c(1, 5, 10, 100, 180, 720, 1000, 1440, 10080), ' minute')
names(duration.labs) <- c(1, 5, 10, 100, 180, 720, 1000, 1440, 10080)
xdf <- expand_grid(area = area, AEP = AEP, duration = duration)
p <- xdf %>%
rowwise() %>%
mutate(ARF_calc = ARF(area, duration, aep = AEP, region = 'Northern Coastal', params = params)) %>%
ggplot(aes(x = 1-AEP, y = ARF_calc, colour = factor(area))) +
geom_line(size = 0.2) +
scale_x_continuous(name = 'AEP(%)', trans = 'probit', breaks = 1-aep_breaks, labels = aep.labs) +
scale_y_continuous(name = 'Areal Reduction Factor') +
scale_colour_hue(name = expression(paste("Area ", km^2))) +
facet_wrap(~duration, labeller = labeller(duration = duration.labs)) +
theme_grey(base_size = 5) +
theme(legend.key.size = unit(1.5, 'mm'), legend.position = c(0.8,0.12))
p
ggsave(file.path('figure','ARF-AEP_facet.png'), plot = p, width = 4, height = 3)
# What combination of area and duration provides the greatest change in ARF as AEP goes from 50% to 0.05%
# Note that Short duration ARFs are not a function of region.
#https://stackoverflow.com/questions/44224761/getting-error-using-optim-function-with-bounds-in-r
# Uses penalties to enforce bounds to keep area between 0 and 1000
# and duration between 1 and 720
# The objective is the difference between the ARF at a low AEP and the ARF and a high AEP
# Find where this is a maximum
Obj <- function(par){
area <- par[1]
duration <- par[2]
if(area > 1000) return(100 + sum(par^2))
if(area < 0) return(100 + sum(par^2))
if(duration > 720) return(100 + sum(par^2))
if(duration < 1)return(100 + sum(par^2))
-ARF_short(area, duration, aep = 0.5) +
ARF_short(area, duration, aep = 0.005)
}
par.init <- (c(10, 10))
Obj(par.init)
optimx(par = par.init, fn = Obj, method = 'Nelder-Mead')
# I tried this using several different starting points and got similar results.
# > optimx(par = par.init, fn = Obj, method = 'Nelder-Mead')
# p1 p2 value fevals gevals niter convcode kkt1 kkt2 xtime
# Nelder-Mead 1000 183.4142 -0.1640775 161 NA NA 0 FALSE FALSE 0.001
# i.e. 1000 km-squared catchment and a duration of 183.4 min
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment