Hydstra rating plots explained. Code to accompany the blog at https://tonyladson.wordpress.com/2016/11/29/hydstra-rating-plots-explained/

Hydstra_rating.R

library(tidyverse)
library(stringr)
library(RColorBrewer)
library(scales)
library(lubridate)
library(readxl)


# Avon River at Stratford, gauge 225201
# Stream gauging and rating data is available at http://www.bom.gov.au/waterdata/

# Load current rating table from the Bureau of Meteorology site
# I've saved some example data to dropbox

rating_bom <- read_csv('https://dl.dropboxusercontent.com/u/10963448/Avon_ratingtable_BOM.csv', skip = 4)


# Load current rating table from the Victorian data warehouse
# http://data.water.vic.gov.au/monitoring.htm?ppbm=225201&rs&1&rspf_org


rating_warehouse <- read_csv('https://dl.dropboxusercontent.com/u/10963448/Avon-ratingtable_warehouse.csv', skip = 0)


# Load Avon Gaugings from the BoM
# http://www.bom.gov.au/waterdata/

gaugings <- read_csv('https://dl.dropboxusercontent.com/u/10963448/Avon_gaugings.csv', skip = 2)


# Function to define the transformation 
# of the gauge heights 
# gh = gauge height
# CTF = gauge height at cease to flow 

y_trans <- function(gh, CTF){
  log10(gh - CTF)
}


# Prepare and tidy data
# Convert discharge in cumec to ML/d
# Remove any levels < 1.225
# Transform the levels

gaugings <- gaugings %>% 
  mutate(Date = as.Date(Date, format = '%d/%m/%Y'),
         discharge_MLd = discharge_cumec*86.4) %>% 
  filter(level_m > 1.225) %>% 
  mutate(y_Hydstra = y_trans(level_m, CTF = 1.225))

rating_bom <- rating_bom %>% 
  filter(level_m > 1.225) %>% 
  mutate(discharge_MLd = discharge_cumec*86.4,
         y_Hydstra = y_trans(level_m, CTF = 1.225))

rating_warehouse <- rating_warehouse %>% 
  filter(GH > 1.225) %>% 
  mutate(discharge_MLd = as.numeric(flow),
         y_Hydstra = y_trans(GH, CTF = 1.225))
  


# Labeller function to label the date of gaugines
# see http://stackoverflow.com/questions/21311489/scatter-plot-with-ggplot2-colored-by-dates/39283368#39283368

as.Date_origin <- function(x){
  format(as.Date(x, origin = '1970-01-01'), format = '%Y')
}

my.x.breaks <- 10^(-2:6)
my.x.labels <- formatC(my.x.breaks, format = 'fg', digits = 1, big.mark = ',')

my.y.breaks <- -3:1
my.y.labels <- 10^(my.y.breaks) + 1.225



Hydstra_plot <- gaugings %>% 
  ggplot(aes(x = discharge_MLd, y = y_Hydstra)) + 
  geom_point(aes(colour = as.integer(Date))) +
  scale_x_log10(name = 'Discharge ML/d', breaks = my.x.breaks, labels = my.x.labels ) +
  scale_y_continuous(name = 'Stream water level (m)', breaks = my.y.breaks, labels = my.y.labels, limits = c(-3,1)) +
  scale_colour_gradientn(name = 'Date', colours=c('black', 'blue', 'red'), labels=as.Date_origin) 
  
# Add BoM rating
  
Hydstra_plot = Hydstra_plot +
  geom_line(data = rating_bom, aes(x = discharge_MLd, y = y_Hydstra)) 

# Add warehouse rating and plot

Hydstra_plot + geom_line(data = rating_warehouse, aes(x = discharge_MLd, y = y_Hydstra), colour = 'blue', alpha = 0.5, size = 2 ) +
  annotate(geom = 'segment', x = 1000, xend = 10000, y = -2.5, yend = -2.5, colour = 'blue', alpha = 0.5, size = 2) +
  annotate(geom = 'text', x = 11000, y = -2.5, label = 'Warehouse rating', hjust = 0) +
  annotate(geom = 'segment', x = 1000, xend = 10000, y = -2.7, yend = -2.7, colour = 'black', size = 1) +
  annotate(geom = 'text', x = 11000, y = -2.7, label = 'Bureau rating', hjust = 0)