From 7302cf843c90ef3384b6eb376b590f3cea04e5b2 Mon Sep 17 00:00:00 2001 From: Valentin Heimhuber Date: Tue, 26 Jun 2018 22:16:51 +1000 Subject: [PATCH] #added code for dealing with OEH buoy data on SST and Hsig --- ...ot_time_series_with_trends_Waves_SST_OEH.R | 165 ++++++++++++++++++ ...ggplot_time_series_with_trends_sea_level.R | 113 ++++++++++++ 2 files changed, 278 insertions(+) create mode 100644 Analysis/Code/ggplot_time_series_with_trends_Waves_SST_OEH.R create mode 100644 Analysis/Code/ggplot_time_series_with_trends_sea_level.R diff --git a/Analysis/Code/ggplot_time_series_with_trends_Waves_SST_OEH.R b/Analysis/Code/ggplot_time_series_with_trends_Waves_SST_OEH.R new file mode 100644 index 0000000..ccc5a13 --- /dev/null +++ b/Analysis/Code/ggplot_time_series_with_trends_Waves_SST_OEH.R @@ -0,0 +1,165 @@ +#R code for creating ggplots of time series with smooth (GAM) and linear term for the OEH buoy data on wave heigth and SST + + +###################### +#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 <- "RICHMOND" +Climvar <- 'tasmean' +ggplotGAM.k <- 7 +###################### + + +###################### +Output.Directory <- paste('./Output/', Case.Study, '/', Estuary,'/Recent_Trends/', sep="") +if (file.exists(Output.Directory)){ + print('output folder already existed and was not created again') +} else { + dir.create(file.path(Output.Directory)) + print('output folder did not exist and was created') +} +###################### + + +###################### +#Set input file paths +###################### +pattern = paste(Estuary, '.*.csv', sep="") +Buoy_CSV_Path <- list.files("./Data/Ocean_Data/Waves_SST/", pattern, full.names=T) +Buoy.name <- substr(list.files("./Data/Ocean_Data/Waves_SST/", pattern, full.names=F), nchar(Estuary)+2, nchar(Estuary)+7) + +SST.df <- data.frame(read.csv(AirT_CSV_Path, colClasses = "character")) +colnames(SST.df) <- as.character(SST.df[9,]) +SST.df = SST.df[-c(1:9), ] + +colnames(SST.df) <- gsub(x=colnames(SST.df), pattern="Date/Time",replacement="Date",fixed=T) +colnames(SST.df) <- gsub(x=colnames(SST.df), pattern=" ",replacement="",fixed=T) +colnames(SST.df) <- gsub(x=colnames(SST.df), pattern="Hsig(m)",replacement="Hsig",fixed=T) +colnames(SST.df) <- gsub(x=colnames(SST.df), pattern="SeaTemp(C)",replacement="SST",fixed=T) +SST.df$Date <- substr(SST.df$Date, 1, 10) + + +###################### +#Hsig +###################### +Hsig.full.TS <- zoo(as.numeric(SST.df$Hsig), order.by= as.Date(SST.df$Date, format = "%d/%m/%Y")) #=daily time series of rainfall for creation of clean, daily TS of ET and Q +Hsig.full.TS <- aggregate(Hsig.full.TS, index(Hsig.full.TS) , FUN=mean) +Hsig.full.TS <- aggregate(Hsig.full.TS , as.yearmon(time(Hsig.full.TS )),FUN = mean) +Hsig.full.TS <- Hsig.full.TS [1:length(Hsig.full.TS)-1] +Hsig.TS <- window(Hsig.full.TS, start=as.Date("1990-01-01"), end=as.Date("2018-01-01")) +Hsig.full.df <- data.frame(Hsig.full.TS) +Hsig.df <- data.frame(Hsig.TS) +str(Hsig.df) +colnames(Hsig.df) <- 'Hsigmean' +colnames(Hsig.full.df) <- 'Hsigmean' +#rid of outliers +#Hsig.full.df$MSL <- replace(Hsig.full.df$MSL, which(Hsig.full.df$MSL > 4), NA) +Hsig.full.df$Julday1 <- seq(1,length(Hsig.full.df[,1]),1) +linear.trend.model_EC_all <- lm(Hsigmean ~ Julday1, Hsig.full.df) +Hsig.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4] +Hsig.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 12 +###################### + +###################### +#SST +###################### +SST.full.TS <- zoo(as.numeric(SST.df$SST), order.by= as.Date(SST.df$Date, format = "%d/%m/%Y")) #=daily time series of rainfall for creation of clean, daily TS of ET and Q +SST.full.TS <- na.omit(SST.full.TS) +SST.full.TS <- aggregate(SST.full.TS, index(SST.full.TS) , FUN=mean) +SST.full.TS <- SST.full.TS[1:length(SST.full.TS)-1] +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) +tail(SST.full.TS) +colnames(SST.df) <- 'SSTmean' +colnames(SST.full.df) <- 'SSTmean' +#rid of outliers +#SST.full.df$MSL <- replace(SST.full.df$MSL, which(SST.full.df$MSL > 4), NA) +SST.full.df$Julday1 <- seq(1,length(SST.full.df[,1]),1) +linear.trend.model_EC_all <- lm(SSTmean ~ Julday1, SST.full.df) +SST.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4] +SST.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 365 +###################### + + +###################### +#Plot +###################### + +##################################### Full Time Period +p1air <- ggplot(Hsig.full.df, aes(y=Hsigmean, x=index(Hsig.full.TS))) + geom_line(alpha=0.5) + + ggtitle(paste(Estuary, " - Linear and smooth trend in significant wave height (OEH Buoy:",Buoy.name, "| lin trend was ", + round(Hsig.lintrend,3), ' m/year with p=', round(Hsig.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=4), se=T, size=0.5, col="red") + + ylab("Significant Wave Height [m]") + xlab("Time") + +#export to png +png.name <- paste(Output.Directory, '/Trends_Significant_Wave_Height_',Buoy.name, '_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() + +#multiple smooths +p1air <- ggplot(Hsig.full.df, aes(y=Hsigmean, x=index(Hsig.full.TS))) + geom_line(alpha=0.5) + + ggtitle(paste(Estuary, " - Linear and smooth trend in monthly mean sea level (OEH Buoy:",Buoy.name, "| lin trend was ", + round(Hsig.lintrend,3), ' m/year with p=', round(Hsig.pvalNCV_ECall,10), sep=" ")) + + theme(plot.title=element_text(face="bold", size=9)) + + stat_smooth(method=gam, formula=y~s(x, k=13), se=T, size=0.5, col="red") + + #stat_smooth(method=gam, formula=y~s(x, k=8), se=T, size=0.5, cor="blue") + + stat_smooth(method=gam, formula=y~s(x, k=5), se=T, size=0.5, col="green") + + ylab("Significant Wave Height [m]") + xlab("Time") + +#export to png +png.name <- paste(Output.Directory, '/Trends_Significant_Wave_Height_',Buoy.name, '_MultiGAM_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() + +##################################### Full Time Period SST +p1air <- 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 SST (OEH Buoy:",Buoy.name, "| lin trend was ", + round(SST.lintrend,3), ' deg 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=4), se=T, size=0.5, col="red") + + ylab("Sea Surface Temperature [C°]") + xlab("Time") + +#export to png +png.name <- paste(Output.Directory, '/Trends_SST_',Buoy.name, '_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() + +#multiple smooths +p1air <- 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 monthly mean sea level (OEH Buoy:",Buoy.name, "| lin trend was ", + round(SST.lintrend,3), ' deg C/year with p=', round(SST.pvalNCV_ECall,10), sep=" ")) + + theme(plot.title=element_text(face="bold", size=9)) + + stat_smooth(method=gam, formula=y~s(x, k=13), se=T, size=0.5, col="red") + + #stat_smooth(method=gam, formula=y~s(x, k=8), se=T, size=0.5, cor="blue") + + stat_smooth(method=gam, formula=y~s(x, k=5), se=T, size=0.5, col="green") + + ylab("Sea Surface Temperature [C°]") + xlab("Time") + +#export to png +png.name <- paste(Output.Directory, '/Trends_SST_',Buoy.name, '_MultiGAM_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() diff --git a/Analysis/Code/ggplot_time_series_with_trends_sea_level.R b/Analysis/Code/ggplot_time_series_with_trends_sea_level.R new file mode 100644 index 0000000..f07567d --- /dev/null +++ b/Analysis/Code/ggplot_time_series_with_trends_sea_level.R @@ -0,0 +1,113 @@ +#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 <- "NADGEE" +Climvar <- 'tasmean' +ggplotGAM.k <- 7 +###################### + + +###################### +Output.Directory <- paste('./Output/', Case.Study, '/', Estuary,'/Recent_Trends/', sep="") +if (file.exists(Output.Directory)){ + print('output folder already existed and was not created again') +} else { + dir.create(file.path(Output.Directory)) + print('output folder did not exist and was created') +} +###################### + + +###################### +#Set input file paths +###################### +AirT_CSV_Path <- paste("./Data/Ocean_Data/BOM_monthly_SL_",Estuary, "_Eden.txt", sep="") + + +dat = readLines(AirT_CSV_Path) +dat = as.data.frame(do.call(rbind, strsplit(dat, split=" {2,10}")), stringsAsFactors=FALSE) +colnames(dat) <-dat[3,] + +if (Estuary=='NADGEE'){ + dat2 = dat[-c(1:11), ] + dat2 = dat2[-(365:381),] + dat2 = dat2[,-1] + } + +if (Estuary=='HUNTER'){ + dat2 = dat[-c(1:3), ] + dat2 = dat2[-(720:726),] + dat2 = dat2[-(1:108),] + dat2 = dat2[,-1] + } + + + +dat2$Date <- as.yearmon(dat2[,1], "%m %Y") + +SeaLev.df <- dat2 +head(SeaLev.df) +SeaLev.df$MSL <- as.numeric(SeaLev.df$Mean) + +#rid of outliers +SeaLev.df$MSL <- replace(SeaLev.df$MSL, which(SeaLev.df$MSL > 4), NA) + +SeaLev.df$Julday1 <- seq(1,length(SeaLev.df[,1]),1) +linear.trend.model_EC_all <- lm(MSL ~ Julday1, SeaLev.df) +SeaLev.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4] +SeaLev.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 12 + + + +###################### +#Plot +###################### + +##################################### Full Time Period +p1air <- ggplot(SeaLev.df, aes(y=MSL, x=Date)) + geom_line(alpha=0.5) + + ggtitle(paste(Estuary, " - Linear and smooth trend in monthly mean sea level (BOM Gauge) | lin trend was ", + round(100*SeaLev.lintrend,3), ' cm/year with p=', round(SeaLev.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=4), se=T, size=0.5, col="red") + + ylab("Monthly Mean Sea Level [m]") + xlab("Time") + +#export to png +png.name <- paste(Output.Directory, '/Trends_MonthlyMeanSeaLevel_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() + +#multiple smooths +p1air <- ggplot(SeaLev.df, aes(y=MSL, x=Date)) + geom_line(alpha=0.5) + + ggtitle(paste(Estuary, " - Linear and smooth trend in monthly mean sea level (BOM Gauge) | lin trend was ", + round(100*SeaLev.lintrend,3), ' cm/year with p=', round(SeaLev.pvalNCV_ECall,10), sep=" ")) + + theme(plot.title=element_text(face="bold", size=9)) + + stat_smooth(method=gam, formula=y~s(x, k=13), se=T, size=0.5, col="red") + + #stat_smooth(method=gam, formula=y~s(x, k=8), se=T, size=0.5, cor="blue") + + stat_smooth(method=gam, formula=y~s(x, k=5), se=T, size=0.5, col="green") + + ylab("Monthly Mean Sea Level [m]") + xlab("Time") + +#export to png +png.name <- paste(Output.Directory, '/Trends_MonthlyMeanSeaLevel_MultiGAM_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()