Antecedent Precipitation Index calculations see https://tonyladson.wordpress.com/2016/06/14/how-wet-is-the-catchment/

API.R

#
#  API  - Antecedent Precipitation Index
#
#
#

library(stringr)
library(readr)
library(dplyr)
library(lubridate)
library(timeDate)
library(ggplot2)
library(forecast)
library(hexbin)
library(repmis)
library(quantreg)
library(splines)
library(devtools)


devtools::install_github("slowkow/ggrepel")  # labelling of ggplot
library(ggrepel)

#remove(list = objects())

# Read in Silo rainfall data
# We just need a data frame of dates and daily rainfalls to run the script
# The data that I used is available at https://drive.google.com/file/d/101bPqP4YkxEwaEirG2-FSgNMPu2hWhPB/view?usp=sharing

my.path <- "/Users/anthonyladson/Dropbox/Data_misc" # Needs to be changed to be revelant for your computer
my.file <- "Creswick.txt"
fname <- str_c(my.path, my.file, sep = '/')
creswick <- read_table(fname, skip = 49, col_names = FALSE)
creswick.rain <- creswick[ , c(1,8)] # retain date and rain i.e. column 1 and column 8
names(creswick.rain) <- c('xDate', 'rain')


rain <- creswick.rain$rain # pull out the daily rainfall for API calculations

# Calculate API
# 
api <- numeric(length(rain)) # set up the api vector
api[1] <- 0.95*rain[1] # set the first api value to the rainfall on day 1


# loop through the vector of rainfalls to calculate API for each day
for(i in 2:length(rain)){
  api[i] <- 0.95*api[i-1] + rain[i]
}


#  Function to calculate date as an angle so we can make a polar plot


# When I ran this code there is an issue with ceiling_date and the way it handles
# dates 1 Jan before 1970
# 


Date_to_Angle <- function(my.date){
  
  start.year <- floor_date(my.date, 'year')
  end.year <- ceiling_date(my.date, 'year')
  length.year <- as.numeric(difftime(end.year, start.year, units = 'sec'))
  time.since.year.start <- as.numeric(difftime(my.date, start.year, units = 'sec'))
  
  2 * pi * time.since.year.start/length.year # date as angle
  
}




# Tidy all results into a data frame
creswick.api <- creswick.rain %>% 
  mutate(API = api) %>% 
  mutate(xDate = parse_date_time(xDate, orders = 'Ymd'))  %>% 
  mutate(my.year = year(xDate)) %>% 
  mutate(my.month = month(xDate)) %>% 
  mutate(my.yday = yday(xDate)) %>%
  mutate(my.angle = Date_to_Angle(xDate)) %>% 
  mutate(my.angle = ifelse(is.na(my.angle), 0, my.angle)) # work around for problem in ceiling_date for dates prior to 1970




# peak API values

creswick.api %>% 
  arrange(desc(API)) %>% 
  View

# Max API per year
creswick.api %>% 
  group_by(my.year) %>% 
  summarise(max.API = max(API)) %>% 
  arrange(desc(max.API)) %>% 
  View



# Preliminary plots



# All data with over-plotting controlled by alpha
# Add quantiles
creswick.api %>% 
  ggplot(aes(x = my.yday, y = API)) + 
  geom_point(size = 2, colour = alpha("black", 1/50)) 

#__________________________________________________________________________________________________________

# Figure 1
# see https://gist.github.com/TonyLadson/2d63ca70eef92583001dece607127759
# https://tonyladson.wordpress.com/2016/06/20/fitting-non-linear-models/



#__________________________________________________________________________________________________________

# Figure 2


# 2010 API from Aug to Nov

# creswick.api %>% 
#   filter(my.year == 2010) %>% 
#   filter(my.month %in% c(8:10)) %>% 
#   View


month.starts <- yday(seq(ymd('2010-01-01'),ymd('2010-12-31'), by='months'))
month.starts


creswick.api %>% 
  filter(my.year == 2010) %>% 
  filter(my.month %in% c(8:10)) %>% 
  ggplot(aes(x = my.yday, y = API)) + 
  geom_point() +
  geom_line() +
  geom_step(aes(y = rain), colour = 'Blue') +
  ylab('Rainfall, API (mm)') +
  scale_x_continuous(name = '', breaks = month.starts, labels = month.abb) +
  annotate("text", x = yday('2010-08-01'), y = 30, label = 'Rainfall', colour = 'Blue',  fontface = 'bold', hjust = 'left') +
  annotate("text", x = yday('2010-08-01'), y = 90, label = 'API', colour = 'Black', fontface = 'bold', hjust = 'left') +
  theme_bw() + 
  theme(
    panel.background = element_rect(fill="grey98"),
    axis.title.x = element_text(colour="grey20", size=14, margin=margin(20,0,0,0)),
    axis.text.x = element_text(colour="grey20",size=12),
    axis.title.y = element_text(colour="grey20",size=14, margin = margin(0,20,0,0)),
    axis.text.y = element_text(colour="grey20",size=14),
    legend.title = element_text(colour="grey20",size=12),
    plot.margin = unit(c(0.5, 0.5, 1, 0.5), "cm")) # top, right, bottom, left


#__________________________________________________________________________________________________________

# Figure 3


# Plot for 2010 including background API values


api.2010 <- creswick.api %>% 
  filter(my.year %in% c(2010)) 


creswick.api%>% 
  ggplot(aes(x = my.yday, y = API)) + 
  geom_point(colour = alpha('grey80', 1/5)) +
  geom_line(data = api.2010, aes(x = my.yday,  y = API)) +
  geom_step(data = api.2010, aes(y = rain), colour = 'Blue') +
  scale_x_continuous(name = '', breaks = month.starts, labels = month.abb) +
  ylab('Rainfall, API (mm)') +
  annotate("text", x = yday('2010-01-01'), y = 20, label = 'Rainfall', colour = 'Blue',  fontface = 'bold', hjust = 'left', size = 3) +
  annotate("text", x = yday('2010-01-01'), y = 55, label = 'API', colour = 'Black', fontface = 'bold', hjust = 'left', size = 3) +
  theme_bw() + 
  theme(
    panel.background = element_rect(fill="white"),
    axis.title.x = element_text(colour="grey20", size=14, margin=margin(20,0,0,0)),
    axis.text.x = element_text(colour="grey20",size=12),
    axis.title.y = element_text(colour="grey20",size=14, margin = margin(0,20,0,0)),
    axis.text.y = element_text(colour="grey20",size=14),
    legend.title = element_text(colour="grey20",size=12),
    plot.margin = unit(c(0.5, 0.5, 1, 0.5), "cm")) # top, right, bottom, left




#__________________________________________________________________________________________________________

# Figure 4


# Polar plot

# Angles for month start

month.start <- yday(seq(as.Date('2010-01-01'),as.Date('2010-12-31'), by='1 months'))
month.start.radians <- 2*pi*month.start/365.25 


lab.df <- creswick.api %>% 
  group_by(my.year) %>%   # grab only 1 from each year
  summarise(max.API = max(API), max.date = xDate[which.max(API)], max.angle = my.angle[which.max(API)]) %>% 
  arrange(desc(max.API)) %>% # sort in order of API from largest to smallest
  slice(1:7) 




creswick.api %>% 
  ggplot(aes(x=my.angle, y=API)) + 
  coord_polar('x', start = 0, direction = 1) + 
  geom_point(colour = alpha("grey10", 1/2)) +
  scale_x_continuous(name = '', breaks = month.start.radians, labels = month.name, limits = c(0, 2*pi)) +
  scale_y_continuous(name = 'API (mm)') +
  geom_text_repel(data = lab.df, aes(x = max.angle, y = max.API + 20), label = lab.df$my.year, segment.size = 0, colour = 'Blue') +
  theme_bw() + 
  theme(
    panel.background = element_rect(fill="white"),
    #axis.title.x = element_text(colour="grey20", size=8, margin=margin(20,0,0,0)),
    axis.text.x = element_text(colour="grey30",size=10),
    axis.title.y = element_text(colour="grey20",size=14, margin = margin(0,20,0,0)),
    axis.text.y = element_text(colour="grey20",size=14),
    #legend.title = element_text(colour="grey20",size=12),
    plot.margin = unit(c(0.5, 0.5, 1, 0.5), "cm")) # top, right, bottom, left



#__________________________________________________________________________________________________________

# Figure 5



# 2010 and 2011
# Plot from 1 July 2010 to 1 Jun 2011


yday('2010-07-01') # day number of 1 July
#182

month.starts <- yday(seq(ymd('2010-01-01'),ymd('2010-12-31'), by='months'))
month.starts.shift <- month.starts - 182
month.starts.shift <- ifelse(month.starts.shift < 0, month.starts.shift + 365, month.starts.shift)
month.starts.shift <- sort(month.starts.shift)

month.labels.shift <- month.abb[c(7:12, 1:6)]

creswick.api.2010_2011 <- creswick.api %>% 
  filter( (my.year == 2011 & my.month %in% c(1:6)) | (my.year == 2010 & my.month %in% c(7:12)) ) %>% 
  mutate(my.yday.shift =  my.yday - 182) %>% 
  mutate(my.yday.shift = ifelse(my.yday.shift < 0, my.yday.shift + 365, my.yday.shift)) 




creswick.api %>% 
  mutate(my.yday.shift =  my.yday - 182) %>% 
  mutate(my.yday.shift = ifelse(my.yday.shift < 0, my.yday.shift + 365, my.yday.shift)) %>% 
  ggplot(aes(x = my.yday.shift, y = API)) + 
  geom_point(colour = alpha('grey80', 1/5)) +
  scale_x_continuous(name = '', breaks = month.starts.shift, labels = month.labels.shift) +
  geom_line(data = creswick.api.2010_2011, aes(x = my.yday.shift, y = API)) +
  geom_point(data = creswick.api.2010_2011, aes(x = my.yday.shift, y = API)) +
  geom_step(data = creswick.api.2010_2011, aes(x = my.yday.shift, y = rain), colour = 'Blue') +
  ylab('Rainfall, API (mm)') +
  annotate("text", x = 0, y = 20, label = 'Rainfall', colour = 'Blue',  fontface = 'bold', hjust = 'left', size = 4) +
  annotate("text", x = 0, y = 70, label = 'API', colour = 'Black', fontface = 'bold', hjust = 'left', size = 4) +
  theme_bw() + 
  theme(
    panel.background = element_rect(fill="white"),
    axis.title.x = element_text(colour="grey20", size=14, margin=margin(20,0,0,0)),
    axis.text.x = element_text(colour="grey20",size=12),
    axis.title.y = element_text(colour="grey20",size=14, margin = margin(0,20,0,0)),
    axis.text.y = element_text(colour="grey20",size=14),
    legend.title = element_text(colour="grey20",size=12),
    plot.margin = unit(c(0.5, 0.5, 1, 0.5), "cm")) # top, right, bottom, left



# Rain in Jan 2011

creswick.api %>% 
  filter(my.year == 2011) %>%
  filter(my.yday %in% c(11:14)) %>% 
  summarise(sum(rain))


# 4 day rainfal totals

creswick.api %>% 
  mutate(fday  = rain + lag(rain, 1) + lag(rain, 2) +  lag(rain, 3)  ) %>% 
  arrange(desc(fday))



#__________________________________________________________________________________________________________

# Figure 6

# Spells plot


creswick.api %>%
  filter(my.year > 1995) %>% 
  mutate(high = ifelse(API > 90, 1, 0)) %>% # We are interested in flows above 100
  ggplot(aes(x = my.yday, y = my.year)) + 
  geom_tile(aes(fill = factor(high, levels = c(1,0) ))) +
  scale_fill_manual(name = "High API", values=c('Red', 'light blue'), labels = c(' > 90 mm', ' < 90 mm')) +
  scale_x_continuous(name = '', breaks = month.starts, labels = month.abb) +
  scale_y_continuous(name = '', breaks = 1990:2016, labels = 1990:2016) +
  theme_bw() + 
  theme(
    panel.background = element_rect(fill="white"),
    axis.title.x = element_text(colour="grey20", size=14, margin=margin(20,0,0,0)),
    axis.text.x = element_text(colour="grey20",size=12),
    axis.title.y = element_text(colour="grey20",size=14, margin = margin(0,20,0,0)),
    axis.text.y = element_text(colour="grey20",size=14),
    legend.title = element_text(colour="grey20",size=12),
    plot.margin = unit(c(0.5, 0.5, 1, 0.5), "cm")) # top, right, bottom, left