Compare commits
No commits in common. '9eda19058ec68183bb67cd0f92826be077638274' and 'c5f2f4f743a590b155e2aaebd3de661556b454c7' have entirely different histories.
9eda19058e
...
c5f2f4f743
@ -1,227 +0,0 @@
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
Created on Thu Jun 28 12:15:17 2018
|
||||
|
||||
@author: z5025317
|
||||
"""
|
||||
|
||||
# -*- coding: utf-8 -*-
|
||||
#==========================================================#
|
||||
#Last Updated - March 2018
|
||||
#@author: z5025317 Valentin Heimhuber
|
||||
#code for creating future climate variability deviation plots for NARCLIM variables.
|
||||
#This code is derived from P1_NARCliM_First_Pass_variab_deviation_plots.py to deal with variables for which we don't have NARCLIM deltas.
|
||||
#These variables are Sea Level, Acidity and...
|
||||
# it uses CSIRO CC in Australia CMIP-5 Deltas instead
|
||||
|
||||
|
||||
#==========================================================#
|
||||
#Load packages
|
||||
#==========================================================#
|
||||
import numpy as np
|
||||
import os
|
||||
import pandas as pd
|
||||
import glob
|
||||
import matplotlib
|
||||
import matplotlib.pyplot as plt
|
||||
from datetime import datetime
|
||||
from datetime import timedelta
|
||||
from matplotlib.backends.backend_pdf import PdfPages
|
||||
from ggplot import *
|
||||
import csv
|
||||
matplotlib.style.use('ggplot')
|
||||
# Load my own functions
|
||||
os.chdir('C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdoc/Projects/Paper#1/Analysis/Code')
|
||||
import climdata_fcts as fct
|
||||
import silo as sil
|
||||
import re
|
||||
#==========================================================#
|
||||
|
||||
#==========================================================#
|
||||
# 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/')
|
||||
#==========================================================#
|
||||
|
||||
|
||||
#==========================================================#
|
||||
#set input parameters
|
||||
#==========================================================#
|
||||
Case_Study_Name = 'CASESTUDY2'
|
||||
Casestudy2_csv_path = "Data/NARCLIM_Site_CSVs/CASESTUDY2/CASESTDUY2_NARCLIM_Point_Sites.csv"
|
||||
|
||||
Estuary = 'HUNTER' # 'Belongil'
|
||||
Base_period_start = '1970-01-01' #Start of interval for base period of climate variability
|
||||
Base_period_DELTA_start = '1990-01-01' #Start of interval used for calculating mean and SD for base period (Delta base period)
|
||||
Base_period_end = '2009-01-01' #use last day that's not included in period as < is used for subsetting
|
||||
Clim_var_type = 'Acidity' #Name of climate variable in NARCLIM models '*' will create pdf for all variables in folder
|
||||
[minplotDelta, maxplotDelta]=[0,1]
|
||||
|
||||
#Source for present day climate data (historic time series) can be either: 'Station' or 'SILO'
|
||||
Location = 'Estuary' # pick locaiton for extracting the SILO data can be: Estuary, Catchment, or Ocean
|
||||
|
||||
presentdaybar = False #include a bar for present day variability in the plots?
|
||||
present_day_plot = 'yes' #save a time series of present day
|
||||
Version = "V1"
|
||||
Stats = 'mindaily' #'maxdaily' #maximum takes the annual max Precipitation instead of the sum
|
||||
ALPHA_figs = 1 #Set alpha of figure background (0 makes everything around the plot panel transparent)
|
||||
#==========================================================#
|
||||
|
||||
|
||||
|
||||
#==========================================================#
|
||||
#set directory path for output files
|
||||
output_directory = 'Output/' + Case_Study_Name + '/' + Estuary + '/' + '/Clim_Deviation_Plots/'
|
||||
#output_directory = 'J:/Project wrl2016032/NARCLIM_Raw_Data/Extracted'
|
||||
if not os.path.exists(output_directory):
|
||||
os.makedirs(output_directory)
|
||||
print('-------------------------------------------')
|
||||
print("output directory folder didn't exist and was generated")
|
||||
print('-------------------------------------------')
|
||||
print('-------------------')
|
||||
#==========================================================#
|
||||
|
||||
|
||||
#==========================================================#
|
||||
#read ensemble change CSV file
|
||||
#==========================================================#
|
||||
|
||||
clim_var_csv_path = './Data/Non_NARCLIM/AUZ_CC_CSIRO_Deltas/Brisbane_Ocean_Variables.csv'
|
||||
df = pd.read_csv(clim_var_csv_path, index_col=0, parse_dates = True)
|
||||
df = df.filter(regex= 'RCP 8.5|Variable')
|
||||
df = df[df.Variable=='Ocean pH']
|
||||
Ensemble_Delta_df = pd.DataFrame(index=['Median', '10th', '90th'], columns =['near', 'far'])
|
||||
Ensemble_Delta_df['near'] = re.split('\(|\)|\ to |', df.iloc[0,1])[0:3]
|
||||
Ensemble_Delta_df['far'] = re.split('\(|\)|\ to |', df.iloc[0,2])[0:3]
|
||||
|
||||
Ensemble_Delta_df = Ensemble_Delta_df.astype(float).fillna(0.0)
|
||||
#==========================================================#
|
||||
|
||||
|
||||
#==========================================================#
|
||||
#load present day climate variable time series
|
||||
#==========================================================#
|
||||
#Present_day_df,minplotDelta, maxplotDelta = fct.import_present_day_climdata_csv(Estuary, Clim_var_type)
|
||||
Present_day_Var_CSV = glob.glob('./Data/Non_NARCLIM/Hunter_Lower_Acidity2.csv')
|
||||
Present_day_df = pd.read_csv(Present_day_Var_CSV[0],parse_dates=True, index_col=0)
|
||||
Present_day_df.columns = ['PH']
|
||||
df = Present_day_df
|
||||
df.index = df.index.to_period('D')
|
||||
Present_day_df = df
|
||||
Present_day_df = df.groupby(df.index).mean().reset_index()
|
||||
Present_day_df.columns = ['Date','PH']
|
||||
Present_day_df.index = Present_day_df['Date']
|
||||
Present_day_df= Present_day_df['PH']
|
||||
#Present_day_df.index = pd.to_datetime(Present_day_df.index,format='%Y%m%d')
|
||||
|
||||
if Stats == 'mindaily':
|
||||
Present_day_df_annual = Present_day_df.resample('A').min()
|
||||
Present_day_df_annual = Present_day_df_annual.replace(0, np.nan)
|
||||
if Stats == 'mean':
|
||||
Present_day_df_annual = Present_day_df.resample('A').mean()
|
||||
Present_day_df_annual = Present_day_df_annual.replace(0, np.nan)
|
||||
|
||||
Present_Day_ref_df = Present_day_df_annual
|
||||
|
||||
#Subset to present day variability period
|
||||
Present_Day_ref_df = pd.DataFrame(Present_Day_ref_df.loc[(Present_Day_ref_df.index >= Base_period_start) & (Present_Day_ref_df.index <= Base_period_end)])
|
||||
#Subset to present day variability delta base period (this is the statistical baseline used for present day conditions)
|
||||
Present_Day_Delta_ref_df = pd.DataFrame(Present_Day_ref_df.loc[(Present_Day_ref_df.index >= Base_period_DELTA_start) & (Present_Day_ref_df.index <= Base_period_end)])
|
||||
Present_Day_Mean = np.percentile(Present_Day_Delta_ref_df, 50)
|
||||
Present_Day_SD = np.std(Present_Day_Delta_ref_df)
|
||||
#create data frame for floating stacked barplots
|
||||
index=['-2std', '-1std', 'Med', '1std', '2std']
|
||||
columns = ['present','near future', 'far future']
|
||||
Plot_in_df = pd.DataFrame(index=index, columns =columns)
|
||||
#
|
||||
Plot_in_df['present'] = [float(Present_Day_Mean-2*Present_Day_SD),float(Present_Day_Mean-Present_Day_SD), float(Present_Day_Mean),
|
||||
float(Present_Day_Mean+Present_Day_SD), float(Present_Day_Mean+2*Present_Day_SD)]
|
||||
Plot_in_df['near future'] = [float(Present_Day_Mean + Ensemble_Delta_df.near['10th']),np.NaN, float(Present_Day_Mean + Ensemble_Delta_df.near['Median']),
|
||||
np.NaN, float(Present_Day_Mean + Ensemble_Delta_df.near['90th'])]
|
||||
Plot_in_df['far future'] = [float(Present_Day_Mean + Ensemble_Delta_df.far['10th']),np.NaN, float(Present_Day_Mean + Ensemble_Delta_df.far['Median']),
|
||||
np.NaN, float(Present_Day_Mean + Ensemble_Delta_df.far['90th'])]
|
||||
#Create a second data frame that has the values in a way that they can be stacked up in bars with the correct absolute values
|
||||
Plot_in_df2 = pd.DataFrame(index=index, columns =columns )
|
||||
#
|
||||
index=['-2std', '-1std', 'Med', '1std', '2std']
|
||||
columns = ['near future', 'far future']
|
||||
Plot_in_df2 = pd.DataFrame(index=index, columns =columns )
|
||||
|
||||
Plot_in_df2['near future'] = [float(Present_Day_Mean + Ensemble_Delta_df.near['10th']),np.NaN, float(Ensemble_Delta_df.near['Median']-Ensemble_Delta_df.near['10th']),
|
||||
np.NaN, float(Ensemble_Delta_df.near['90th']-Ensemble_Delta_df.near['Median'])]
|
||||
Plot_in_df2['far future'] = [float(Present_Day_Mean + Ensemble_Delta_df.far['10th']),np.NaN, float(Ensemble_Delta_df.far['Median']-Ensemble_Delta_df.far['10th']),
|
||||
np.NaN, float(Ensemble_Delta_df.far['90th']-Ensemble_Delta_df.far['Median'])]
|
||||
#transpose the data frame
|
||||
Plot_in_df_tp = Plot_in_df2.transpose()
|
||||
#do the individual plots
|
||||
xmin = int(min(Plot_in_df.min(axis=1))-minplotDelta)
|
||||
xmax = int(max(Plot_in_df.max(axis=1))+maxplotDelta)
|
||||
#define colour scheme
|
||||
#likert_colors = ['none', 'firebrick','firebrick','lightcoral','lightcoral']
|
||||
likert_colors = ['none', 'darkblue', 'darkblue','cornflowerblue','cornflowerblue']
|
||||
uni_colors = ['none', 'cornflowerblue', 'cornflowerblue','cornflowerblue','cornflowerblue']
|
||||
|
||||
#plot the stacked barplot
|
||||
fig = plt.figure(figsize=(14,8))
|
||||
ax=plt.subplot(2,4,3)
|
||||
Plot_in_df_tp.plot.bar(stacked=True, color=uni_colors, edgecolor='none', legend=False, ax=ax, width=0.5)
|
||||
df = pd.DataFrame(Plot_in_df.iloc[2,[1,2]])
|
||||
index2=['Med']
|
||||
columns2 = ['near future', 'far future']
|
||||
Plot_in_df3 = pd.DataFrame(index=index2, columns =columns2)
|
||||
Plot_in_df3['near future'] = df.iloc[1,0]
|
||||
Plot_in_df3['far future'] = df.iloc[0,0]
|
||||
Plot_in_df3 = Plot_in_df3.transpose()
|
||||
plt.plot(Plot_in_df3['Med'], linestyle="", markersize=52,
|
||||
marker="_", color='darkblue', label="Median")
|
||||
z = plt.axhline(float(Present_Day_Mean-2*Present_Day_SD), linestyle='-', color='black', alpha=.5)
|
||||
z.set_zorder(-1)
|
||||
z = plt.axhline(float(Present_Day_Mean+2*Present_Day_SD), linestyle='-', color='black', alpha=.5)
|
||||
z.set_zorder(-1)
|
||||
z = plt.axhline(float(Present_Day_Mean-Present_Day_SD), linestyle='--', color='black', alpha=.5)
|
||||
z.set_zorder(-1)
|
||||
z = plt.axhline(float(Present_Day_Mean+Present_Day_SD), linestyle='--', color='black', alpha=.5)
|
||||
z.set_zorder(-1)
|
||||
z = plt.axhline(float(Present_Day_Mean), linestyle='--', color='red', alpha=.5)
|
||||
z.set_zorder(-1)
|
||||
plt.ylim(xmin, xmax)
|
||||
plt.title(Clim_var_type)
|
||||
ax.grid(False)
|
||||
for tick in ax.get_xticklabels():
|
||||
tick.set_rotation(0)
|
||||
fig.tight_layout()
|
||||
fig.patch.set_alpha(ALPHA_figs)
|
||||
|
||||
#plot the present day time series
|
||||
ax=plt.subplot(2,2,1)
|
||||
Present_Day_ref_df.plot(legend=False, ax=ax)
|
||||
z = plt.axhline(float(Present_Day_Mean-2*Present_Day_SD), linestyle='-', color='black', alpha=.5)
|
||||
z.set_zorder(-1)
|
||||
z = plt.axhline(float(Present_Day_Mean+2*Present_Day_SD), linestyle='-', color='black', alpha=.5)
|
||||
z.set_zorder(-1)
|
||||
z = plt.axhline(float(Present_Day_Mean-Present_Day_SD), linestyle='--', color='black', alpha=.5)
|
||||
z.set_zorder(-1)
|
||||
z = plt.axhline(float(Present_Day_Mean+Present_Day_SD), linestyle='--', color='black', alpha=.5)
|
||||
z.set_zorder(-1)
|
||||
z = plt.axhline(float(Present_Day_Mean), linestyle='--', color='red', alpha=.5)
|
||||
z.set_zorder(-1)
|
||||
#fig.patch.set_facecolor('deepskyblue')
|
||||
fig.tight_layout()
|
||||
fig.patch.set_alpha(ALPHA_figs)
|
||||
plt.title(Clim_var_type + ' ' + Stats + ' annual present day')
|
||||
plt.ylim(xmin, xmax)
|
||||
ax.grid(False)
|
||||
|
||||
if presentdaybar == False:
|
||||
out_file_name = Estuary + '_' + Clim_var_type + '_' + Stats + '_' + Base_period_start + '_' + Base_period_end + Version + '_' + '_NPDB.png'
|
||||
out_path = output_directory + '/' + out_file_name
|
||||
fig.savefig(out_path)
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@ -1,165 +0,0 @@
|
||||
#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()
|
||||
Loading…
Reference in New Issue