library(httr) #
library(jsonlite) #
library(tidyverse) #
library(lubridate) #
library(imputeTS) #
library(forecast) #
library(zoo) #
library(xts)
library(ggmap)
library(gmapsdistance)
library(gridExtra)
#Stations Locations
# SOUTH
# Princes 37001 = 56914
# Batemans 08052 = 56413
# Cooma 08082 = 56424
# NORTH
# Pacific 76001 = 57456
# Port 6128-PR = 99990008
# Byron 6116 = 15190087
# WEST
# BlueMt 99914 = 58847
# Bathurst 99001 = 58691
# Yass 94002 = 58049
# Find my stations
stations_api<- GET("https://api.transport.nsw.gov.au/v1/roads/spatial?format=geojson&q=select%20*%20from%20road_traffic_counts_station_reference%20",
verbose(),
encode="json",
add_headers(`Authorization` = "apikey fUa8N1LC42AYtVDKIt6jbAzQXFPcf9b31GYv"))
# Extract a clean stations dataframe from the raw API output
stations_raw <- rawToChar(stations_api$content)
stations_clean <- fromJSON(stations_raw)
stations_df <- as.data.frame(stations_clean[[2]])
stations <- as.data.frame(stations_df$properties)
# Recover
saveRDS(stations,file="stations.rds")
stations <- read_rds("/Users/RohanDanisCox/STDS/stations.rds")
# Write a 'for' loop to get my data and clean it
Station_code <- c(56914,56413,56424,57456,99990008,15190087,58847,58691,58049)
Station_data <-data.frame()
for(i in Station_code) {
url<- paste("https://api.transport.nsw.gov.au/v1/roads/spatial?format=geojson&q=select%20*%20from%20road_traffic_counts_hourly_permanent%20where%20station_key%20%3D%20",i,"%20",sep="")
get_api <- GET(url,
verbose(),
encode="json",
add_headers(`Authorization` = "apikey fUa8N1LC42AYtVDKIt6jbAzQXFPcf9b31GYv"))
take_raw <- rawToChar(get_api$content)
take_clean <- fromJSON(take_raw)
take_df <- as.data.frame(take_clean[[2]])
Station_data <- rbind(Station_data,take_df$properties)
}
Station_data$date <- gsub('.{10}$', '', Station_data$date)
Station_data$date <- as.Date(Station_data$date, format="%Y-%m-%d")
# Recover
saveRDS(Station_data,file="Station_data.rds")
Station_data <- read_rds("/Users/RohanDanisCox/STDS/Station_data.rds")
##############################
### PRINCES 56914 - SOUTH ####
##############################
# Filter to required station
Princes56914 <- Station_data %>%
filter(station_key == 56914) %>%
filter(cardinal_direction_seq == 5)
# Clean
Princes56914_dates <- seq.Date(as.Date("2015-10-18"),as.Date("2016-01-23"),by=1)
Princes56914 <- filter(Princes56914, !date %in% Princes56914_dates)
Princes56914 <- filter(Princes56914, daily_total>20000)
# Select the relevant variables
Princes56914 <- Princes56914 %>%
select(date,daily_total)
Princes56914 <- Princes56914[order(as.Date(Princes56914$date, format="%Y-%m-%d")),]
# Add in NAs for missing dates
Princes56914_dates <- seq.Date(as.Date(Princes56914[1,1]),as.Date(Princes56914[nrow(Princes56914),1]),by=1)
Princes56914_dates <-as.data.frame(Princes56914_dates)
names(Princes56914_dates)[names(Princes56914_dates)=="Princes56914_dates"] <- "date"
Princes56914 <- left_join(Princes56914_dates,Princes56914,by="date")
# Calculate NAs
Princes56914_rows <- nrow(Princes56914)
Princes56914_NAs <- sum(is.na(Princes56914$daily_total))/Princes56914_rows
# Show NAs
Princes56914_NA <-zoo(Princes56914$daily_total,Princes56914$date)
plotNA.distribution(Princes56914_NA, main="Distribution of NAs on Princes Highway near Heathcote")
# Address NAs
Princes56914$year <- format(as.Date(Princes56914$date, format="%d/%m/%Y"),"%Y")
Princes56914$day <- weekdays(as.Date(Princes56914$date))
Princes56914 <- Princes56914 %>%
group_by(year,day) %>%
mutate(average = round(mean(daily_total,na.rm = TRUE),0))
Princes56914$daily_total[is.na(Princes56914$daily_total)] <- Princes56914$average[is.na(Princes56914$daily_total)]
# Aggregate to Weekly
Princes56914_zoo <- zoo(Princes56914$daily_total,Princes56914$date)
Princes56914_week <- apply.weekly(Princes56914_zoo, mean)
Princes56914_ts <- msts(Princes56914_week,seasonal.periods=52.1785,start=2006)
Princes56914_tbats_full <- tbats(Princes56914_ts)
plot(forecast(Princes56914_tbats_full))
plot(forecast(Princes56914_tbats_full,h=104),include=104)
# Visualise Weekly data
Princes56914_plot <- window(Princes56914_ts,start=2006,end=2016)
ggseasonplot(Princes56914_plot)
# Decompose
plot(Princes56914_tbats_full)
Princes56914_season_tbats <- tbats.components(Princes56914_tbats_full)[,'season']
autoplot(Princes56914_season_tbats)
#or
Princes56914_decom <- decompose(Princes56914_ts, type="multiplicative")
autoplot(Princes56914_decom)
Princes56914_season_decom <- Princes56914_decom$seasonal
autoplot(Princes56914_season_decom)
# Model and test
Princes56914_train <- window(Princes56914_ts, end=2016+180/365.25)
Princes56914_test <- window(Princes56914_ts, start=2016+180/365.25)
Princes56914_tbats <- tbats(Princes56914_train)
autoplot(forecast(Princes56914_tbats,h=45),ylab="Heathcote Count") +
autolayer(Princes56914_test, series="Test Data",alpha=.8) +
theme_bw(base_size=20)
autoplot(forecast(Princes56914_tbats,h=52),include=104) +
autolayer(Princes56914_test, series="Test Data",alpha=.7)
accuracy(forecast(Princes56914_tbats),Princes56914_test)
# Test for ACF and PACF
Princes56914_tbats_random <- (tbats.components(Princes56914_tbats)[,'observed'])-(tbats.components(Princes56914_tbats)[,'level'])-(tbats.components(Princes56914_tbats)[,'season'])
plot(Princes56914_tbats_random)
par(mfrow = c(1, 2))
Acf(Princes56914_tbats_random,lag.max = 20,main="Heathcote ACF"
Pacf(Princes56914_tbats_random,lag.max = 20,main="Heathcote PACF")
###############################
### BATEMANS 56413 - SOUTH ####
###############################
# Filter to required station
Batemans56413 <- Station_data %>%
filter(station_key == 56413) %>%
filter(cardinal_direction_seq == 5)
# Select the relevant variables
Batemans56413 <- Batemans56413 %>%
select(date,daily_total)
Batemans56413 <- Batemans56413[order(as.Date(Batemans56413$date, format="%Y-%m-%d")),]
# Add in NAs for missing dates
Batemans56413_dates <- seq.Date(as.Date(Batemans56413[1,1]),as.Date(Batemans56413[nrow(Batemans56413),1]),by=1)
Batemans56413_dates <-as.data.frame(Batemans56413_dates)
names(Batemans56413_dates)[names(Batemans56413_dates)=="Batemans56413_dates"] <- "date"
Batemans56413 <- left_join(Batemans56413_dates,Batemans56413,by="date")
# Calculate NAs
Batemans56413_rows <- nrow(Batemans56413)
Batemans56413_NAs <- sum(is.na(Batemans56413$daily_total))/Batemans56413_rows
# Show NAs
Batemans56413_NA <-zoo(Batemans56413$daily_total,Batemans56413$date)
plotNA.distribution(Batemans56413_NA,main="Distribution of NAs on Princes Highway near Batemans Bay",cex.axis=1.7,cex.main=1.7)
# Address NAs
Batemans56413$year <- format(as.Date(Batemans56413$date, format="%d/%m/%Y"),"%Y")
Batemans56413$day <- weekdays(as.Date(Batemans56413$date))
Batemans56413 <- Batemans56413 %>%
group_by(year,day) %>%
mutate(average = round(mean(daily_total,na.rm = TRUE),0))
Batemans56413$daily_total[is.na(Batemans56413$daily_total)] <- Batemans56413$average[is.na(Batemans56413$daily_total)]
# Aggregate to Weekly
Batemans56413_zoo <- zoo(Batemans56413$daily_total,Batemans56413$date)
Batemans56413_week <- apply.weekly(Batemans56413_zoo, mean)
Batemans56413_ts <- msts(Batemans56413_week,seasonal.periods=52.1785,start=2008+169/365.25)
Batemans56413_tbats_full <- tbats(Batemans56413_ts)
plot(forecast(Batemans56413_tbats_full))
plot(forecast(Batemans56413_tbats_full,h=104),include=104)
# Visualise Weekly data
Batemans56413_plot <- window(Batemans56413_ts,start=2008+169/365.25,end=2016)
ggseasonplot(Batemans56413_plot)
# Decompose
plot(Batemans56413_tbats_full)
Batemans56413_season_tbats <- tbats.components(Batemans56413_tbats_full)[,'season']
autoplot(Batemans56413_season_tbats)
# or
Batemans56413_decom <- decompose(Batemans56413_ts, type="multiplicative")
autoplot(Batemans56413_decom)
Batemans56413_season_decom <- Batemans56413_decom$seasonal
autoplot(Batemans56413_season_decom)
# Model and test
Batemans56413_train <- window(Batemans56413_ts, end=2016+180/365.25)
Batemans56413_test <- window(Batemans56413_ts, start=2016+180/365.25)
Batemans56413_tbats <- tbats(Batemans56413_train)
autoplot(forecast(Batemans56413_tbats,h=45),ylab="Batemans Bay Count") +
autolayer(Batemans56413_test, series="Test Data",alpha=.8) +
theme_bw(base_size=20)
autoplot(forecast(Batemans56413_tbats,h=52),include=104) +
autolayer(Batemans56413_test, series="Test Data",alpha=.7)
accuracy(forecast(Batemans56413_tbats),Batemans56413_test)
# Test for ACF and PACF
Batemans56413_tbats_random <- (tbats.components(Batemans56413_tbats)[,'observed'])-(tbats.components(Batemans56413_tbats)[,'level'])-(tbats.components(Batemans56413_tbats)[,'season'])
plot(Batemans56413_tbats_random)
Acf(Batemans56413_tbats_random,lag.max = 20,main="Batemans Bay ACF")
Pacf(Batemans56413_tbats_random,lag.max = 20,main="Batemans Bay PACF")
############################
### COOMA 56424 - SOUTH ####
############################
# Filter to required stations
Cooma56424 <- Station_data %>%
filter(station_key == 56424 ) %>%
filter(cardinal_direction_seq == 5)
# Select the relevant variables
Cooma56424 <- Cooma56424 %>%
select(date,daily_total)
Cooma56424 <- Cooma56424[order(as.Date(Cooma56424$date, format="%Y-%m-%d")),]
# Add in NAs for missing dates
Cooma56424_dates <- seq.Date(as.Date(Cooma56424[1,1]),as.Date(Cooma56424[nrow(Cooma56424),1]),by=1)
Cooma56424_dates <-as.data.frame(Cooma56424_dates)
names(Cooma56424_dates)[names(Cooma56424_dates)=="Cooma56424_dates"] <- "date"
Cooma56424 <- left_join(Cooma56424_dates,Cooma56424,by="date")
# Calculate NAs
Cooma56424_rows <- nrow(Cooma56424)
Cooma56424_NAs <- sum(is.na(Cooma56424$daily_total))/Cooma56424_rows
# Show NAs
Cooma56424_NA <-zoo(Cooma56424$daily_total,Cooma56424$date)
plotNA.distribution(Cooma56424_NA,main="Distribution of NAs on Monaro Highway near Cooma")
# Address NAs
Cooma56424$year <- format(as.Date(Cooma56424$date, format="%d/%m/%Y"),"%Y")
Cooma56424$day <- weekdays(as.Date(Cooma56424$date))
Cooma56424 <- Cooma56424 %>%
group_by(year,day) %>%
mutate(average = round(mean(daily_total,na.rm = TRUE),0))
Cooma56424$daily_total[is.na(Cooma56424$daily_total)] <- Cooma56424$average[is.na(Cooma56424$daily_total)]
# Aggregate to Weekly
Cooma56424_zoo <- zoo(Cooma56424$daily_total,Cooma56424$date)
Cooma56424_week <- apply.weekly(Cooma56424_zoo, mean)
Cooma56424_ts <- msts(Cooma56424_week,seasonal.periods=52.1785,start=2008+259/365.25)
Cooma56424_tbats_full <- tbats(Cooma56424_ts)
plot(forecast(Cooma56424_tbats_full))
plot(forecast(Cooma56424_tbats_full,h=104),include=104)
# Visualise Weekly data
Cooma56424_plot <- window(Cooma56424_ts,start=2009,end=2017)
ggseasonplot(Cooma56424_plot)
# Decompose
plot(Cooma56424_tbats_full)
Cooma56424_season_tbats <- tbats.components(Cooma56424_tbats_full)[,'season']
autoplot(Cooma56424_season_tbats)
# or
Cooma56424_decom <- decompose(Cooma56424_ts, type="multiplicative")
autoplot(Cooma56424_decom)
Cooma56424_season_decom <- Cooma56424_decom$seasonal
autoplot(Cooma56424_season_decom)
# Model and test
Cooma56424_train <- window(Cooma56424_ts, end=2016+180/365.25)
Cooma56424_test <- window(Cooma56424_ts, start=2016+180/365.25)
Cooma56424_tbats <- tbats(Cooma56424_train)
autoplot(forecast(Cooma56424_tbats,h=45),ylab="Cooma Count") +
autolayer(Cooma56424_test, series="Test Data",alpha=.8) +
theme_bw(base_size=20)
Cooma_yearly <- autoplot(forecast(Cooma56424_tbats,h=45),ylab="Cooma Count",include=13) +
autolayer(Cooma56424_test, series="Test Data",alpha=.8) +
theme_bw(base_size=20)
Cooma_yearly
accuracy(forecast(Cooma56424_tbats),Cooma56424_test)
##############################
### PACIFIC 57456 - NORTH ####
##############################
# Filter to required stations
Pacific57456 <- Station_data %>%
filter(station_key == 57456 ) %>%
filter(cardinal_direction_seq == 1)
# Clean
Pacific57456 <- Pacific57456 %>%
filter(daily_total > 500 & daily_total<6000) %>%
filter(date<"2015-06-28")
# Select the relevant variables
Pacific57456 <- Pacific57456 %>%
select(date,daily_total)
Pacific57456 <- Pacific57456[order(as.Date(Pacific57456$date, format="%Y-%m-%d")),]
# Add in NAs for missing dates
Pacific57456_dates <- seq.Date(as.Date(Pacific57456[1,1]),as.Date(Pacific57456[nrow(Pacific57456),1]),by=1)
Pacific57456_dates <-as.data.frame(Pacific57456_dates)
names(Pacific57456_dates)[names(Pacific57456_dates)=="Pacific57456_dates"] <- "date"
Pacific57456 <- left_join(Pacific57456_dates,Pacific57456,by="date")
# Calculate NAs
Pacific57456_rows <- nrow(Pacific57456)
Pacific57456_NAs <- sum(is.na(Pacific57456$daily_total))/Pacific57456_rows
# Show NAs
Pacific57456_NA <-zoo(Pacific57456$daily_total,Pacific57456$date)
plotNA.distribution(Pacific57456_NA)
# Address NAs
Pacific57456$year <- format(as.Date(Pacific57456$date, format="%d/%m/%Y"),"%Y")
Pacific57456$day <- weekdays(as.Date(Pacific57456$date))
Pacific57456 <- Pacific57456 %>%
group_by(year,day) %>%
mutate(average = round(mean(daily_total,na.rm = TRUE),0))
Pacific57456$daily_total[is.na(Pacific57456$daily_total)] <- Pacific57456$average[is.na(Pacific57456$daily_total)]
# Aggregate to Weekly
Pacific57456_zoo <- zoo(Pacific57456$daily_total,Pacific57456$date)
Pacific57456_week <- apply.weekly(Pacific57456_zoo, mean)
Pacific57456_ts <- msts(Pacific57456_week,seasonal.periods=52.1785,start=2008+259/365.25)
Pacific57456_tbats_full <- tbats(Pacific57456_ts)
plot(forecast(Pacific57456_tbats_full))
plot(forecast(Pacific57456_tbats_full,h=104),include=104)
# Visualise Weekly data
Pacific57456_plot <- window(Pacific57456_ts,start=2009,end=2016)
ggseasonplot(Pacific57456_plot)
# Decompose
plot(Pacific57456_tbats_full)
Pacific57456_season_tbats <- tbats.components(Pacific57456_tbats_full)[,'season']
autoplot(Pacific57456_season_tbats)
# or
Pacific57456_decom <- decompose(Pacific57456_ts, type="multiplicative")
autoplot(Pacific57456_decom)
Pacific57456_season_decom <- Pacific57456_decom$seasonal
autoplot(Pacific57456_season_decom)
# Model and test
Pacific57456_train <- window(Pacific57456_ts, end=2014+180/365.25)
Pacific57456_test <- window(Pacific57456_ts, start=2014+180/365.25)
Pacific57456_tbats <- tbats(Pacific57456_train)
autoplot(forecast(Pacific57456_tbats,h=50),ylab="Hornsby Count") +
autolayer(Pacific57456_test, series="Test Data",alpha=.7) +
theme_bw(base_size=20)
autoplot(forecast(Pacific57456_tbats,h=52),include=104) +
autolayer(Pacific57456_test, series="Test Data",alpha=.7)
accuracy(forecast(Pacific57456_tbats),Pacific57456_test)
# Test for ACF and PACF
Pacific57456_tbats_random <- (tbats.components(Pacific57456_tbats)[,'observed'])-(tbats.components(Pacific57456_tbats)[,'level'])-(tbats.components(Pacific57456_tbats)[,'season'])
plot(Pacific57456_tbats_random)
Acf(Pacific57456_tbats_random,lag.max = 20,main="Hornsby ACF")
Pacf(Pacific57456_tbats_random,lag.max = 20,main="Hornsby PACF")
##############################
### PORT 99990008 - NORTH ####
##############################
# Filter to required stations
Port99990008 <- Station_data %>%
filter(station_key == 99990008 ) %>%
filter(cardinal_direction_seq == 1)
# Clean
Port99990008 <- Port99990008 %>%
group_by(date) %>%
summarise(daily_total = sum(daily_total)) %>%
ungroup()
# Select the relevant variables
Port99990008 <- Port99990008 %>%
select(date,daily_total)
Port99990008 <- Port99990008[order(as.Date(Port99990008$date, format="%Y-%m-%d")),]
# Add in NAs for missing dates
Port99990008_dates <- seq.Date(as.Date("2015-04-01"),as.Date("2017-04-29"),by=1)
Port99990008_dates <-as.data.frame(Port99990008_dates)
names(Port99990008_dates)[names(Port99990008_dates)=="Port99990008_dates"] <- "date"
Port99990008 <- left_join(Port99990008_dates,Port99990008,by="date")
# Calculate NAs
Port99990008_rows <- nrow(Port99990008)
Port99990008_NAs <- sum(is.na(Port99990008$daily_total))/Port99990008_rows
# Show NAs
Port99990008_NA <-zoo(Port99990008$daily_total,Port99990008$date)
plotNA.distribution(Port99990008_NA)
# Address NAs
Port99990008$year <- format(as.Date(Port99990008$date, format="%d/%m/%Y"),"%Y")
Port99990008$day <- weekdays(as.Date(Port99990008$date))
Port99990008 <- Port99990008 %>%
group_by(year,day) %>%
mutate(average = round(mean(daily_total,na.rm = TRUE),0))
Port99990008$daily_total[is.na(Port99990008$daily_total)] <- Port99990008$average[is.na(Port99990008$daily_total)]
# Aggregate to Weekly
Port99990008_zoo <- zoo(Port99990008$daily_total,Port99990008$date)
Port99990008_week <- apply.weekly(Port99990008_zoo, mean)
Port99990008_ts <- msts(Port99990008_week,seasonal.periods=52.1785,start=2015+90/365.25)
Port99990008_tbats_full <- tbats(Port99990008_ts)
plot(forecast(Port99990008_tbats_full))
plot(forecast(Port99990008_tbats_full,h=104),include=104)
# Visualise Weekly data
Port99990008_plot <- window(Port99990008_ts,start=2015+90/365.25,end=2017)
ggseasonplot(Port99990008_plot)
# Decompose
plot(Port99990008_tbats_full)
Port99990008_season_tbats <- tbats.components(Port99990008_tbats_full)[,'season']
autoplot(Port99990008_season_tbats)
# or
Port99990008_decom <- decompose(Port99990008_ts, type="multiplicative")
autoplot(Port99990008_decom)
Port99990008_season_decom <- Port99990008_decom$seasonal
autoplot(Port99990008_season_decom)
# Model and test
Port99990008_train <- window(Port99990008_ts, end=2017+91/365.25)
Port99990008_test <- window(Port99990008_ts, start=2017+91/365.25)
Port99990008_tbats <- tbats(Port99990008_train)
autoplot(forecast(Port99990008_tbats,h=26),ylab="Port Macquarie Count") +
autolayer(Port99990008_test, series="Test Data",alpha=.8) +
theme_bw(base_size=20)
accuracy(forecast(Port99990008_tbats),Port99990008_test)
######################
### BYRON 15190087 ###
######################
# Filter to required stations
Byron15190087 <- Station_data %>%
filter(station_key == 15190087 ) %>%
filter(cardinal_direction_seq == 1) #%>%
# Clean
Byron15190087 <- Byron15190087 %>%
group_by(date) %>%
summarise(daily_total = sum(daily_total)) %>%
ungroup()
# Select the relevant variables
Byron15190087 <- Byron15190087 %>%
select(date,daily_total)
Byron15190087 <- Byron15190087[order(as.Date(Byron15190087$date, format="%Y-%m-%d")),]
# Add in NAs for missing dates
Byron15190087_dates <- seq.Date(as.Date("2015-04-01"),as.Date("2017-04-29"),by=1)
Byron15190087_dates <-as.data.frame(Byron15190087_dates)
names(Byron15190087_dates)[names(Byron15190087_dates)=="Byron15190087_dates"] <- "date"
Byron15190087 <- left_join(Byron15190087_dates,Byron15190087,by="date")
# Calculate NAs
Byron15190087_rows <- nrow(Byron15190087)
Byron15190087_NAs <- sum(is.na(Byron15190087$daily_total))/Byron15190087_rows
# Show NAs
Byron15190087_NA <-zoo(Byron15190087$daily_total,Byron15190087$date)
plotNA.distribution(Byron15190087_NA)
# Address NAs
Byron15190087$year <- format(as.Date(Byron15190087$date, format="%d/%m/%Y"),"%Y")
Byron15190087$day <- weekdays(as.Date(Byron15190087$date))
Byron15190087 <- Byron15190087 %>%
group_by(year,day) %>%
mutate(average = round(mean(daily_total,na.rm = TRUE),0))
Byron15190087$daily_total[is.na(Byron15190087$daily_total)] <- Byron15190087$average[is.na(Byron15190087$daily_total)]
# Aggregate to Weekly
Byron15190087_zoo <- zoo(Byron15190087$daily_total,Byron15190087$date)
Byron15190087_week <- apply.weekly(Byron15190087_zoo, mean)
Byron15190087_ts <- msts(Byron15190087_week,seasonal.periods=52.1785,start=2015+90/365.25)
Byron15190087_tbats_full <- tbats(Byron15190087_ts)
plot(forecast(Byron15190087_tbats_full))
plot(forecast(Byron15190087_tbats_full,h=104),include=104)
# Visualise Weekly data
Byron15190087_plot <- window(Byron15190087_ts,start=2015+90/365.25,end=2017)
ggseasonplot(Byron15190087_plot)
# Decompose
plot(Byron15190087_tbats_full)
Byron15190087_season_tbats <- tbats.components(Byron15190087_tbats_full)[,'season']
autoplot(Byron15190087_season_tbats)
# or
Byron15190087_decom <- decompose(Byron15190087_ts, type="multiplicative")
autoplot(Byron15190087_decom)
Byron15190087_season_decom <- Byron15190087_decom$seasonal
autoplot(Byron15190087_season_decom)
# Model and test
Byron15190087_train <- window(Byron15190087_ts, end=2017+91/365.25)
Byron15190087_test <- window(Byron15190087_ts, start=2017+91/365.25)
Byron15190087_tbats <- tbats(Byron15190087_train)
autoplot(forecast(Byron15190087_tbats,h=26),ylab="Byron Bay Count") +
autolayer(Byron15190087_test, series="Test Data",alpha=.8) +
theme_bw(base_size=20)
accuracy(forecast(Byron15190087_tbats),Byron15190087_test)
#####################################
### BLUE MOUNTAINS 58847 - WEST ####
#####################################
# Filter to required stations
BlueMt58847 <- Station_data %>%
filter(station_key == 58847 ) %>%
filter(cardinal_direction_seq == 7)
# Clean
BlueMt58847 <- BlueMt58847 %>%
filter(date<"2013-04-01")
# Select the relevant variables
BlueMt58847 <- BlueMt58847 %>%
select(date,daily_total)
BlueMt58847 <- BlueMt58847[order(as.Date(BlueMt58847$date, format="%Y-%m-%d")),]
# Add in NAs for missing dates
BlueMt58847_dates <- seq.Date(as.Date(BlueMt58847[1,1]),as.Date(BlueMt58847[nrow(BlueMt58847),1]),by=1)
BlueMt58847_dates <-as.data.frame(BlueMt58847_dates)
names(BlueMt58847_dates)[names(BlueMt58847_dates)=="BlueMt58847_dates"] <- "date"
BlueMt58847 <- left_join(BlueMt58847_dates,BlueMt58847,by="date")
# Calculate NAs
BlueMt58847_rows <- nrow(BlueMt58847)
BlueMt58847_NAs <- sum(is.na(BlueMt58847$daily_total))/BlueMt58847_rows
# Show NAs
BlueMt58847_NA <-zoo(BlueMt58847$daily_total,BlueMt58847$date)
plotNA.distribution(BlueMt58847_NA)
# Address NAs
BlueMt58847$year <- format(as.Date(BlueMt58847$date, format="%d/%m/%Y"),"%Y")
BlueMt58847$day <- weekdays(as.Date(BlueMt58847$date))
BlueMt58847 <- BlueMt58847 %>%
group_by(year,day) %>%
mutate(average = round(mean(daily_total,na.rm = TRUE),0))
BlueMt58847$daily_total[is.na(BlueMt58847$daily_total)] <- BlueMt58847$average[is.na(BlueMt58847$daily_total)]
# Aggregate to Weekly
BlueMt58847_zoo <- zoo(BlueMt58847$daily_total,BlueMt58847$date)
BlueMt58847_week <- apply.weekly(BlueMt58847_zoo, mean)
BlueMt58847_ts <- msts(BlueMt58847_week,seasonal.periods=52.1785,start=2008+268/365.25)
BlueMt58847_tbats_full <- tbats(BlueMt58847_ts)
plot(forecast(BlueMt58847_tbats_full))
plot(forecast(BlueMt58847_tbats_full,h=104),include=104)
# Visualise Weekly data
BlueMt58847_plot <- window(BlueMt58847_ts,start=2009,end=2013)
ggseasonplot(BlueMt58847_plot)
# Decompose
plot(BlueMt58847_tbats_full)
BlueMt58847_season_tbats <- tbats.components(BlueMt58847_tbats_full)[,'season']
autoplot(BlueMt58847_season_tbats)
# or
BlueMt58847_decom <- decompose(BlueMt58847_ts, type="multiplicative")
autoplot(BlueMt58847_decom)
BlueMt58847_season_decom <- BlueMt58847_decom$seasonal
autoplot(BlueMt58847_season_decom)
# Model and test
BlueMt58847_train <- window(BlueMt58847_ts, end=2012+100/365.25)
BlueMt58847_test <- window(BlueMt58847_ts, start=2012+100/365.25)
BlueMt58847_tbats <- tbats(BlueMt58847_train)
autoplot(forecast(BlueMt58847_tbats,h=51),ylab="Blue Mountains Count") +
autolayer(BlueMt58847_test, series="Test Data",alpha=.8)+
theme_bw(base_size=20)
accuracy(forecast(BlueMt58847_tbats),BlueMt58847_test)
###############################
### BATHURST 58691 - WEST ####
###############################
# Filter to required stations
Bathurst58691 <- Station_data %>%
filter(station_key == 58691 ) %>%
filter(cardinal_direction_seq == 7)
# Clean
Bathurst58691 <- Bathurst58691 %>%
filter(date<"2015-08-01")
# Select the relevant variables
Bathurst58691 <- Bathurst58691 %>%
select(date,daily_total)
Bathurst58691 <- Bathurst58691[order(as.Date(Bathurst58691$date, format="%Y-%m-%d")),]
# Add in NAs for missing dates
Bathurst58691_dates <- seq.Date(as.Date(Bathurst58691[1,1]),as.Date(Bathurst58691[nrow(Bathurst58691),1]),by=1)
Bathurst58691_dates <-as.data.frame(Bathurst58691_dates)
names(Bathurst58691_dates)[names(Bathurst58691_dates)=="Bathurst58691_dates"] <- "date"
Bathurst58691 <- left_join(Bathurst58691_dates,Bathurst58691,by="date")
# Calculate NAs
Bathurst58691_rows <- nrow(Bathurst58691)
Bathurst58691_NAs <- sum(is.na(Bathurst58691$daily_total))/Bathurst58691_rows
# Show NAs
Bathurst58691_NA <-zoo(Bathurst58691$daily_total,Bathurst58691$date)
plotNA.distribution(Bathurst58691_NA,main="Distribution of NAs on Great Western Highway near Bathurst")
# Address NAs
Bathurst58691$year <- format(as.Date(Bathurst58691$date, format="%d/%m/%Y"),"%Y")
Bathurst58691$day <- weekdays(as.Date(Bathurst58691$date))
Bathurst58691 <- Bathurst58691 %>%
group_by(year,day) %>%
mutate(average = round(mean(daily_total,na.rm = TRUE),0))
Bathurst58691$daily_total[is.na(Bathurst58691$daily_total)] <- Bathurst58691$average[is.na(Bathurst58691$daily_total)]
# Aggregate to Weekly
Bathurst58691_zoo <- zoo(Bathurst58691$daily_total,Bathurst58691$date)
Bathurst58691_week <- apply.weekly(Bathurst58691_zoo, mean)
Bathurst58691_ts <- msts(Bathurst58691_week,seasonal.periods=52.1785,start=2008+289/365.25)
Bathurst58691_tbats_full <- tbats(Bathurst58691_ts)
plot(forecast(Bathurst58691_tbats_full))
plot(forecast(Bathurst58691_tbats_full,h=104),include=104)
# Visualise Weekly data
Bathurst58691_plot <- window(Bathurst58691_ts,start=2009,end=2015)
ggseasonplot(Bathurst58691_plot)
# Decompose
plot(Bathurst58691_tbats_full)
Bathurst58691_season_tbats <- tbats.components(Bathurst58691_tbats_full)[,'season']
autoplot(Bathurst58691_season_tbats)
# or
Bathurst58691_decom <- decompose(Bathurst58691_ts, type="multiplicative")
autoplot(Bathurst58691_decom)
Bathurst58691_season_decom <- Bathurst58691_decom$seasonal
autoplot(Bathurst58691_season_decom)
# Model and test
Bathurst58691_train <- window(Bathurst58691_ts, end=2014+170/365.25)
Bathurst58691_test <- window(Bathurst58691_ts, start=2014+170/365.25)
Bathurst58691_tbats <- tbats(Bathurst58691_train)
autoplot(forecast(Bathurst58691_tbats,h=51),ylab="Bathurst Count") +
autolayer(Bathurst58691_test, series="Test Data",alpha=.8)+
theme_bw(base_size=20)
accuracy(forecast(Bathurst58691_tbats),Bathurst58691_test)
# test for ACF and PACF
Bathurst58691_tbats_random <- (tbats.components(Bathurst58691_tbats)[,'observed'])-(tbats.components(Bathurst58691_tbats)[,'level'])-(tbats.components(Bathurst58691_tbats)[,'season'])
plot(Bathurst58691_tbats_random)
Acf(Bathurst58691_tbats_random,lag.max = 20,main="Bathurst ACF")
Pacf(Bathurst58691_tbats_random,lag.max = 20,main="Bathurst PACF")
###########################
### YASS 58049 - WEST ####
###########################
# Filter to required stations
Yass58049 <- Station_data %>%
filter(station_key == 58049 ) %>%
filter(cardinal_direction_seq == 7)
# Clean
Yass58049 <- Yass58049 %>%
filter(date<"2016-05-01" & date>"2010-05-01")
# Select the relevant variables
Yass58049 <- Yass58049 %>%
select(date,daily_total)
Yass58049 <- Yass58049[order(as.Date(Yass58049$date, format="%Y-%m-%d")),]
# Add in NAs for missing dates
Yass58049_dates <- seq.Date(as.Date(Yass58049[1,1]),as.Date(Yass58049[nrow(Yass58049),1]),by=1)
Yass58049_dates <-as.data.frame(Yass58049_dates)
names(Yass58049_dates)[names(Yass58049_dates)=="Yass58049_dates"] <- "date"
Yass58049 <- left_join(Yass58049_dates,Yass58049,by="date")
# Calculate NAs
Yass58049_rows <- nrow(Yass58049)
Yass58049_NAs <- sum(is.na(Yass58049$daily_total))/Yass58049_rows
# Show NAs
Yass58049_NA <-zoo(Yass58049$daily_total,Yass58049$date)
plotNA.distribution(Yass58049_NA)
# Address NAs
Yass58049$year <- format(as.Date(Yass58049$date, format="%d/%m/%Y"),"%Y")
Yass58049$day <- weekdays(as.Date(Yass58049$date))
Yass58049 <- Yass58049 %>%
group_by(year,day) %>%
mutate(average = round(mean(daily_total,na.rm = TRUE),0))
Yass58049$daily_total[is.na(Yass58049$daily_total)] <- Yass58049$average[is.na(Yass58049$daily_total)]
# Aggregate to Weekly
Yass58049_zoo <- zoo(Yass58049$daily_total,Yass58049$date)
Yass58049_week <- apply.weekly(Yass58049_zoo, mean)
Yass58049_ts <- msts(Yass58049_week,seasonal.periods=52.1785,start=2010+250/365.25)
Yass58049_tbats_full <- tbats(Yass58049_ts)
plot(forecast(Yass58049_tbats_full))
plot(forecast(Yass58049_tbats_full,h=104),include=104)
# Visualise Weekly data
Yass58049_plot <- window(Yass58049_ts,start=2011,end=2016)
ggseasonplot(Yass58049_plot)
# Decompose
plot(Yass58049_tbats_full)
Yass58049_season_tbats <- tbats.components(Yass58049_tbats_full)[,'season']
autoplot(Yass58049_season_tbats)
# or
Yass58049_decom <- decompose(Yass58049_ts, type="multiplicative")
autoplot(Yass58049_decom)
Yass58049_season_decom <- Yass58049_decom$seasonal
autoplot(Yass58049_season_decom)
# Model and test
Yass58049_train <- window(Yass58049_ts, end=2015+180/365.25)
Yass58049_test <- window(Yass58049_ts, start=2015+180/365.25)
Yass58049_tbats <- tbats(Yass58049_train)
autoplot(forecast(Yass58049_tbats,h=45),ylab="Yass Count") +
autolayer(Yass58049_test, series="Test Data",alpha=.8) +
theme_bw(base_size=20)
accuracy(forecast(Yass58049_tbats),Yass58049_test)
################
#### GRAPHS ####
################
# Map of NSW
stations_map <- stations %>%
filter(station_key %in% Station_code) %>%
select(station_key,wgs84_latitude,wgs84_longitude)
NSW <- ggmap(get_googlemap(c(lon=147.5, lat=-32.8), zoom=6, maptype="roadmap", size=c(600,500)))
NSW +
geom_point(aes(x=wgs84_longitude, y=wgs84_latitude), data=stations_map,colour="red",size=4, alpha=0.8) +
labs(x="Longitude",y="Latitiude") +
ggtitle("Traffic Counter Stations") +
#theme(plot.title = element_text(family = "Trebuchet MS", color="#666666", face="bold", size=22, hjust=0)) +
#theme(axis.title = element_text(family = "Trebuchet MS", color="#666666", face="bold", size=12)) +
theme_bw(base_size = 20)
# Calculate distance to Sydney
get.distance <- function(x,y,z) {
results = gmapsdistance(
origin = paste(x,y,sep=","),
destination = "-33.8704512,151.2058792",
mode = "driving",
traffic_model = "best_guess",
shape = "long",
key = "AIzaSyCwoq2tqj0RX54SrVnnRYccCs6AV1kHj34")
results.df<-data.frame(results)
results.df$station_key <- z
results.df
}
stations_table <- stations %>%
filter(station_key %in% Station_code) %>%
select(station_key,wgs84_latitude,wgs84_longitude,road_name,lga,rms_region)
# Calculate distance for each station
stations_distance <- get.distance(stations_table$wgs84_latitude,stations_table$wgs84_longitude, stations_table$station_key)
# Establish dataframe of distances and merge to stations
stations_distance <- select(stations_distance,station_key,Distance.Distance)
stations_table <- inner_join(stations_table,stations_distance,by="station_key")
stations_table <- stations_table %>%
mutate(Distance_Sydney_km=Distance.Distance/1000) %>%
select(road_name,lga,rms_region,Distance_Sydney_km)
# Combine seasonal charts
Usage_Patterns <- ts.union(Princes56914_season_decom,
Batemans56413_season_decom,
Cooma56424_season_decom,
Pacific57456_season_decom,
Bathurst58691_season_decom,
BlueMt58847_season_decom,
Yass58049_season_decom,
Port99990008_season_decom,
Byron15190087_season_decom,
dframe = FALSE)
colnames(Usage_Patterns) <- c("Heathcote","Batemans", "Cooma", "Hornsby","Bathurst","Blue Mountains","Yass","Port Macquarie","Byron Bay")
autoplot(Usage_Patterns,facets=TRUE)+
scale_x_continuous(limits= c(2006,2018),breaks=c(2006,2007,2008,2009,2010,2011,2012,2013,2014,2015,2016,2017,2018)) +
theme_bw(base_size = 25)
#### USE TO SHOW ALL THE OF THE FORECASTS
plot1 <- autoplot(forecast(Princes56914_tbats,h=45),ylab="Heathcote Count") +
autolayer(Princes56914_test, series="Test Data",alpha=.9) +
theme_bw(base_size=30)
plot2 <- autoplot(forecast(Batemans56413_tbats,h=45),ylab="Batemans Bay Count") +
autolayer(Batemans56413_test, series="Test Data",alpha=.9) +
theme_bw(base_size=30)
plot3 <- autoplot(forecast(Cooma56424_tbats,h=45),ylab="Cooma Count") +
autolayer(Cooma56424_test, series="Test Data",alpha=.9) +
theme_bw(base_size=30)
plot4 <- autoplot(forecast(Pacific57456_tbats,h=50),ylab="Hornsby Count") +
autolayer(Pacific57456_test, series="Test Data",alpha=.9) +
theme_bw(base_size=30)
plot5 <- autoplot(forecast(Port99990008_tbats,h=50),ylab="Port Macquarie Count") +
autolayer(Port99990008_test, series="Test Data",alpha=.9) +
theme_bw(base_size=30)
plot6 <- autoplot(forecast(Byron15190087_tbats,h=50),ylab="Byron Bay Count") +
autolayer(Byron15190087_test, series="Test Data",alpha=.9) +
theme_bw(base_size=30)
plot7 <- autoplot(forecast(BlueMt58847_tbats,h=51),ylab="Blue Mountains Count") +
autolayer(BlueMt58847_test, series="Test Data",alpha=.9)+
theme_bw(base_size=30)
plot8 <- autoplot(forecast(Bathurst58691_tbats,h=51),ylab="Bathurst Count") +
autolayer(Bathurst58691_test, series="Test Data",alpha=.9)+
theme_bw(base_size=30)
plot9 <- autoplot(forecast(Yass58049_tbats,h=45),ylab="Yass Count") +
autolayer(Yass58049_test, series="Test Data",alpha=.9) +
theme_bw(base_size=30)
grid.arrange(plot1, plot2, plot3, nrow=3)
grid.arrange(plot4, plot5, plot6,nrow=3)
grid.arrange(plot7, plot8, plot9,nrow=3)
#### USE TO SHOW SMALL AMOUNT OF DATA & FORECAST 2 YEARS WITH PARTIAL DATA COMPARISON
plot10 <- autoplot(forecast(Princes56914_tbats,h=104),ylab="Heathcote Count",include=13) +
autolayer(Princes56914_test, series="Test Data",alpha=.9) +
theme_bw(base_size=30)
plot11 <- autoplot(forecast(Batemans56413_tbats,h=104),ylab="Batemans Bay Count",include=13) +
autolayer(Batemans56413_test, series="Test Data",alpha=.9) +
theme_bw(base_size=30)
plot12 <- autoplot(forecast(Cooma56424_tbats,h=104),ylab="Cooma Count",include=13) +
autolayer(Cooma56424_test, series="Test Data",alpha=.9) +
theme_bw(base_size=30)
plot13 <- autoplot(forecast(Pacific57456_tbats,h=104),ylab="Hornsby Count",include=13) +
autolayer(Pacific57456_test, series="Test Data",alpha=.9) +
theme_bw(base_size=30)
plot14 <- autoplot(forecast(Port99990008_tbats,h=104),ylab="Port Macquarie Count",include=13) +
autolayer(Port99990008_test, series="Test Data",alpha=.9) +
theme_bw(base_size=30)
plot15 <- autoplot(forecast(Byron15190087_tbats,h=104),ylab="Byron Bay Count",include=13) +
autolayer(Byron15190087_test, series="Test Data",alpha=.9) +
theme_bw(base_size=30)
plot16 <- autoplot(forecast(BlueMt58847_tbats,h=104),ylab="Blue Mountains Count",include=13) +
autolayer(BlueMt58847_test, series="Test Data",alpha=.9)+
theme_bw(base_size=30)
plot17 <- autoplot(forecast(Bathurst58691_tbats,h=104),ylab="Bathurst Count",include=13) +
autolayer(Bathurst58691_test, series="Test Data",alpha=.9)+
theme_bw(base_size=30)
plot18 <- autoplot(forecast(Yass58049_tbats,h=104),ylab="Yass Count",include=13) +
autolayer(Yass58049_test, series="Test Data",alpha=.9) +
theme_bw(base_size=30)
grid.arrange(plot10, plot11, plot12, nrow=3)
grid.arrange(plot13, plot14, plot15,nrow=3)
grid.arrange(plot16, plot17, plot18,nrow=3)
plot11
# Graphing comparison in Cooma
grid.arrange(Cooma_weekly, Cooma_yearly, nrow=2)
# Graphing scaled estimate of seasonality from TBATS models
Princes56914_tbats_yearseason <- tbats.components(Princes56914_tbats)[,'season']
Princes56914_tbats_yearseason <- window(Princes56914_tbats_yearseason,start=2015,end=2015+365/365.25)
Princes56914_tbats_yearseason <- scale(Princes56914_tbats_yearseason)
Princes56914_tbats_yearseason <- msts(Princes56914_tbats_yearseason,seasonal.periods=52.1785,start=0)
Batemans56413_tbats_yearseason <- tbats.components(Batemans56413_tbats)[,'season']
Batemans56413_tbats_yearseason <- window(Batemans56413_tbats_yearseason,start=2015,end=2015+365/365.25)
Batemans56413_tbats_yearseason <- scale(Batemans56413_tbats_yearseason)
Batemans56413_tbats_yearseason <- msts(Batemans56413_tbats_yearseason,seasonal.periods=52.1785,start=0)
Cooma56424_tbats_yearseason <- tbats.components(Cooma56424_tbats)[,'season']
Cooma56424_tbats_yearseason <- window(Cooma56424_tbats_yearseason,start=2015,end=2015+365/365.25)
Cooma56424_tbats_yearseason <- scale(Cooma56424_tbats_yearseason)
Cooma56424_tbats_yearseason <- msts(Cooma56424_tbats_yearseason,seasonal.periods=52.1785,start=0)
Pacific57456_tbats_yearseason <- tbats.components(Pacific57456_tbats)[,'season']
Pacific57456_tbats_yearseason <- window(Pacific57456_tbats_yearseason,start=2012,end=2012+365/365.25)
Pacific57456_tbats_yearseason <- scale(Pacific57456_tbats_yearseason)
Pacific57456_tbats_yearseason <- msts(Pacific57456_tbats_yearseason,seasonal.periods=52.1785,start=0)
Port99990008_tbats_yearseason <- tbats.components(Port99990008_tbats)[,'season']
Port99990008_tbats_yearseason1 <- window(Port99990008_tbats_yearseason,start=2015+90/365.25,end=2015+365/365.25)
Port99990008_tbats_yearseason2 <- window(Port99990008_tbats_yearseason,start=2016,end=2016+90/365.25)
Port99990008_tbats_yearseason <- c(Port99990008_tbats_yearseason1,Port99990008_tbats_yearseason2)
Port99990008_tbats_yearseason <- scale(Port99990008_tbats_yearseason)
Port99990008_tbats_yearseason <- msts(Port99990008_tbats_yearseason,seasonal.periods=52.1785,start=0)
Byron15190087_tbats_yearseason <- tbats.components(Byron15190087_tbats)[,'season']
Byron15190087_tbats_yearseason1 <- window(Byron15190087_tbats_yearseason,start=2015+90/365.25,end=2015+365/365.25)
Byron15190087_tbats_yearseason2 <- window(Byron15190087_tbats_yearseason,start=2016,end=2016+90/365.25)
Byron15190087_tbats_yearseason <- c(Byron15190087_tbats_yearseason1,Byron15190087_tbats_yearseason2)
Byron15190087_tbats_yearseason <- scale(Byron15190087_tbats_yearseason)
Byron15190087_tbats_yearseason <- msts(Byron15190087_tbats_yearseason,seasonal.periods=52.1785,start=0)
BlueMt58847_tbats_yearseason <- tbats.components(BlueMt58847_tbats)[,'season']
BlueMt58847_tbats_yearseason <- window(BlueMt58847_tbats_yearseason,start=2011,end=2011+365/365.25)
BlueMt58847_tbats_yearseason <- scale(BlueMt58847_tbats_yearseason)
BlueMt58847_tbats_yearseason <- msts(BlueMt58847_tbats_yearseason,seasonal.periods=52.1785,start=0)
Bathurst58691_tbats_yearseason <- tbats.components(Bathurst58691_tbats)[,'season']
Bathurst58691_tbats_yearseason <- window(Bathurst58691_tbats_yearseason,start=2013,end=2013+365/365.25)
Bathurst58691_tbats_yearseason <- scale(Bathurst58691_tbats_yearseason)
Bathurst58691_tbats_yearseason <- msts(Bathurst58691_tbats_yearseason,seasonal.periods=52.1785,start=0)
Yass58049_tbats_yearseason <- tbats.components(Yass58049_tbats)[,'season']
Yass58049_tbats_yearseason <- window(Yass58049_tbats_yearseason,start=2013,end=2013+365/365.25)
Yass58049_tbats_yearseason <- scale(Yass58049_tbats_yearseason)
Yass58049_tbats_yearseason <- msts(Yass58049_tbats_yearseason,seasonal.periods=52.1785,start=0)
Monthbreaks <- c(0/365,31/365,59/365,90/365,120/365,151/365,181/365,212/365,243/365,273/365,304/365,334/365,365/365)
Monthlabels <- c("1-Jan","1-Feb","1-Mar","1-Apr","1-May","1-Jun","1-Jul","1-Aug","1-Sept","1-Oct","1-Nov","1-Dec","31-Dec")
# All
autoplot(Princes56914_tbats_yearseason,ylab="Scaled Count") +
autolayer(Princes56914_tbats_yearseason,series = "Heathcote") +
autolayer(Batemans56413_tbats_yearseason,series = "Batemans Bay") +
autolayer(Cooma56424_tbats_yearseason,series = "Cooma") +
autolayer(Pacific57456_tbats_yearseason,series = "Hornsby") +
autolayer(Port99990008_tbats_yearseason,series = "Port Macquarie") +
autolayer(Byron15190087_tbats_yearseason,series = "Byron Bay") +
autolayer(BlueMt58847_tbats_yearseason,series = "Blue Mountains") +
autolayer(Bathurst58691_tbats_yearseason,series = "Bathurst") +
autolayer(Yass58049_tbats_yearseason,series = "Yass") +
scale_x_continuous(breaks=Monthbreaks,labels=Monthlabels,minor_breaks = NULL) +
theme_bw(base_size = 20)
## Travelling South
autoplot(Princes56914_tbats_yearseason,ylab="Scaled Count",main="Travelling South") +
autolayer(Cooma56424_tbats_yearseason,series = "Cooma") +
autolayer(Batemans56413_tbats_yearseason,series = "Batemans Bay") +
autolayer(Princes56914_tbats_yearseason,series = "Heathcote") +
scale_x_continuous(breaks=Monthbreaks,labels=Monthlabels,minor_breaks = NULL) +
theme_bw(base_size = 20)
## Travelling North
autoplot(Pacific57456_tbats_yearseason,ylab="Scaled Count",main="Travelling North") +
autolayer(Byron15190087_tbats_yearseason,series = "Byron Bay") +
autolayer(Port99990008_tbats_yearseason,series = "Port Macquarie") +
autolayer(Pacific57456_tbats_yearseason,series = "Hornsby") +
scale_x_continuous(breaks=Monthbreaks,labels=Monthlabels,minor_breaks = NULL) +
theme_bw(base_size = 20)
## Travelling West
autoplot(Yass58049_tbats_yearseason,ylab="Scaled Count",main="Travelling West") +
autolayer(Yass58049_tbats_yearseason,series = "Yass") +
autolayer(Bathurst58691_tbats_yearseason,series = "Bathurst") +
autolayer(BlueMt58847_tbats_yearseason,series = "Blue Mountains") +
scale_x_continuous(breaks=Monthbreaks,labels=Monthlabels,minor_breaks = NULL) +
theme_bw(base_size = 20)
### Travelling West -> Easter and April School Holidays
autoplot(Bathurst58691_tbats_yearseason,ylab="Scaled Count",main="Group 2 - Using Easter & April School Holidays to Travel West") +
autolayer(BlueMt58847_tbats_yearseason,series = "Blue Mountains") +
autolayer(Bathurst58691_tbats_yearseason,series = "Bathurst") +
scale_x_continuous(breaks=Monthbreaks,labels=Monthlabels,minor_breaks = NULL) +
theme_bw(base_size = 20)
## Winter polarises
autoplot(Byron15190087_tbats_yearseason,ylab="Scaled Count",main="Group 3 - Winter Related Travel") +
autolayer(Byron15190087_tbats_yearseason,series = "Byron Bay") +
autolayer(Port99990008_tbats_yearseason,series = "Port Macquarie") +
autolayer(Cooma56424_tbats_yearseason,series = "Cooma") +
scale_x_continuous(breaks=Monthbreaks,labels=Monthlabels,minor_breaks = NULL) +
theme_bw(base_size = 20)
## North or West for Sept/Oct School holidays - Ignored
autoplot(Port99990008_tbats_yearseason,ylab="Scaled Count",main="Travelling North or West for Sept/Oct School Holidays") +
autolayer(Port99990008_tbats_yearseason,series = "Port Macquarie") +
autolayer(Byron15190087_tbats_yearseason,series = "Byron Bay") +
autolayer(Bathurst58691_tbats_yearseason,series = "Bathurst") +
autolayer(Yass58049_tbats_yearseason,series = "Yass") +
scale_x_continuous(breaks=Monthbreaks,labels=Monthlabels,minor_breaks = NULL) +
theme_bw(base_size = 20)
## Before Christmas
autoplot(Princes56914_tbats_yearseason,ylab="Scaled Count",main="Group 4 - Using December Holidays to Escape Sydney") +
autolayer(Pacific57456_tbats_yearseason,series = "Hornsby") +
autolayer(Princes56914_tbats_yearseason,series = "Heathcote") +
scale_x_continuous(breaks=Monthbreaks,labels=Monthlabels,minor_breaks = NULL) +
theme_bw(base_size = 20)
## South After New Years
autoplot(Batemans56413_tbats_yearseason,ylab="Scaled Count",main="Group 1 - Using the January Holidays to Travel South") +
autolayer(Yass58049_tbats_yearseason,series = "Yass") +
autolayer(Batemans56413_tbats_yearseason,series = "Batemans Bay") +
scale_x_continuous(breaks=Monthbreaks,labels=Monthlabels,minor_breaks = NULL) +
theme_bw(base_size = 20)