|
|
#R code for creating ggplots of time series with smooth (GAM) and linear term
|
|
|
|
|
|
|
|
|
######################
|
|
|
#Import Libraries and set working directory
|
|
|
######################
|
|
|
library(zoo)
|
|
|
library(hydroTSM) #you need to install these packages first before you can load them here
|
|
|
library(lubridate)
|
|
|
library(mgcv)
|
|
|
library(ggplot2)
|
|
|
library(gridExtra)
|
|
|
library(scales)
|
|
|
options(scipen=999)
|
|
|
setwd("C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdoc/Projects/Paper#1/")
|
|
|
######################
|
|
|
|
|
|
######################
|
|
|
#Set inputs
|
|
|
######################
|
|
|
Case.Study <- "CASESTUDY2"
|
|
|
Estuary <- "HUNTER"
|
|
|
Climvar <- 'tasmean'
|
|
|
ggplotGAM.k <- 14
|
|
|
######################
|
|
|
|
|
|
|
|
|
######################
|
|
|
#Set input file paths
|
|
|
######################
|
|
|
AirT_CSV_Path <- "./Data/SILO/CASESTUDY2/HUNTER/SILO_climdata_HUNTER_Catchment_-32.162479_150.5335812.csv"
|
|
|
RivT_CSV_Path <- "./Data/River_Gauge_Data/HUNTER_210064_20180615/210064_Temp.csv"
|
|
|
SST_CSV_Path <- './Data/NARCLIM_Site_CSVs/CASESTUDY2/HUNTER_32.751_151.690/sstmean_NNRP_HUNTER_33.034_152.156_NARCliM_summary.csv'
|
|
|
######################
|
|
|
|
|
|
|
|
|
######################
|
|
|
#Analyse
|
|
|
######################
|
|
|
|
|
|
############tasmean
|
|
|
#Load a daily (no gaps) time series as a time serie baseline for other time series used here
|
|
|
AirT.df <- data.frame(read.csv(AirT_CSV_Path))
|
|
|
AirT.full.TS <- zoo((AirT.df$max_temp_Celsius + AirT.df$max_temp_Celsius)/2, order.by= as.Date(AirT.df[,"date"], format = "%Y-%m-%d")) #=daily time series of rainfall for creation of clean, daily TS of ET and Q
|
|
|
AirT.TS <- window(AirT.full.TS, start=as.Date("1990-01-01"), end=as.Date("2018-01-01"))
|
|
|
AirT.full.df <- data.frame(AirT.full.TS)
|
|
|
AirT.df <- data.frame(AirT.TS)
|
|
|
colnames(AirT.df) <- 'tasmean'
|
|
|
colnames(AirT.full.df) <- 'tasmean'
|
|
|
############
|
|
|
AirT.df$Julday1 <- seq(1,length(AirT.df[,1]),1)
|
|
|
linear.trend.model_EC_all <- lm(tasmean ~ Julday1, AirT.df)
|
|
|
AirT.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4]
|
|
|
AirT.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 356
|
|
|
############
|
|
|
AirT.full.df$Julday1 <- seq(1,length(AirT.full.df[,1]),1)
|
|
|
linear.trend.model_EC_all <- lm(tasmean ~ Julday1, AirT.full.df)
|
|
|
AirT.full.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4]
|
|
|
AirT.full.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 356
|
|
|
############tasmean
|
|
|
|
|
|
|
|
|
############River temp
|
|
|
#Load a daily (no gaps) time series as a time serie baseline for other time series used here
|
|
|
RivT.df <- data.frame(read.csv(RivT_CSV_Path))
|
|
|
RivT.full.TS <- zoo(RivT.df$Temp, order.by= as.Date(RivT.df[,"Date"], format = "%d/%m/%Y")) #=daily time series of rainfall for creation of clean, daily TS of ET and Q
|
|
|
RivT.TS <- window(RivT.full.TS, start=as.Date("1995-01-01"), end=as.Date("2018-01-01"))
|
|
|
RivT.full.df <- data.frame(RivT.TS) ### This is only done because
|
|
|
RivT.df <- data.frame(RivT.TS)
|
|
|
colnames(RivT.df) <- 'rivTmean'
|
|
|
colnames(RivT.full.df) <- 'rivTmean'
|
|
|
############
|
|
|
RivT.df$Julday1 <- seq(1,length(RivT.df[,1]),1)
|
|
|
linear.trend.model_EC_all <- lm(rivTmean ~ Julday1, RivT.df)
|
|
|
RivT.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4]
|
|
|
RivT.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 356
|
|
|
RivT.full.df$Julday1 <- seq(1,length(RivT.full.df[,1]),1)
|
|
|
linear.trend.model_EC_all <- lm(rivTmean ~ Julday1, RivT.full.df)
|
|
|
RivT.full.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4]
|
|
|
RivT.full.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 356
|
|
|
############River temp
|
|
|
|
|
|
############ SST
|
|
|
#Load a daily (no gaps) time series as a time serie baseline for other time series used here
|
|
|
SST.df <- data.frame(read.csv(SST_CSV_Path))
|
|
|
SST.full.TS <- zoo(SST.df$NNRP_R1_1950-273.15, order.by= as.Date(SST.df[,"X"], format = "%Y-%m-%d")) #=daily time series of rainfall for creation of clean, daily TS of ET and Q
|
|
|
SST.TS <- window(SST.full.TS, start=as.Date("1990-01-01"), end=as.Date("2018-01-01"))
|
|
|
SST.full.df <- data.frame(SST.full.TS)
|
|
|
SST.df <- data.frame(SST.TS)
|
|
|
str(SST.df)
|
|
|
colnames(SST.df) <- 'SSTmean'
|
|
|
colnames(SST.full.df) <- 'SSTmean'
|
|
|
############
|
|
|
SST.full.df$Julday1 <- seq(1,length(SST.full.df[,1]),1)
|
|
|
linear.trend.model_EC_all <- lm(SSTmean ~ Julday1, SST.full.df)
|
|
|
SST.full.pvalNCV_ECall <- summary(linear.trend.model_EC_all)$coefficients[2,4]
|
|
|
SST.full.lintrend <- summary(linear.trend.model_EC_all)$coefficients[2,1] * 356
|
|
|
SST.df$Julday1 <- seq(1,length(SST.df[,1]),1)
|
|
|
linear.trend.model_EC_all2 <- lm(SSTmean ~ Julday1, SST.df)
|
|
|
SST.pvalNCV_ECall <- summary(linear.trend.model_EC_all2)$coefficients[2,4]
|
|
|
SST.lintrend <- summary(linear.trend.model_EC_all2)$coefficients[2,1] * 356
|
|
|
############ SST
|
|
|
|
|
|
|
|
|
|
|
|
######################
|
|
|
#Plot
|
|
|
######################
|
|
|
|
|
|
##################################### Full Time Period
|
|
|
#Air temp Full period
|
|
|
p1air <- ggplot(AirT.full.df, aes(y=tasmean, x=index(AirT.full.TS))) + geom_line(alpha=0.5) +
|
|
|
ggtitle(paste(Estuary, " - Linear and smooth trend in catchment airT (SILO) lin trend was ",
|
|
|
round(AirT.full.lintrend,3), ' C<>/year with p=', round(AirT.full.pvalNCV_ECall,10), sep=" ")) +
|
|
|
theme(plot.title=element_text(face="bold", size=9)) +
|
|
|
geom_smooth(method='lm',fill="green", formula=y~x, colour="darkgreen", size = 0.5) +
|
|
|
stat_smooth(method=gam, formula=y~s(x, k=ggplotGAM.k), se=T, size=0.5) +
|
|
|
ylab("Air Temperature [C<>]") + xlab("Time")
|
|
|
|
|
|
#Riv temp Full period
|
|
|
p1riv <- ggplot(RivT.full.df, aes(y=rivTmean, x=index(RivT.TS))) + geom_line(alpha=0.5) +
|
|
|
ggtitle(paste(Estuary, " - Linear and smooth trend in river temperature (GAUGE) lin trend was ",
|
|
|
round(RivT.full.lintrend,3), ' C<>/year with p=', round(RivT.full.pvalNCV_ECall,10), sep=" ")) +
|
|
|
theme(plot.title=element_text(face="bold", size=9)) +
|
|
|
geom_smooth(method='lm',fill="green", formula=y~x, colour="darkgreen", size = 0.5) +
|
|
|
stat_smooth(method=gam, formula=y~s(x, k=ggplotGAM.k), se=T, size=0.5) +
|
|
|
ylab("River Temperature [C<>]") + xlab("Time")
|
|
|
|
|
|
#Sea temp Full period
|
|
|
p1sst <- ggplot(SST.full.df, aes(y=SSTmean, x=index(SST.full.TS))) + geom_line(alpha=0.5) +
|
|
|
ggtitle(paste(Estuary, " - Linear and smooth trend in sea surface temperature (NNRP) lin trend was ",
|
|
|
round(SST.full.lintrend,3), ' C<>/year with p=', round(SST.full.pvalNCV_ECall,10), sep=" ")) +
|
|
|
theme(plot.title=element_text(face="bold", size=9)) +
|
|
|
geom_smooth(method='lm',fill="green", formula=y~x, colour="darkgreen", size = 0.5) +
|
|
|
stat_smooth(method=gam, formula=y~s(x, k=ggplotGAM.k), se=T, size=0.5) +
|
|
|
ylab("Sea Surface Temperature [C<>]") + xlab("Time")
|
|
|
|
|
|
#export to png
|
|
|
png.name <- paste('./Output/', Case.Study, '/', Estuary, '/Trends_tasmean_full_period_', Sys.Date(),".png", sep="")
|
|
|
png(file = png.name, width = 10.5, height = 7, units='in', res=500)
|
|
|
grid.arrange(p1air,ncol=1)
|
|
|
dev.off()
|
|
|
png.name <- paste('./Output/', Case.Study, '/', Estuary, '/Trends_rivTmean_full_period_', Sys.Date(),".png", sep="")
|
|
|
png(file = png.name, width = 10.5, height = 7, units='in', res=500)
|
|
|
grid.arrange(p1riv,ncol=1)
|
|
|
dev.off()
|
|
|
png.name <- paste('./Output/', Case.Study, '/', Estuary, '/Trends_SSTmean_full_period_', Sys.Date(),".png", sep="")
|
|
|
png(file = png.name, width = 10.5, height = 7, units='in', res=500)
|
|
|
grid.arrange(p1sst,ncol=1)
|
|
|
dev.off()
|
|
|
##################################### Full Time Period
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
######################### 1990-present
|
|
|
combined.TS <- window(merge(AirT.full.TS, window(RivT.full.TS, start=as.Date("1995-01-01"), end=end(RivT.full.TS)), SST.full.TS, all=T), start=as.Date("1990-01-01"), end=end(AirT.full.TS))
|
|
|
combined.df <- data.frame(combined.TS)
|
|
|
colnames(combined.df) <- c('tasmean','rivTmean', 'SSTmean')
|
|
|
|
|
|
#Air temp
|
|
|
p1air <- ggplot(combined.df, aes(y=tasmean, x=index(combined.TS))) + geom_line(alpha=0.5) +
|
|
|
ggtitle(paste(Estuary, " - Linear and smooth trend in catchment airT (SILO) lin trend was ",
|
|
|
round(AirT.lintrend,3), ' C<>/year with p=', round(AirT.pvalNCV_ECall,10), sep=" ")) +
|
|
|
theme(plot.title=element_text(face="bold", size=9)) +
|
|
|
geom_smooth(method='lm',fill="green", formula=y~x, colour="darkgreen", size = 0.5) +
|
|
|
stat_smooth(method=gam, formula=y~s(x, k=ggplotGAM.k), se=T, size=0.5) +
|
|
|
ylab("Air Temperature [C<>]") + xlab("Time")
|
|
|
|
|
|
#Riv temp
|
|
|
p1riv <- ggplot(combined.df, aes(y=rivTmean, x=index(combined.TS))) + geom_line(alpha=0.5) +
|
|
|
ggtitle(paste(Estuary, " - Linear and smooth trend in river temperature (GAUGE) lin trend was ",
|
|
|
round(RivT.lintrend,3), ' C<>/year with p=', round(RivT.pvalNCV_ECall,10), sep=" ")) +
|
|
|
theme(plot.title=element_text(face="bold", size=9)) +
|
|
|
geom_smooth(method='lm',fill="green", formula=y~x, colour="darkgreen", size = 0.5) +
|
|
|
stat_smooth(method=gam, formula=y~s(x, k=ggplotGAM.k), se=T, size=0.5) +
|
|
|
ylab("River Temperature [C<>]") + xlab("Time")
|
|
|
|
|
|
#Sea temp
|
|
|
p1sst <- ggplot(combined.df, aes(y=SSTmean, x=index(combined.TS))) + geom_line(alpha=0.5) +
|
|
|
ggtitle(paste(Estuary, " - Linear and smooth trend in sea surface temperature (NNRP) lin trend was ",
|
|
|
round(SST.lintrend,3), ' C<>/year with p=', round(SST.pvalNCV_ECall,10), sep=" ")) +
|
|
|
theme(plot.title=element_text(face="bold", size=9)) +
|
|
|
geom_smooth(method='lm',fill="green", formula=y~x, colour="darkgreen", size = 0.5) +
|
|
|
stat_smooth(method=gam, formula=y~s(x, k=ggplotGAM.k), se=T, size=0.5) +
|
|
|
ylab("Sea Surface Temperature [C<>]") + xlab("Time")
|
|
|
|
|
|
#export to png
|
|
|
png.name <- paste('./Output/', Case.Study, '/', Estuary, '/Trends_tasmean_1990-present_', Sys.Date(),".png", sep="")
|
|
|
png(file = png.name, width = 10.5, height = 7, units='in', res=500)
|
|
|
grid.arrange(p1air,ncol=1)
|
|
|
dev.off()
|
|
|
png.name <- paste('./Output/', Case.Study, '/', Estuary, '/Trends_rivTmean_1990-present_', Sys.Date(),".png", sep="")
|
|
|
png(file = png.name, width = 10.5, height = 7, units='in', res=500)
|
|
|
grid.arrange(p1riv,ncol=1)
|
|
|
dev.off()
|
|
|
png.name <- paste('./Output/', Case.Study, '/', Estuary, '/Trends_SSTmean_1990-present_', Sys.Date(),".png", sep="")
|
|
|
png(file = png.name, width = 10.5, height = 7, units='in', res=500)
|
|
|
grid.arrange(p1sst,ncol=1)
|
|
|
dev.off()
|
|
|
|
|
|
#export an overview plot as png
|
|
|
gA <- ggplotGrob(p1air)
|
|
|
gB <- ggplotGrob(p1riv)
|
|
|
gC <- ggplotGrob(p1sst)
|
|
|
gA$widths <- gB$widths
|
|
|
gC$widths <- gB$widths
|
|
|
png.name <- paste('./Output/', Case.Study, '/', Estuary, '/Trends_Summary_1990-present_', Sys.Date(),".png", sep="")
|
|
|
png(file = png.name, width = 10.5, height = 7, units='in', res=500)
|
|
|
grid.arrange(gA,gB ,gC,ncol=1)
|
|
|
dev.off()
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|