Record Australian Rainfalls (and an envelope curve) https://tonyladson.wordpress.com/2016/02/08/envelop-curve-for-record-australian-rainfall/

AusRecordRain.R

# See https://tonyladson.wordpress.com/2016/02/08/envelop-curve-for-record-australian-rainfall/
# remove(list = objects())

library(XML)
library(dplyr)
library(MASS)
library(xtable)



# Australian Rainfall Records
my.url <- 'http://www.bom.gov.au/water/designRainfalls/rainfallEvents/ausRecordRainfall.shtml'


# Read in the tables for each jurisdiction
au.rain <- readHTMLTable(my.url, 
                         header = FALSE, 
                         colClasses = c('character', 'character', 'numeric', 'numeric', 'character'),
                         encoding = "UTF-8",  )

# set header = FALSE so column names are ignored, there are small inconsistencies
# in column names which was causing trouble
# set encoding to "UTF-8" which seem to be what the Bureau is using

# name the juridication for each data frame
names(au.rain) <- c('ACT', 'NSW', 'NT', 'QLD', 'SA', 'Tas', 'Vic', 'WA')


# bind into a single data frame
au.rain <- do.call('rbind', au.rain) 

# add row names as a column (so we can find where values come from)

au.rain <- au.rain %>% mutate(source.row = row.names(au.rain))

# name columns

names(au.rain)[1:5] = c('location', 'date', 'duration.min', 'rainfall.mm', 'source')


# get maximum rainfall for each duration

au.rain.max <- au.rain %>%
  group_by(duration.min) %>%
  slice(which.max(rainfall.mm)) 

# Calculate duration in hours

au.rain.max <- au.rain.max %>%
  mutate(duration.hour = as.vector(duration.min/60))

# create an increasing series
# start with the first value and only accept the next value if its larger

au.rain.max <- au.rain.max %>%
  ungroup() %>%
  mutate(running.max = cummax(rainfall.mm)) %>%
  filter(rainfall.mm >= running.max)



# create a HTML table in include in blog

au.rain.max
rain.table <- au.rain.max[, 1:4]

rep('c', 5)
rain.table <- xtable(au.rain.max[, 1:4 ], digits = 0, align = rep('c',5))
print(rain.table, type = 'html')

# Calculate Envelope Curve
# We'll use robust regression

m.1 <- MASS::rlm(log(rainfall.mm) ~ log(duration.hour), data = au.rain.max)
coef(m.1)
# (Intercept) log(duration.hour) 
# 5.2959970          0.5696566 

# So the slope we need is 0.5696566
# The line needs to pass through the 5th data point which determines the intercept
au.rain.max[5, ]

# location        date duration.min rainfall.mm source source.row duration.hr duration.hour running.max
# 1    Nobby 11 Mar 1939           20         203            QLD.10   0.3333333     0.3333333         203

env.intercept <- as.vector(203/((1/3)^coef(m.1)[2]))
env.intercept
# [1] 379.5695

# Envelope curve
#380*d^0.57






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

#remove(list = objects())

# World Rainfall Records and Envelope curve
# download CSV file extracted from Annex II of WMO (2009)

# WMO (2009) Manual on estimation of probable maximum precipitation (PMP). WMO-No. 1045.  
# World Meteorological Organization. http://www.wmo.int/pages/prog/hwrp/publications/PMP/WMO%201045%20en.pdf


library(repmis)
library(dplyr)

my.url <- 'https://dl.dropboxusercontent.com/u/10963448/RecordRainfall-World.csv'

world.rain <- tbl_df(repmis::source_data(my.url, skip = 3, sep = ','))

world.rain <- world.rain %>%
  mutate(duration.hour = `Duration (min)`/60) # work in hours



# World Envelope Curve
# WMO (2009) 

# R = 491 D^0.452

env.world <- function(d) {
  491 * d^0.452
}

# Range of durations

# from 0.167 hours to 120 hours
d.seq.world <- seq(from = min(world.rain$duration.hour), to = max(world.rain$duration.hour), length.out = 100) 



# New Zealand rainfall records
# from Griffiths, G. A., A. I. McKerchar and C. P. Pearson (2014). 
#"Towards prediction of extreme rainfalls in New Zealand." 
# Journal of Hydrology (NZ) 53(1): 41-52.

nz.rain <- data_frame(
  rainfall = c(34, 134, 428, 566, 758, 1049, 2927, 18405),
  duration.hour = c(0.167, 1, 8.5, 12, 24, 48, 730, 8760)
)


# Envelope curve for NZ

# Tomlinson envelope curve
# The Griffith curve (below) seems to have superceeded the Tomlinson curve, so that will not be used here
# Tomlinson, A. I. (1980)
# The Frequency of high intensity rainfalls in New Zealand, Part 1.
# Water and Soil Technical Pulbication No. 19. Water and Soil Division
# Ministry of Workd and Development, Wellington New Zealand, 36 p

# env.tomlinson <- function(d) {
#   112 * d^0.55
# }
# 
# # from 0.167 hours to 120 hours
# d.seq.tomlinson <- seq(from = 0.167, to = 120, length.out = 100) 


# Griffiths envelope curve
# from Griffiths, G. A., A. I. McKerchar and C. P. Pearson (2014)

env.griffiths <- function(d) {
  145 * d^0.55
}

# from 0.167 hours to 8760 hours (1 years)
d.seq.griffiths <- seq(from = 0.167, to = 8760, length.out = 100) 

# Australian data and envelope curve 

# loading these in from dropbox

my.url <- 'https://dl.dropboxusercontent.com/u/10963448/RecordRainfall-au.csv'
au.rain.max <- tbl_df(repmis::source_data(my.url, skip = 3, sep = ','))


au.rain.max <- au.rain.max %>%
  mutate(duration.hour = `Duration (min)`/60) # work in hours


# Australian Envelope Curve

# R = 380 D^0.57

env.au <- function(d) {
  380 * d^0.57
}

# from 0.033 hours to 192 hours 
d.seq.au <- seq(from = 0.033, to = 192, length.out = 100) 



# Plot

par(oma = c(0,3,0,0))

plot(`Depth (mm)` ~ duration.hour, data = world.rain, 
     log = 'xy',
     xlab = 'Duration (hr)',
     ylab = '',
     xaxt = 'n',
     yaxt = 'n',
     pch = 21,
     bg = 'red')

my.labels = format(axTicks(1), big.mark = ',', scientific = FALSE, digits = 0, trim = TRUE)

axis(side = 1, at = axTicks(1), labels = my.labels)
axis(side = 2, at = axTicks(2), labels = format(axTicks(2), big.mark = ','), las = 2)

mtext(side = 2, text = 'Rainfall (mm)', line = 4 )


points(rainfall ~ duration.hour, data = nz.rain, pch = 21, bg = 'grey', col = 'black') # add NZ records

points(`Depth (mm)` ~ duration.hour, data = au.rain, pch = 21, bg = 'green', col = 'black') # add Australian records


#lines(d.seq.tomlinson, env.tomlinson(d.seq.tomlinson), lty = 2, col = 'blue', lwd = 2) # Tomlinson envelope curve (NZ)
lines(d.seq.griffiths, env.griffiths(d.seq.griffiths), lty = 3, col = 'grey', lwd = 2) # Griffiths envelope curve (NZ)
lines(d.seq.world, env.world(d.seq.world), lty = 2, col = 'red', lwd = 2) # World envelope curve
lines(d.seq.au, env.au(d.seq.au), lty = 2, col = 'green', lwd = 2) # Australian envelope curve

legend('topleft', 
       lty = 2, 
       pch = 21,
       col = c('red', 'green', 'grey'), 
       legend = c('World', 'Australia', 'New Zealand'), 
       pt.bg = c('red', 'green','grey'),
       bty = 'n'
)

text(0.5, 1000, col = 'red', lwd = 2, expression(paste('R = 491 ',D^0.452))) # World
text(50, 7000, col = 'dark green',  lwd = 2, expression(paste('R = 380 ',D^0.57))) # Australia
text(200, 500, col = grey(0.2),  lwd = 2, expression(paste('R = 145 ',D^0.55))) # New Zealand





# Review Australian envelope curve
# Use the slope of the world envelope curve
#  The critical point that the envelope curve passes through is from row 19 of the Australian data
#
au.rain.max[19, ]
# location                   date duration.min rainfall.mm source source.row duration.hour running.max
# 1 Bellenden Ker Top 5 Jan 1979         4320        2517    BoM     QLD.41            72        2517

env.intercept <- as.vector(2517/(72^0.452))
env.intercept
# [1] 364.2243

# Plot data, curve and annotate graph



par(oma = c(6,3,0,0))
plot(`Depth (mm)` ~ duration.hour, data = au.rain.max, 
     xlim = c(0.01, 200),
     ylim = c(1, 5000),
     log = 'xy',
     xlab = '',
     ylab = '',
     xaxt = 'n',
     yaxt = 'n',
     pch = 21, 
     bg = 'grey')



min.values <- c(1, 5, 10, 30, 60, 120, 6*60, 12*60, 24*60, 48*60, 96*60, 4*24*60, 7*24*60 )
min.ticks <- min.values/60
min.labels <- format(min.values, big.mark = ',', trim = TRUE)
axis(side = 1, at = min.ticks, labels = min.labels)  # label with minutes
mtext(side = 1, text = 'min', line = 1, adj = c(0,1))

hr.labels <- format(min.values/60, big.mark = ',', trim = TRUE, digits = 0)
hr.labels[as.numeric(hr.labels) < 1] <- NA
mtext(side = 1, text = 'hour', line = 3, adj = c(0, 1))

axis(side = 1, at = min.ticks, labels = hr.labels, line = 3) # label with hours

day.labels <- format(min.values/(60*24), big.mark = ',', trim = TRUE, digits = 0)
day.labels[as.numeric(day.labels) < 1] <- NA

axis(side = 1, at = min.ticks, labels = day.labels, line = 6) # label with hours
mtext(side = 1, text = 'day', line = 6, adj = c(0,1))

mtext(side = 1, text = 'Duration', line = 8 )

axis(side = 2, at = axTicks(2), labels = format(axTicks(2), big.mark = ','), las = 2)

mtext(side = 2, text = 'Rainfall (mm)', line = 4 )



# Add envelope curve

# env.au <- function(d) {
#   env.intercept * d^coef(m.1)[2]
# }

# Try with world slope

env.au <- function(d) {
  364.22 * d^0.452
 }


d.seq.au <- seq(from = 0.033, to = 200, length.out = 100)

lines(d.seq.au, env.au(d.seq.au), lty = 2, col = 'blue', lwd = 2)
text(0.06, 300, expression(paste('R = 364.2 ',D^0.452, ' (D in hours)')))