#added code for dealing with OEH buoy data on SST and Hsig
Valentin Heimhuber
7 years ago
parent
c5f2f4f743
commit
7302cf843c
@ -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()
|
||||
@ -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()
|
||||
Loading…
Reference in New Issue