From 6baeb9c237dbce64531c37c0950862c72a41640c Mon Sep 17 00:00:00 2001 From: tinoheimhuber Date: Thu, 21 Jun 2018 19:11:51 +1000 Subject: [PATCH] #just backup --- Analysis/Code/P1_SILO_Batch_Download.py | 35 +-- Analysis/Code/climdata_fcts.pyc | Bin 0 -> 4926 bytes .../Code/ggplot_time_series_with_trends.R | 220 ++++++++++++++++++ ...ggplot_time_series_with_trends_sea_levle.R | 90 +++++++ 4 files changed, 331 insertions(+), 14 deletions(-) create mode 100644 Analysis/Code/climdata_fcts.pyc create mode 100644 Analysis/Code/ggplot_time_series_with_trends.R create mode 100644 Analysis/Code/ggplot_time_series_with_trends_sea_levle.R diff --git a/Analysis/Code/P1_SILO_Batch_Download.py b/Analysis/Code/P1_SILO_Batch_Download.py index 4bc5450..dfb5379 100644 --- a/Analysis/Code/P1_SILO_Batch_Download.py +++ b/Analysis/Code/P1_SILO_Batch_Download.py @@ -20,14 +20,12 @@ from matplotlib.backends.backend_pdf import PdfPages from ggplot import * matplotlib.style.use('ggplot') import csv -from datetime import datetime - # import own modules # Set working direcotry (where postprocessed NARClIM data is located) os.chdir('C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdoc/Projects/Paper#1/Analysis/Code') -import climdata_fcts as fct +#import climdata_fcts as fct import silo as sil - +#==========================================================# #==========================================================# # Set working direcotry (where postprocessed NARClIM data is located) @@ -39,9 +37,13 @@ os.chdir('C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdo #set input parameters Case_Study_Name = 'CASESTUDY2' Casestudy2_csv_path = "Data/NARCLIM_Site_CSVs/CASESTUDY2/CASESTDUY2_NARCLIM_Point_Sites.csv" -Silo_variables = ['daily_rain', "max_temp", "min_temp", 'et_short_crop'] +Silo_variables = ['daily_rain', "max_temp", "min_temp", 'et_short_crop', 'evap_syn'] +Location = 'Catchment' +startdate= '19600101' +enddate= '20180101' #==========================================================# + #==========================================================# #set directory path for output files output_directory = 'Data/SILO/' + Case_Study_Name + '/' @@ -52,19 +54,25 @@ if not os.path.exists(output_directory): print("output directory folder didn't exist and was generated") print('-------------------------------------------') #==========================================================# + - + +#==========================================================# #read the CSV to extract the lat long and case stuy sites with open(Casestudy2_csv_path, mode='r') as infile: reader = csv.reader(infile) next(reader, None) with open('coors_new.csv', mode='w') as outfile: writer = csv.writer(outfile) - mydict = dict((rows[0],[rows[1],rows[2]]) for rows in reader) - + if Location == 'Estuary': + mydict = dict((rows[0],[rows[1],rows[2]]) for rows in reader) + if Location == 'Ocean': + mydict = dict((rows[0],[rows[3],rows[4]]) for rows in reader) + if Location == 'Catchment': + mydict = dict((rows[0],[rows[5],rows[6]]) for rows in reader) + for Estuary in mydict: - print Estuary - print str(mydict[Estuary][0]) + print Estuary, mydict[Estuary][0], mydict[Estuary][1] #==========================================================# #set directory path for output files @@ -77,11 +85,10 @@ for Estuary in mydict: print('-------------------------------------------') #==========================================================# - output_csv = output_directory_internal + 'SILO_climdata_' + Estuary +'_' + mydict[Estuary][0] + '_' + mydict[Estuary][1] + '.csv' - silo_df = sil.pointdata(Silo_variables, 'Okl9EDxgS2uzjLWtVNIBM5YqwvVcCxOmpd3nCzJh','19900101','20100101', + output_csv = output_directory_internal + 'SILO_climdata_' + Estuary +'_'+ Location +'_' + mydict[Estuary][0] + '_' + mydict[Estuary][1] + '2.csv' + silo_df = sil.pointdata(Silo_variables, 'Okl9EDxgS2uzjLWtVNIBM5YqwvVcCxOmpd3nCzJh',startdate, enddate, None, mydict[Estuary][0], mydict[Estuary][1], True, output_csv) - - +#==========================================================# diff --git a/Analysis/Code/climdata_fcts.pyc b/Analysis/Code/climdata_fcts.pyc new file mode 100644 index 0000000000000000000000000000000000000000..d6785bb0032c5461c8921092c616e9c93fb9b561 GIT binary patch literal 4926 zcmai2TW{Og5gt;qY{_;UJBcrwY{Jd%gIjE~wz5f+Al+^d!$8x-Eu?H0aO<+PcxY3h zD4C=hc_9G>5_Gpf-}{65lloG$58Lk>>SCo%oS4INX3or+nKPF|75=qYIobX6_YE2T zv+@2KPx==CUt}N965p1{lBO;D1)UZoDoU#;`+zgDKO@#~QAze^M+|eaUlzY8Zy}{3 z{*3q~d27l3y!f+? z+h01{_aEJ@KHA!Ls$1I+DwVIjZg$*G9y#CL+p6BXyZ!S@_#^)9ZPKc2jut^k7Pd@Q zCI$oFM48thy@;oT=P91_7$6gATJlE3;>f01Z!8(uN+?Jl(ts=gX^mDCX=WeT65$)m z9u!8a)rn_T7gaJmQQxFu>0N+M)Xq|;>Vz%t$fQnedYO~89pg2QopyF?lFyy=*h|9r z$id9QXn7qS93|dq_#M$ShYXN29-i)P&X86B=tP+kg7B!Dm~58IyNA?5 z(P1k+ZogKU9BgtIB(Sv^_9p* zebWL<`UyZ{ng&WZ?QRx@CP|$zc04EUwmPSdm%wKe`@F+6_tuY6lx9`IM5{6YT~^p%`|4h%(T=&p~BupZmK}?MDtrz zP(5U+hemP0JV_q{;NTCwkbxy|MiW!@g)j1|D2b~6gOEhLP3Qg2I?>`)*CPzi*un5> zw7YSZA1-&gvsJr}mm+F6Fskvy-XGT~>&fFIgOu$gU){z;buZ--)Ezj3ERQ0SssAP- z-`VE;%wZHw+~ZOl6-#-VQYaVHON<1j@l!nMZvk+;aMaGhFX^egq9L?%XW_)55F!ow z4c%wU>Kn1vax|-zSnE9bSq257psvT9(l?EE{AZJ74ofP6fQ4#YU8$4rou=yLFov(S z=+uCgj72?2+O2lthp`u_Pr^8>@kBD7MG&>UOoLiW(r^RKNIzCzvT|P(PAoq|*6K*H z&wv{`?0Ftmi@=N>_HU>FZy}lxu^h0558D~4JUv*DraG%R!}<+s?w34Ynk6EKxr3Ex zg#}qf*1jspMvnH2veBI5GUOMU#bdnAVJm^Hc!rIuW``^rO|;cn$tsXyAwM?;(5NU) z^i>E)OWhct326eK;X$+gFi%M`B)TfKF?=;IQBkd}z?W6oi6&)Oa-mcx4L_otSEOGO ze_kAaLDu^442Eku#>IeR7&DlY{@lp@7`)^D^8d*q4+M{Q$VCx5e=c@G`VayR46p)S z8v4yhz%I;4*wg6JzRG9}3j}x=4@fHQN_i`N6?S4zqZQQ^{gw0`ohRnV^uqt*-l53; zCIc~JgG-FAbtylW(oS%v(I9WMqi!qC+fK*V(z=5NvzvG=Q(K_~rcGi}lesFT_A??i z(ZV25ab9xN@_LGUnRJZ~!l>OC6T`Ts%J4FPt51lw#$gafnMu?E2r78tN$u(LU%hbs zAiIj`t{=GD54X15>ehoTH_dha=sg6N%`!Iw?DUG4anBcC4<##bTW#NsVN;F_<8zS4 zpGNW!{Gh7N^#d=8hNBjW2FryX(CXpH_!rfM8(>LIgu8$X3DDPf!9Uv4MD%O=Qvk7V zTGy8f!hVdS8&<`+byB{KbJW^{K%N75?M zds0A2-8pzz$TJuWY`5QN=py#9NH5Ju+#|H)lEkPQkcE1JWq)VQ#$QuQ64IWffGLx+ zp+;1u*$T?`&uO;#ci527i_>fn$0R>Z$DVZ#~H z9yI}0vADu=m@~(#EGh#qR__**cIRMC21sQT31l#Ay($i=CkNKIxaj9=Q5wiSuW&(n z$4vd5gFndtg>FVzMz2YK5gxM!mvjr%I|%9HGXzzXuT|MtA1tw0aG#|yc|Ip&)%xHv z%L&Psm3;CKsDO%u!$Cz}#t~ex+3npu2~Qi z9~5Yma&tgOW!$zy8x9+}raEOShrmU-iYGC%&X z%rSXd=FMCt{DW?$fim%r_P6@C{wwo5j61BSF|RPEPS8$(i@RRrV3Zm!ZR2v{@B;?&c(w8P#n8cp+Do06buxgG z*H-1zi_FfmU46fJ@U2!(*Br9^0zi?tO5$s!#;d)~ z>m~_JT3hCvDs@cKz&VAHQP-iRZL->GmwckWXRl&|tPjtPb9Mw?6jN=!8VWAD2TAH4 z)9PsJeSwZeBwO@TZNu83@)nJLW$-~W2U=`i>rA=%O^$pM>Biwn*RZc~9aYYE+S_$a zoSN1ll%8OB&tSRVS#CAm56~q8-ESo4%w%_0ZO$ovust=ao>}Mz7q!7P>^@ExS0$z{ zl-}K~J_jj%2_S`Q=zEKJ7Vu@hVBI(uqDNX4tioEmgz=)ag87QRm+fo7UoVvHmlpop z*RA|I#6WRfzujA$}wn{{WE>YrLb>sLtdQK|kB`03kKc$~+|bom5*T x_=Mmm1ivQWFC`PN=;-_}Egtq`9B4N(Ur~C>U@O^6==AVJv0TP8S1y-t{SS9H$1?x` literal 0 HcmV?d00001 diff --git a/Analysis/Code/ggplot_time_series_with_trends.R b/Analysis/Code/ggplot_time_series_with_trends.R new file mode 100644 index 0000000..986a51e --- /dev/null +++ b/Analysis/Code/ggplot_time_series_with_trends.R @@ -0,0 +1,220 @@ +#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() + + + + + + + + + diff --git a/Analysis/Code/ggplot_time_series_with_trends_sea_levle.R b/Analysis/Code/ggplot_time_series_with_trends_sea_levle.R new file mode 100644 index 0000000..58eee02 --- /dev/null +++ b/Analysis/Code/ggplot_time_series_with_trends_sea_levle.R @@ -0,0 +1,90 @@ +#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 <- 7 +###################### + + +###################### +#Set input file paths +###################### +AirT_CSV_Path <- "./Data/Ocean_Data/BOM_monthly_SL_Hunter_Newcastle.txt" + + + + + + +dat = readLines(AirT_CSV_Path) +dat = as.data.frame(do.call(rbind, strsplit(dat, split=" {2,10}")), stringsAsFactors=FALSE) +colnames(dat) <-dat[3,] + +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) +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/', Case.Study, '/', Estuary, '/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/', Case.Study, '/', Estuary, '/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()