Calculation and graphical display of TQmean. Code to reproduce the figures in the blog https://tonyladson.wordpress.com/2018/08/20/tqmean-a-measure-of-the-impact-of-urbanisation-on-flow/

TQmean

library(tidyverse)
library(lubridate)
library(naniar)

# Reference

# TQmean is defined in Booth et al. (2004)
# Booth, D. B., J. R. Karr, S. Schauman, C. P. Konrad, S. A. Morley, M. G. Larson and S. J. Burges (2004). 
# "Reviving urban streams: land use, hydrology, biology and human behaviour." 
# Journal of the American Water Resources Association 40(5): 1351-1364.



# Constants
#_________________________________________________________________________

# Day of the year that is the first day of each month

month_starts <- yday(as.Date('2001-01-01') %m+% months(0:11))


# Functions
#_________________________________________________________________________

# Function to read and process a flow file
# This is based on file files from the BoM www.bom.gov.au

Read_flow_file <- function(fname) {
  flow <- read_csv(file = fname, skip = 9, 
                   col_types = cols(
                     `#Timestamp` = col_datetime(format = ""),
                     Value = col_double(),
                     `Interpolation Type` = col_integer(),
                     `Quality Code` = col_integer()
                   ))
  
  
  flow %>% 
    select(xdate = `#Timestamp`, # Rename
           flow_cumec = Value,
           interp_type = `Interpolation Type`,
           qcode = `Quality Code`) %>% 
    mutate(xdate = as.Date(xdate),  # Add columns we use later
           xyear = year(xdate),
           xmonth = month(xdate),
           jday = yday(xdate))
  
}

# Function to set the breaks on a graph at each year

Year_break = function (...){
  function(x) {
    minx         = floor(min(x)) - 1;
    maxx         = ceiling(max(x)) + 1;
    major_breaks = seq(minx, maxx, by=1)
    return(major_breaks)
  }
}

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


# Data files on Google Drive

# Brushy Creek
#https://drive.google.com/open?id=165YuTKX-CLl2ucZgoN_uPrm27sR1n2kg

# Olinda Ck
#https://drive.google.com/open?id=1S3ECfYZMna_KsI6gr7ABxc0Be4mWt5k5

# read in data
#____________________________________________________________________________________________________________
id <- "165YuTKX-CLl2ucZgoN_uPrm27sR1n2kg" # Brushy
Brushy <- Read_flow_file(fname = sprintf("https://docs.google.com/uc?id=%s&export=download", id   ))


id <- "1S3ECfYZMna_KsI6gr7ABxc0Be4mWt5k5"
Olinda <- Read_flow_file(fname = sprintf("https://docs.google.com/uc?id=%s&export=download", id   ))

# Check and adjust data
#____________________________________________________________________________________________________________

# Remove 2 ML/d from Olinda
# See https://www.melbournewater.com.au/waterdata/waterwaydiversionstatus/Documents/Olinda_Creek_catchment.pdf
# Silvan Reservoir provides 2 ML/d to Olinda Ck

Olinda <- Olinda %>% 
  mutate(flow_cumec = pmax(0, flow_cumec - 2/86.4)) 

# Plot missing

Brushy  %>% 
  ggplot(aes(x = xdate, y = flow_cumec)) +
  geom_miss_point()


Olinda  %>% 
  ggplot(aes(x = xdate, y = flow_cumec)) +
  geom_miss_point()


# Zero flows

Olinda %>% 
  summarise(mean(near(flow_cumec, 0), na.rm = TRUE))

Brushy %>% 
  summarise(mean(near(flow_cumec, 0), na.rm = TRUE))


# TQ mean calculation


# Including the years with missing data (just taking the mean of the flows that are available in that year)
Olinda %>% 
  group_by(xyear) %>% 
  summarise(flow_mean = mean(flow_cumec, na.rm = TRUE), TQmean = mean(flow_cumec > flow_mean, na.rm = TRUE)) %>%  
  summarise(mean(TQmean, na.rm = TRUE))
# 0.37



# Excluding years with missing data
Olinda %>% 
  group_by(xyear) %>% 
  summarise(flow_mean = mean(flow_cumec, na.rm = TRUE), TQmean = mean(flow_cumec > flow_mean)) %>% 
  summarise(mean(TQmean, na.rm = TRUE))
# 0.35



# Including the years with missing data (just taking the mean of the flows that are available in that year)
Brushy %>% 
  group_by(xyear) %>% 
  summarise(flow_mean = mean(flow_cumec, na.rm = TRUE), TQmean = mean(flow_cumec > flow_mean, na.rm = TRUE))   %>% 
  summarise(mean(TQmean, na.rm = TRUE))
# 0.209



# Excluding years with missing data
Brushy %>% 
  group_by(xyear) %>% 
  summarise(flow_mean = mean(flow_cumec, na.rm = TRUE), TQmean = mean(flow_cumec > flow_mean)) %>% 
  summarise(mean(TQmean, na.rm = TRUE))
# 0.207


# Tile plot of Tqmean


Olinda %>% 
  filter(xyear %in% c(1988:2016)) %>% 
  group_by(xyear) %>% 
  mutate(flow_mean = mean(flow_cumec, na.rm = TRUE),
         Above_mean = flow_cumec > flow_mean) %>% 
  ggplot(aes(x = jday, y = xyear, fill = Above_mean)) +
  geom_tile() +
  scale_fill_discrete(name = NULL, 
                      breaks = c(TRUE, FALSE, NA), 
                      labels = c("Above mean", "Below mean", "Missing"), 
                      l = 40) +
  scale_x_continuous(name = '', breaks = month_starts, labels= month.abb, expand = c(0,0)) +
  scale_y_continuous(name = '', breaks = Year_break(), expand = c(0,0) ) +
  theme_grey() +
  theme(legend.position = 'none') +
  annotate(geom = 'label', x = 8, y = 2015.5, label = 'A')


Brushy %>% 
  filter(xyear %in% c(1988:2016)) %>% 
  group_by(xyear) %>% 
  mutate(flow_mean = mean(flow_cumec, na.rm = TRUE),
         Above_mean = flow_cumec > flow_mean) %>% 
  ggplot(aes(x = jday, y = xyear, fill = Above_mean)) +
  geom_tile() +
  scale_fill_discrete(name = NULL, 
                      breaks = c(TRUE, FALSE, NA), 
                      labels = c("Above mean", "Below mean", "Missing"),
                      l = 40) +
  scale_x_continuous(name = '', breaks = month_starts, labels= month.abb, expand = c(0,0)) +
  scale_y_continuous(name = '', breaks = Year_break(), expand = c(0,0) ) +
  theme_grey() +
  theme(legend.position = 'bottom') +
  annotate(geom = 'label', x = 8, y = 2015.5, label = 'B')