valentinh
/
CC_Est_NARC
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Compare commits

base: valentinh:c5f2f4f743a590b155e2aaebd3de661556b454c7
Branches Tags
valentinh:Development1
valentinh:master
..
compare: valentinh:9eda19058ec68183bb67cd0f92826be077638274
Branches Tags
valentinh:Development1
valentinh:master

2 Commits
c5f2f4f743 ... 9eda19058e

Author SHA1 Message Date
Valentin Heimhuber 9eda19058e #several major imporvements on all fronts 7 years ago
Valentin Heimhuber 7302cf843c #added code for dealing with OEH buoy data on SST and Hsig 7 years ago
9 changed files with 1142 additions and 517 deletions
Show Stats

6
Analysis/Code/NARCLIM_Download/R_Preparing_BASH_script_for_NARCLIM_batch_download.R
Unescape Escape View File

@ -17,12 +17,12 @@
# 'wssmean' Surface wind speed standard_name: air_velocity units: m s-1
# 'sstmean' Sea surface temperatuer daily mean
Clim_Var <- 'sstmean'
Clim_Var <- 'rldsmean'
Datatype <- 'T_GCMS' #T_GCMS for GCM forcing, T_NNRP for reanalysis (only 1950-2009)
Biasboolean <- 'False' #use bias corrected data?
Directory <- 'C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdoc/Projects/Paper#1/Data/NARCLIM_Site_CSVs/'
Filename <- 'NARCLIM_Point_Sites.csv'
Directory <- 'C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdoc/Projects/Paper#1/Data/NARCLIM_Site_CSVs/CASESTUDY2/'
Filename <- 'CASESTDUY2_NARCLIM_Point_Sites.csv'
#Load CSV with location names and lat lon coordinates
Location.df <- data.frame(read.csv(paste(Directory, Filename, sep=""), header=T, fileEncoding="UTF-8-BOM"))

551
Analysis/Code/P1_NARCliM_First_Pass_variab_deviation_plots.py
Unescape Escape View File

@ -1,13 +1,14 @@
# -*- coding: utf-8 -*-
#####################################----------------------------------
#==========================================================#
#Last Updated - March 2018
#@author: z5025317 Valentin Heimhuber
#code for creating future climate variability deviation plots for NARCLIM variables.
#Inputs: Uses CSV files that containthe deltas of all 12 NARCLIM models for 1 grid cell at the site of interest, generated with P1_NARCLIM_plots_Windows.py
#This code is used only for the NARCLIM variables - a separate code is used for ocean variables etc that are not part of the NARCLIM ensemble
#####################################----------------------------------
#==========================================================#
#Load packages
#####################################----------------------------------
#==========================================================#
import numpy as np
import os
import pandas as pd
@ -18,234 +19,342 @@ from datetime import datetime
from datetime import timedelta
from matplotlib.backends.backend_pdf import PdfPages
from ggplot import *
import csv
matplotlib.style.use('ggplot')
#plt.rcParams.update(plt.rcParamsDefault)
#
# 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
#==========================================================#
#==========================================================#
# 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
Base_period_start = '1986-01-01'
Base_period_end = '2005-01-01' #use last day that's not included in period as < is used for subsetting
Estuary = 'Nadgee' # 'Belongil'
Clim_var_type = "*" # '*' will create pdf for all variables in folder
Clim_var_type = "pracc*" # '*' will create pdf for all variables in folder
Present_Day_Clim_Var = 'Rainfall' #MaxT, MinT, Rainfall, ET Wind
present_day_plot = 'yes'
Version = "V2"
Stats = 'dailymax' #maximum takes the annual max Precipitation instead of the sum
#####################################----------------------------------
#set directory path for output files
output_directory = 'Output/Clim_Deviation_Plots/'+ Estuary
#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('-------------------')
Clim_Var_CSVs = glob.glob('./Output/' + Estuary + '/' + Estuary + '_' + Clim_var_type[:-1] + '_' + Stats + '*')
#read CSV file
clim_var_csv_path = Clim_Var_CSVs[0]
Filename = os.path.basename(os.path.normpath(clim_var_csv_path))
Clim_var_type = Filename.split('_', 2)[1]
print(Clim_var_type)
Ensemble_Delta_full_df = pd.read_csv(clim_var_csv_path, index_col=0, parse_dates = True)
#Ensemble_Delta_full_df = pd.to_numeric(Ensemble_Delta_full_df)
#
#load present day climate data for the same variable
if Clim_var_type == 'evspsblmean': #ET time series that we have is not in the same format as the other variables, hence the different treatment
Present_day_Var_CSV = glob.glob('./Data/Wheather_Station_Data/**/' + Estuary + '_' + Present_Day_Clim_Var + '*csv')
Present_day_df = pd.read_csv(Present_day_Var_CSV[0])
Dates = pd.to_datetime(Present_day_df.Date)
Present_day_df.index = Dates
Present_day_df = Present_day_df.iloc[:,1]
Present_day_df = Present_day_df.replace(r'\s+', np.nan, regex=True)
Present_day_df = pd.to_numeric(Present_day_df)
Present_day_df.index = Dates
[minplotDelta, maxplotDelta]=[50,50]
#for tasmean, observed min and max T need to be converted into mean T
elif Clim_var_type == 'tasmean':
Present_day_Var_CSV = glob.glob('./Data/Wheather_Station_Data/**/' + Estuary + '_MaxT*csv')
Present_day_df = pd.read_csv(Present_day_Var_CSV[0])
Dates = pd.to_datetime(Present_day_df.Year*10000+Present_day_df.Month*100+Present_day_df.Day,format='%Y%m%d')
Present_day_df.index = Dates
Present_day_MaxT_df = Present_day_df.iloc[:,5]
Present_day_Var_CSV = glob.glob('./Data/Wheather_Station_Data/**/' + Estuary + '_MinT*csv')
Present_day_df = pd.read_csv(Present_day_Var_CSV[0])
Dates = pd.to_datetime(Present_day_df.Year*10000+Present_day_df.Month*100+Present_day_df.Day,format='%Y%m%d')
Present_day_df.index = Dates
Present_day_MinT_df = Present_day_df.iloc[:,5]
Present_day_df = (Present_day_MaxT_df + Present_day_MinT_df)/2
[minplotDelta, maxplotDelta]=[1,2]
elif Clim_var_type == 'tasmax':
[minplotDelta, maxplotDelta]=[1,2]
elif Clim_var_type == 'wssmean' or Clim_var_type == 'wss1Hmaxtstep':
Present_day_Var_CSV = glob.glob('./Data/Wheather_Station_Data/**/' + Estuary + '_' + Present_Day_Clim_Var + '*csv')
Present_day_df = pd.read_csv(Present_day_Var_CSV[0])
Present_day_df.index = Present_day_df[['Year', 'Month', 'Day', 'Hour']].apply(lambda s : datetime(*s),axis = 1)
Present_day_df = Present_day_df.filter(regex= 'm/s')
Present_day_df = Present_day_df.replace(r'\s+', np.nan, regex=True)
Present_day_df['Wind speed measured in m/s'] = Present_day_df['Wind speed measured in m/s'].convert_objects(convert_numeric=True)
[minplotDelta, maxplotDelta]=[1, 1.5]
else:
Present_day_Var_CSV = glob.glob('./Data/Wheather_Station_Data/**/' + Estuary + '_' + Present_Day_Clim_Var + '*csv')
Present_day_df = pd.read_csv(Present_day_Var_CSV[0])
Dates = pd.to_datetime(Present_day_df.Year*10000+Present_day_df.Month*100+Present_day_df.Day,format='%Y%m%d')
Present_day_df.index = Dates
Present_day_df = Present_day_df.iloc[:,5]
[minplotDelta, maxplotDelta]=[50,50]
#create seasonal sums etc.
if (Clim_var_type == 'pracc' or Clim_var_type == 'evspsblmean' or Clim_var_type == 'potevpmean'
or Clim_var_type == 'pr1Hmaxtstep' or Clim_var_type == 'wss1Hmaxtstep'):
if Stats == 'dailymax':
Present_day_df_annual = Present_day_df.resample('A').max()
Present_day_df_annual = Present_day_df_annual.replace(0, np.nan)
else:
Present_day_df_annual = Present_day_df.resample('A').sum()
Present_day_df_annual = Present_day_df_annual.replace(0, np.nan)
Fdf_Seas_means = Present_day_df.resample('Q-NOV').sum() #seasonal means
Fdf_Seas_means = Fdf_Seas_means.replace(0, np.nan)
else:
Present_day_df_annual = Present_day_df.resample('A').mean()
Present_day_df_monthly = Present_day_df.resample('M').mean()
Present_day_df_weekly = Present_day_df.resample('W').mean()
Fdf_Seas_means = Present_day_df.resample('Q-NOV').mean() #seasonal means
#Loop through annual and seasons and create a deviation plot for each.
times = ['annual', 'DJF', 'MAM', 'JJA','SON']
fig = plt.figure(figsize=(14,8))
delta_all_df = pd.DataFrame()
i=1
for temp in times:
temp = 'annual'
#subset the ensemble dataframe for the period used:
if temp == 'annual':
Ensemble_Delta_df = Ensemble_Delta_full_df.iloc[:,range(0,2)]
#
Present_Day_ref_df = Present_day_df_annual
else:
Ensemble_Delta_df = Ensemble_Delta_full_df.filter(regex= temp)
Ensemble_Delta_df.columns = ['near', 'far']
if temp == 'DJF':
Mean_df = Fdf_Seas_means[Fdf_Seas_means.index.quarter==1]
Column_names = ['DJF_near', 'DJF_far']
if temp == 'MAM':
Mean_df = Fdf_Seas_means[Fdf_Seas_means.index.quarter==2]
Column_names = ['MAM_near', 'MAM_far']
if temp == 'JJA':
Mean_df = Fdf_Seas_means[Fdf_Seas_means.index.quarter==3]
Column_names = ['JJA_near', 'JJA_far']
if temp == 'SON':
Mean_df = Fdf_Seas_means[Fdf_Seas_means.index.quarter==4]
Present_Day_ref_df = Mean_df
print(Ensemble_Delta_df.columns.values)
#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)])
Present_Day_Mean = np.percentile(Present_Day_ref_df, 50)
Present_Day_SD = np.std(Present_Day_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[-3:][0]),np.NaN, float(Present_Day_Mean + Ensemble_Delta_df.near[-3:][1]),
np.NaN, float(Present_Day_Mean + Ensemble_Delta_df.near[-3:][2])]
Plot_in_df['far future'] = [float(Present_Day_Mean + Ensemble_Delta_df.far[-3:][0]),np.NaN, float(Present_Day_Mean + Ensemble_Delta_df.far[-3:][1]),
np.NaN, float(Present_Day_Mean + Ensemble_Delta_df.far[-3:][2])]
#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 )
#==========================================================#
#loop through case study estuaries
Estuaries = ['HUNTER', 'RICHMOND', 'NADGEE', 'SHOALHAVEN', 'GEORGES','CATHIE']
Estuaries = ['HUNTER']
for Est in Estuaries:
Plot_in_df2['present'] = [float(Present_Day_Mean-2*Present_Day_SD),float(Present_Day_SD), float(Present_Day_SD),
float(Present_Day_SD), float(Present_Day_SD)]
Plot_in_df2['near future'] = [float(Present_Day_Mean + Ensemble_Delta_df.near[-3:][0]),np.NaN, float(Ensemble_Delta_df.near[-3:][1]-Ensemble_Delta_df.near[-3:][0]),
np.NaN, float(Ensemble_Delta_df.near[-3:][2]-Ensemble_Delta_df.near[-3:][1])]
Plot_in_df2['far future'] = [float(Present_Day_Mean + Ensemble_Delta_df.far[-3:][0]),np.NaN, float(Ensemble_Delta_df.far[-3:][1]-Ensemble_Delta_df.far[-3:][0]),
np.NaN, float(Ensemble_Delta_df.far[-3:][2]-Ensemble_Delta_df.far[-3:][1])]
#transpose the data frame
Plot_in_df_tp = Plot_in_df2.transpose()
#do the individual plots
if temp == 'annual':
xmin = int(min(Plot_in_df.min(axis=1))-minplotDelta)
xmax = int(max(Plot_in_df.max(axis=1))+maxplotDelta)
else:
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']
#plot the stacked barplot
ax=plt.subplot(2,3,i)
Plot_in_df_tp.plot.bar(stacked=True, color=likert_colors, edgecolor='none', 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)
plt.ylim(xmin, xmax)
plt.title(Clim_var_type + ' ' + temp )
ax.grid(False)
for tick in ax.get_xticklabels():
tick.set_rotation(0)
fig.tight_layout()
fig.patch.set_alpha(0)
#reset i to i+1 for next step
if temp == 'MAM':
i=i+2
else:
i=i+1
print(i)
plt.show()
out_file_name = Estuary + '_' + Clim_var_type + '_' + Stats + '_' + '_CC_prio_plot' + Version + '.png'
out_path = output_directory + '/' + out_file_name
fig.savefig(out_path)
#==========================================================#
#set input parameters
#==========================================================#
Case_Study_Name = 'CASESTUDY2'
Casestudy2_csv_path = "Data/NARCLIM_Site_CSVs/CASESTUDY2/CASESTDUY2_NARCLIM_Point_Sites.csv"
#Estuary = 'HUNTER' # 'Belongil'
Estuary = Est # 'Belongil'
print Estuary
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 = 'tasmax' #Name of climate variable in NARCLIM models '*' will create pdf for all variables in folder
PD_Datasource = 'SILO' #Source for present day climate data (historic time series) can be either: 'Station' or 'SILO'
SILO_Clim_var = ['max_temp'] #select the SILO clim variable to be used for base period. - see silo.py for detailed descriptions
Location = 'Estuary' # pick locaiton for extracting the SILO data can be: Estuary, Catchment, or Ocean
if present_day_plot == 'yes':
#print present day climate data
fig = plt.figure(figsize=(5,4))
ax = fig.add_subplot(1, 1, 1)
if temp == 'annual':
xmin = int(min(Plot_in_df.min(axis=1))-minplotDelta)
xmax = int(max(Plot_in_df.max(axis=1))+maxplotDelta)
else:
xmin = int(min(Plot_in_df.min(axis=1))-minplotDelta)
xmax = int(max(Plot_in_df.max(axis=1))+maxplotDelta)
presentdaybar = False #include a bar for present day variability in the plots?
present_day_plot = 'no' #save a time series of present day
Version = "V1"
Stats = 'days_h_35' #'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)
#==========================================================#
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.patch.set_alpha(0)
plt.ylim(13, xmax)
plt.show()
out_file_name = 'C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdoc/Projects/OEH_Coastal_Node_Deliverables/Technical_Report_1/Figures/tasmean_present_day_manual_backgroundTP.png'
out_path = output_directory + '/' + out_file_name
fig.savefig(out_file_name)
# use transparent=True if you want the whole figure with a transparent background
#==========================================================#
#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 CSV file
#==========================================================#
Clim_Var_CSVs = glob.glob('./Output/' + Case_Study_Name + '/' + Estuary + '/' + Estuary + '_' + Clim_var_type + '_' + Stats + '*')
clim_var_csv_path = Clim_Var_CSVs[0]
Filename = os.path.basename(os.path.normpath(clim_var_csv_path))
Clim_var_type = Filename.split('_', 2)[1]
print(Clim_var_type)
Ensemble_Delta_full_df = pd.read_csv(clim_var_csv_path, index_col=0, parse_dates = True)
if Clim_var_type == 'rsdsmean':
Clim_Var_CSVs = glob.glob('./Output/' + Case_Study_Name + '/' + Estuary + '/' + Estuary + '_rsdsmean_' + Stats + '*')
clim_var_csv_path = Clim_Var_CSVs[0]
Filename = os.path.basename(os.path.normpath(clim_var_csv_path))
Clim_var_type = Filename.split('_', 2)[1]
print(Clim_var_type)
Ensemble_Delta_full_df1 = pd.read_csv(clim_var_csv_path, index_col=0, parse_dates = True)
Clim_Var_CSVs = glob.glob('./Output/' + Case_Study_Name + '/' + Estuary + '/' + Estuary + '_rldsmean_' + Stats + '*')
clim_var_csv_path = Clim_Var_CSVs[0]
Filename = os.path.basename(os.path.normpath(clim_var_csv_path))
Clim_var_type = Filename.split('_', 2)[1]
print(Clim_var_type)
Ensemble_Delta_full_df2 = pd.read_csv(clim_var_csv_path, index_col=0, parse_dates = True)
Ensemble_Delta_full_df = Ensemble_Delta_full_df2 + Ensemble_Delta_full_df1
#Ensemble_Delta_full_df = pd.to_numeric(Ensemble_Delta_full_df)
#==========================================================#
#==========================================================#
#load present day climate variable time series
#==========================================================#
if PD_Datasource == 'Station':
Present_day_df,minplotDelta, maxplotDelta = fct.import_present_day_climdata_csv(Estuary, Clim_var_type)
if PD_Datasource == 'SILO':
#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)
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)
if Clim_var_type == 'tasmean':
silo_df = sil.pointdata(["max_temp", "min_temp"], 'Okl9EDxgS2uzjLWtVNIBM5YqwvVcCxOmpd3nCzJh', Base_period_start.replace("-", ""), Base_period_end.replace("-", ""),
None, mydict[Estuary][0], mydict[Estuary][1], False, None)
#take mean of daily min and max temp
Present_day_df = silo_df.iloc[:,[2,4]].mean(axis=1)
else:
silo_df = sil.pointdata(SILO_Clim_var, 'Okl9EDxgS2uzjLWtVNIBM5YqwvVcCxOmpd3nCzJh', Base_period_start.replace("-", ""), Base_period_end.replace("-", ""),
None, mydict[Estuary][0], mydict[Estuary][1], False, None)
Present_day_df = silo_df.iloc[:,[2]]
#set the x and y limit deltas - for plotting only
if Clim_var_type in ['evspsblmean', 'potevpmean']: #ET time series that we have is not in the same format as the other variables, hence the different treatment
[minplotDelta, maxplotDelta]=[50,50]
#for tasmean, observed min and max T need to be converted into mean T
elif Clim_var_type == 'tasmean':
[minplotDelta, maxplotDelta]=[0.2,1]
elif Clim_var_type == 'tasmax':
[minplotDelta, maxplotDelta]=[1,2]
elif Clim_var_type == 'wssmean' or Clim_var_type == 'wss1Hmaxtstep':
[minplotDelta, maxplotDelta]=[1, 1.5]
elif Clim_var_type == 'pracc':
[minplotDelta, maxplotDelta]=[50,100]
elif Clim_var_type in ['rsdsmean', 'rldsmean']:
[minplotDelta, maxplotDelta]=[100,100]
#[minplotDelta, maxplotDelta]=[1,1]
#==========================================================#
#substract a constant from all values to convert from kelvin to celcius (temp)
if Clim_var_type == 'sstmean':
Present_day_df = Present_day_df.iloc[:,0:(len(Present_day_df)-1)]-273.15
Present_day_df.index = Present_day_df.index.to_period('D')
#==========================================================#
#create seasonal sums etc.
#==========================================================#
if Stats == 'maxdaily':
Present_day_df_annual = Present_day_df.resample('A').max()
Present_day_df_annual = Present_day_df_annual.replace(0, np.nan)
Fdf_Seas_means = Present_day_df.resample('Q-NOV').max() #seasonal means
Fdf_Seas_means = Fdf_Seas_means.replace(0, np.nan)
if(Stats[:4] =='days'):
Threshold = int(Stats[-2:])
agg = ('>'+ str(Threshold) + '_count', lambda x: x.gt(Threshold).sum()),
Present_day_df_annual = Present_day_df.resample('A').agg(agg)
Fdf_Seas_means = Present_day_df.resample('Q-NOV').agg(agg) #seasonal means
else:
if Clim_var_type in ['pracc' ,'evspsblmean' ,'potevpmean' ,
'pr1Hmaxtstep' ,'wss1Hmaxtstep',
'rsdsmean', 'rldsmean']:
Present_day_df_annual = Present_day_df.resample('A').sum()
Present_day_df_annual = Present_day_df_annual.replace(0, np.nan)
Fdf_Seas_means = Present_day_df.resample('Q-NOV').sum() #seasonal means
Fdf_Seas_means = Fdf_Seas_means.replace(0, np.nan)
elif Clim_var_type in ['tasmean' ,'tasmax' ,'sstmean']:
Present_day_df_annual = Present_day_df.resample('A').mean()
Present_day_df_monthly = Present_day_df.resample('M').mean()
Present_day_df_weekly = Present_day_df.resample('W').mean()
Fdf_Seas_means = Present_day_df.resample('Q-NOV').mean() #seasonal means
#==========================================================#
#==========================================================#
#Loop through annual and seasons and create a deviation plot for each.
#==========================================================#
times = ['annual', 'DJF', 'MAM', 'JJA','SON']
fig = plt.figure(figsize=(14,8))
delta_all_df = pd.DataFrame()
i=1
for temp in times:
#temp = 'annual'
#subset the ensemble dataframe for the period used:
if temp == 'annual':
Ensemble_Delta_df = Ensemble_Delta_full_df.iloc[:,range(0,2)]
Present_Day_ref_df = Present_day_df_annual
else:
Ensemble_Delta_df = Ensemble_Delta_full_df.filter(regex= temp)
Ensemble_Delta_df.columns = ['near', 'far']
if temp == 'DJF':
Mean_df = Fdf_Seas_means[Fdf_Seas_means.index.quarter==1]
Column_names = ['DJF_near', 'DJF_far']
if temp == 'MAM':
Mean_df = Fdf_Seas_means[Fdf_Seas_means.index.quarter==2]
Column_names = ['MAM_near', 'MAM_far']
if temp == 'JJA':
Mean_df = Fdf_Seas_means[Fdf_Seas_means.index.quarter==3]
Column_names = ['JJA_near', 'JJA_far']
if temp == 'SON':
Mean_df = Fdf_Seas_means[Fdf_Seas_means.index.quarter==4]
Present_Day_ref_df = Mean_df
#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[-3:][0]),np.NaN, float(Present_Day_Mean + Ensemble_Delta_df.near[-3:][1]),
np.NaN, float(Present_Day_Mean + Ensemble_Delta_df.near[-3:][2])]
Plot_in_df['far future'] = [float(Present_Day_Mean + Ensemble_Delta_df.far[-3:][0]),np.NaN, float(Present_Day_Mean + Ensemble_Delta_df.far[-3:][1]),
np.NaN, float(Present_Day_Mean + Ensemble_Delta_df.far[-3:][2])]
#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 )
#
if presentdaybar == False:
index=['-2std', '-1std', 'Med', '1std', '2std']
columns = ['near future', 'far future']
Plot_in_df2 = pd.DataFrame(index=index, columns =columns )
#Plot_in_df2['present'] = [float(0),float(0), float(0),
# float(0), float(0)]
else:
Plot_in_df2['present'] = [float(Present_Day_Mean-2*Present_Day_SD),float(Present_Day_SD), float(Present_Day_SD),
float(Present_Day_SD), float(Present_Day_SD)]
Plot_in_df2['near future'] = [float(Present_Day_Mean + Ensemble_Delta_df.near[-3:][0]),np.NaN, float(Ensemble_Delta_df.near[-3:][1]-Ensemble_Delta_df.near[-3:][0]),
np.NaN, float(Ensemble_Delta_df.near[-3:][2]-Ensemble_Delta_df.near[-3:][1])]
Plot_in_df2['far future'] = [float(Present_Day_Mean + Ensemble_Delta_df.far[-3:][0]),np.NaN, float(Ensemble_Delta_df.far[-3:][1]-Ensemble_Delta_df.far[-3:][0]),
np.NaN, float(Ensemble_Delta_df.far[-3:][2]-Ensemble_Delta_df.far[-3:][1])]
#transpose the data frame
Plot_in_df_tp = Plot_in_df2.transpose()
#do the individual plots
if temp == 'annual':
xmin = int(min(Plot_in_df.min(axis=1))-minplotDelta)
xmax = int(max(Plot_in_df.max(axis=1))+maxplotDelta)
else:
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
if temp == 'annual':
ax=plt.subplot(2,4,3)
else:
ax=plt.subplot(2,4,i)
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 + ' ' + temp )
ax.grid(False)
for tick in ax.get_xticklabels():
tick.set_rotation(0)
fig.tight_layout()
fig.patch.set_alpha(ALPHA_figs)
if temp == 'annual':
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 + ' ' + temp +' present day')
plt.ylim(xmin, xmax)
ax.grid(False)
#if temp == 'MAM':
i=i+4
else:
i=i+1
#plt.show()
if presentdaybar == False:
out_file_name = Estuary + '_' + Clim_var_type + '_' + Stats + '_' + PD_Datasource + '_' + SILO_Clim_var[0] + Version + '_' + '_NPDB.png'
else:
out_file_name = Estuary + '_' + Clim_var_type + '_' + Stats + '_' + PD_Datasource + '_' + SILO_Clim_var[0] + Version + '_' + '.png'
out_path = output_directory + '/' + out_file_name
fig.savefig(out_path)
if present_day_plot == 'yes':
#print present day climate data
fig = plt.figure(figsize=(5,4))
ax = fig.add_subplot(1, 1, 1)
if temp == 'annual':
xmin = int(min(Plot_in_df.min(axis=1))-minplotDelta)
xmax = int(max(Plot_in_df.max(axis=1))+maxplotDelta)
else:
xmin = int(min(Plot_in_df.min(axis=1))-minplotDelta)
xmax = int(max(Plot_in_df.max(axis=1))+maxplotDelta)
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.patch.set_alpha(0)
plt.ylim(13, xmax)
#export plot to png
out_file_name = Estuary + '_' + Clim_var_type + '_' + Stats + '_' + Base_period_start + '_' + Base_period_end + Version + 'Present_Day_Period.png'
out_path = output_directory + '/' + out_file_name
fig.savefig(out_path)
# use transparent=True if you want the whole figure with a transparent background

227
Analysis/Code/P1_NARCliM_First_Pass_variab_deviation_plots_nonNARCLIM.py
Unescape Escape View File

@ -0,0 +1,227 @@
# -*- 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)

2
Analysis/Code/P1_NARCliM_NC_to_CSV_CCRC_SS_BASH_script_readme.txt
Unescape Escape View File

@ -10,7 +10,7 @@ Variables available from NARCLIM (output):
'wss1Hmaxtstep' Max. 1-hour time-window moving averaged surface wind speed units: m s-1 maximum 1-hour time-window moving averaged values from point values 60.0 second
'wssmax' Surface wind speed standard_name: air_velocity units: m s-1 height: 10 m
'wssmean' Surface wind speed standard_name: air_velocity units: m s-1
'sstmean' Daily mean sea surface temperature
Sites:
Northern NSW:

560
Analysis/Code/P1_NARCliM_plots_Windows.py
Unescape Escape View File

@ -26,8 +26,6 @@ os.chdir('C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdo
import climdata_fcts as fct
import silo as sil
#==========================================================#
# 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/')
@ -36,275 +34,295 @@ os.chdir('C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdo
#==========================================================#
#set input parameters
Base_period_start = '1990-01-01'
Case_Study_Name = 'CASESTUDY2'
Base_period_end = '2080-01-01' #use last day that's not included in period as < is used for subsetting
Estuary = 'HUNTER' # 'Belongil'
Clim_var_type = "sstmean" # '*' will create pdf for all variables in folder "pracc*|tasmax*"
plot_pdf = 'yes'
delta_csv = 'yes'
Stats = 'maxdaily'
Version = 'V4'
#==========================================================#
#==========================================================#
#set directory path for output files
output_directory = 'Output/' + Case_Study_Name + '/' + Estuary
#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('-------------------------------------------')
#==========================================================#
#==========================================================#
Estuary_Folder = glob.glob('./Data/NARCLIM_Site_CSVs/'+ Case_Study_Name + '/' + Estuary + '*' )
Clim_Var_CSVs = glob.glob(Estuary_Folder[0] + '/' + Clim_var_type + '*')
#==========================================================#
#==========================================================#
#read CSV files and start analysis
#==========================================================#
#for clim_var_csv_path in Clim_Var_CSVs:
clim_var_csv_path = Clim_Var_CSVs[0]
Filename = os.path.basename(os.path.normpath(clim_var_csv_path))
Clim_var_type = Filename.split('_', 1)[0]
print(clim_var_csv_path)
Full_df = pd.read_csv(clim_var_csv_path, index_col=0, parse_dates = True)
#pandas datestamp index to period (we don't need the 12 pm info in the index (daily periods are enough))
Full_df.index = Full_df.index.to_period('D')
Full_df = Full_df.drop(columns=['period'])
Ncols_df = len(Full_df)
#check data types of columns
#Full_df.dtypes
#==========================================================#
#substract a constant from all values to convert from kelvin to celcius (temp)
if Clim_var_type == 'tasmean' or Clim_var_type == 'tasmax' or Clim_var_type == 'sstmean':
Full_df = Full_df.iloc[:,0:(Ncols_df-1)]-273.15
if Clim_var_type == 'evspsblmean' or Clim_var_type == 'potevpmean':
Full_df = Full_df.iloc[:,0:(Ncols_df-1)]*60*60*24
Fdf_1900_2080 = Full_df
#==========================================================#
#Subset the data to the minimum base period and above (used to set the lenght of the present day climate period)
#Fdf_1900_2080 = Full_df.loc[(Full_df.index >= Base_period_start) & (Full_df.index < Base_period_end)] # not necessary if not using reanalysis models for base period
#==========================================================#
#Aggregate daily df to annual time series
if (Clim_var_type == 'pracc' or Clim_var_type == 'evspsblmean' or Clim_var_type == 'potevpmean'
or Clim_var_type == 'pr1Hmaxtstep' or Clim_var_type == 'wss1Hmaxtstep'):
if(Stats == 'maxdaily'):
Fdf_1900_2080_annual = Fdf_1900_2080.resample('A').max()
Fdf_1900_2080_annual = Fdf_1900_2080_annual.replace(0, np.nan)
Fdf_Seas_means = Fdf_1900_2080.resample('Q-NOV').max() #seasonal means
Fdf_Seas_means = Fdf_Seas_means.replace(0, np.nan)
else:
Fdf_1900_2080_annual = Fdf_1900_2080.resample('A').sum()
Fdf_1900_2080_annual = Fdf_1900_2080_annual.replace(0, np.nan)
Fdf_Seas_means = Fdf_1900_2080.resample('Q-NOV').sum() #seasonal means
Fdf_Seas_means = Fdf_Seas_means.replace(0, np.nan)
else:
if(Stats == 'maxdaily'):
Fdf_1900_2080_annual = Fdf_1900_2080.resample('A').max()
Fdf_1900_2080_annual = Fdf_1900_2080_annual.replace(0, np.nan)
Fdf_Seas_means = Fdf_1900_2080.resample('Q-NOV').max() #seasonal means
Fdf_Seas_means = Fdf_Seas_means.replace(0, np.nan)
else:
Fdf_1900_2080_annual = Fdf_1900_2080.resample('A').mean()
Fdf_Seas_means = Fdf_1900_2080.resample('Q-NOV').mean() #seasonal means
Fdf_1900_2080_means = Fdf_1900_2080.mean()
#==========================================================#
Estuaries = ['HUNTER', 'RICHMOND', 'NADGEE', 'SHOALHAVEN', 'GEORGES','CATHIE']
Estuaries = ['HUNTER']
#==========================================================#
#Select the 3 most representative models (min med and max difference betwen far future and present)
dfall, dfmin, dfmax, dfmed, Min_dif_mod_name, Med_dif_mod_name, Max_dif_mod_name = fct.select_min_med_max_dif_model(Fdf_1900_2080)
#==========================================================#
#==========================================================#
#create a dataframe that has 1 column for each of the three representative models
dfa = Fdf_1900_2080_annual.iloc[:,[0]]
dfa1 = Fdf_1900_2080_annual.iloc[:,[0,3,6]].loc[(Fdf_1900_2080_annual.index >= '1990') & (Fdf_1900_2080_annual.index <= '2009')]
dfa1.columns = [Min_dif_mod_name[:-5], Med_dif_mod_name[:-5], Max_dif_mod_name[:-5]]
dfa2 = Fdf_1900_2080_annual.iloc[:,[1,4,7]].loc[(Fdf_1900_2080_annual.index >= '2020') & (Fdf_1900_2080_annual.index <= '2039')]
dfa2.columns = [Min_dif_mod_name[:-5], Med_dif_mod_name[:-5], Max_dif_mod_name[:-5]]
dfa3 = Fdf_1900_2080_annual.iloc[:,[2,5,8]].loc[(Fdf_1900_2080_annual.index >= '2060') & (Fdf_1900_2080_annual.index <= '2079')]
dfa3.columns = [Min_dif_mod_name[:-5], Med_dif_mod_name[:-5], Max_dif_mod_name[:-5]]
dfall_annual = dfa1.append(dfa2).append(dfa3)
#==========================================================#
#==========================================================#
#Create Deltas of average change for annual and seasonal basis
#==========================================================#
delta_all_df = fct.calculate_deltas_NF_FF2(Fdf_1900_2080_annual, Fdf_Seas_means)
#==========================================================#
if delta_csv == 'yes':
out_file_name = Estuary + '_' + Clim_var_type + '_' + Stats + '_NARCliM_ensemble_changes.csv'
out_path = output_directory + '/' + out_file_name
delta_all_df.to_csv(out_path)
#==========================================================#
#==========================================================#
#create a dataframe that has a single column for present day, near and far future for the (3 selected models)
Full_current_df = Fdf_1900_2080.iloc[:,range(0,3)]
Full_current_df = Full_current_df.stack()
#nearfuture
Full_nearfuture_df = Fdf_1900_2080.iloc[:,range(3,6)]
Full_nearfuture_df = Full_nearfuture_df.stack()
#farfuture
Full_farfuture_df = Fdf_1900_2080.iloc[:,range(6,len(Fdf_1900_2080.columns))]
Full_farfuture_df = Full_farfuture_df.stack()
Summarized_df = pd.concat([Full_current_df, Full_nearfuture_df], axis=1, ignore_index=True)
Summarized_df = pd.concat([Summarized_df, Full_farfuture_df], axis=1, ignore_index=True)
Summarized_df.columns = ['present', 'near', 'far']
#==========================================================#
#==========================================================#
#output some summary plot into pdf
#==========================================================#
if plot_pdf == 'yes':
plotcolours36 = ['darkolivegreen','turquoise', 'lightgreen', 'darkgreen', 'lightpink','slateblue', 'slategray', 'orange', 'tomato', 'peru', 'navy', 'teal',
'darkolivegreen','turquoise', 'lightgreen', 'darkgreen', 'lightpink','slateblue', 'slategray', 'orange', 'tomato', 'peru', 'navy', 'teal',
'darkolivegreen','turquoise', 'lightgreen', 'darkgreen', 'lightpink','slateblue', 'slategray', 'orange', 'tomato', 'peru', 'navy', 'teal']
plotcolours36b = ['tomato', 'royalblue', 'mediumpurple' , 'tomato', 'royalblue', 'mediumpurple' , 'tomato', 'royalblue', 'mediumpurple' , 'tomato', 'royalblue', 'mediumpurple' ,
'tomato', 'royalblue', 'mediumpurple' , 'tomato', 'royalblue', 'mediumpurple' , 'tomato', 'royalblue', 'mediumpurple' , 'tomato', 'royalblue', 'mediumpurple' ,
'tomato', 'royalblue', 'mediumpurple' , 'tomato', 'royalblue', 'mediumpurple' , 'tomato', 'royalblue', 'mediumpurple' , 'tomato', 'royalblue', 'mediumpurple' ]
plotcolours12 = ['darkolivegreen','turquoise', 'lightgreen', 'darkgreen', 'lightpink','slateblue', 'slategray', 'orange', 'tomato', 'peru', 'navy', 'teal']
plotcolours15 = ['darkolivegreen','turquoise', 'lightgreen', 'darkgreen', 'lightpink','slateblue', 'slategray', 'orange', 'tomato', 'peru', 'navy', 'teal', 'lightgreen','lightpink','slateblue']
#plt.cm.Paired(np.arange(len(Fdf_1900_2080_means)))
#write the key plots to a single pdf document
pdf_out_file_name = Clim_var_type + '_' + Stats + '_start_' + Base_period_start + '_NARCliM_summary_' + Version + '.pdf'
pdf_out_path = output_directory +'/' + pdf_out_file_name
#open pdf and add the plots
with PdfPages(pdf_out_path) as pdf:
#barplot of model means
plt.title(Clim_var_type + ' - model means - full period')
ymin = min(Fdf_1900_2080_means)
ymax = max(Fdf_1900_2080_means) + 0.008 *min(Fdf_1900_2080_means)
Fdf_1900_2080_means.plot(kind='bar', ylim=(ymin,ymax), color=plotcolours36)
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
#
neardeltadf=delta_all_df['near']
ymin = min(neardeltadf) + 0.1 *min(neardeltadf)
ymax = max(neardeltadf) + 0.1 * max(neardeltadf)
neardeltadf=delta_all_df['far']
ymin2 = min(neardeltadf) + 0.1 *min(neardeltadf)
ymax2 = max(neardeltadf) + 0.1 * max(neardeltadf)
ymin = min(ymin, ymin2)
if (Clim_var_type == 'tasmax' or Clim_var_type == 'tasmean'):
ymin = 0
ymax = max(ymax, ymax2)
#
# delta barplot for report 1#################################
ax=plt.subplot(2,1,1)
plt.title(Clim_var_type + ' - model deltas - near-present')
neardeltadf=delta_all_df['near']
neardeltadf.plot(kind='bar', color=plotcolours15, ylim=(ymin,ymax), ax=ax)
plt.xticks([])
#ax.xaxis.set_ticklabels([])
#pdf.savefig(bbox_inches='tight', ylim=(ymin,ymax), pad_inches=0.4)
#plt.close()
#
ax=plt.subplot(2,1,2)
plt.title(Clim_var_type + ' - model deltas - far-present')
neardeltadf=delta_all_df['far']
neardeltadf.plot(kind='bar', color=plotcolours15, ylim=(ymin,ymax), ax=ax)
ax.xaxis.grid(False)
#fig.patch.set_alpha(0)
#plt.show()
pdf.savefig(bbox_inches='tight', ylim=(ymin,ymax), pad_inches=0.4)
plt.close()
# end delta barplot for report 1#################################
#
#full period density comparison
plt.title(Clim_var_type + ' - density comparison - full period - all models')
Summarized_df.plot.kde()
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
#full period density comparison
plt.title(Clim_var_type + ' - density comparison - full period - max delta model')
xmin = float(max(np.nanpercentile(Fdf_1900_2080.filter(regex= Max_dif_mod_name[:-5]),50) - 4 * np.std(Fdf_1900_2080.filter(regex= Max_dif_mod_name[:-5]))))
xmax = float(max(np.nanpercentile(Fdf_1900_2080.filter(regex= Max_dif_mod_name[:-5]),50) + 4 * np.std(Fdf_1900_2080.filter(regex= Max_dif_mod_name[:-5]))))
Fdf_1900_2080.filter(regex= Max_dif_mod_name[:-5]).plot.kde(xlim=(xmin,xmax))
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
#annual box
plt.title(Clim_var_type + ' - Annual means/sums for max diff model')
Fdf_1900_2080_annual.boxplot(rot=90)
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
#
#daily box
plt.title(Clim_var_type + ' - Daily means/sums')
Fdf_1900_2080.boxplot(rot=90)
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
# time series plot annual ALL models
plt.title(Clim_var_type + ' - Time series - all models')
Mod_order = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,19,20,21,16,17,18,22,23,24,31,32,33,25,26,27,28,29,30,34,35,36]
test = Fdf_1900_2080_annual
Mod_Names = test.columns
New_Mod_Name = []
for i in range(0,len(Mod_Names)):
New_Mod_Name.append(str(Mod_order[i]+10) + '_' + Mod_Names[i])
test.columns = New_Mod_Name
test_sorted = test.reindex_axis(sorted(test.columns), axis=1)
colnamest = test.columns
test_sorted.columns = [w[3:-5] for w in colnamest]
test_sorted.plot(legend=False, color = plotcolours36)
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
# time series plot annual ALL models
plt.title(Clim_var_type + ' - Time series - representative models')
dfall_annual.plot(legend=False)
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
# seasonal mean boxplots
ymin = min(Fdf_Seas_means[Fdf_Seas_means.index.quarter==1].mean())
ymax = max(Fdf_Seas_means[Fdf_Seas_means.index.quarter==1].mean())
plt.title(Clim_var_type + ' - DJF Summer means/sums')
pd.DataFrame(Fdf_Seas_means[Fdf_Seas_means.index.quarter==1].mean()).plot(kind='bar', ylim=(ymin,ymax))
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
plt.title(Clim_var_type + ' - DJF Summer means/sums')
Fdf_Seas_means[Fdf_Seas_means.index.quarter==1].boxplot(rot=90)
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
ymin = min(Fdf_Seas_means[Fdf_Seas_means.index.quarter==2].mean())
ymax = max(Fdf_Seas_means[Fdf_Seas_means.index.quarter==2].mean())
plt.title(Clim_var_type + ' - MAM Autumn means/sums')
pd.DataFrame(Fdf_Seas_means[Fdf_Seas_means.index.quarter==2].mean()).plot(kind='bar', ylim=(ymin,ymax))
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
plt.title(Clim_var_type + ' - MAM Autumn means/sums')
Fdf_Seas_means[Fdf_Seas_means.index.quarter==2].boxplot(rot=90)
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
ymin = min(Fdf_Seas_means[Fdf_Seas_means.index.quarter==3].mean())
ymax = max(Fdf_Seas_means[Fdf_Seas_means.index.quarter==3].mean())
plt.title(Clim_var_type + ' - JJA Winter means/sums')
pd.DataFrame(Fdf_Seas_means[Fdf_Seas_means.index.quarter==3].mean()).plot(kind='bar', ylim=(ymin,ymax))
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
plt.title(Clim_var_type + ' - JJA Winter means/sums')
Fdf_Seas_means[Fdf_Seas_means.index.quarter==3].boxplot(rot=90)
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
ymin = min(Fdf_Seas_means[Fdf_Seas_means.index.quarter==4].mean())
ymax = max(Fdf_Seas_means[Fdf_Seas_means.index.quarter==4].mean())
plt.title(Clim_var_type + ' - SON Spring means/sums')
pd.DataFrame(Fdf_Seas_means[Fdf_Seas_means.index.quarter==4].mean()).plot(kind='bar', ylim=(ymin,ymax))
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
plt.title(Clim_var_type + ' - SON Spring means/sums')
Fdf_Seas_means[Fdf_Seas_means.index.quarter==4].boxplot(rot=90)
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
for Est in Estuaries:
#Estuary = 'HUNTER' # 'Belongil'
Estuary = Est # 'Belongil'
print Estuary
#Clim_var_type = 'potevpmean' # '*' will create pdf for all variables in folder "pracc*|tasmax*"
Clim_var_types = ['tasmax']
for climvar in Clim_var_types:
Clim_var_type = climvar
plot_pdf = 'no'
delta_csv = 'yes'
Stats = 'days_h_35' # 'maxdaily', 'mean'
Version = 'V1'
#==========================================================#
#==========================================================#
#set directory path for output files
output_directory = 'Output/' + Case_Study_Name + '/' + Estuary
#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('-------------------------------------------')
#==========================================================#
#==========================================================#
Estuary_Folder = glob.glob('./Data/NARCLIM_Site_CSVs/'+ Case_Study_Name + '/' + Estuary + '*' )
Clim_Var_CSVs = glob.glob(Estuary_Folder[0] + '/' + Clim_var_type + '*')
#==========================================================#
#==========================================================#
#read CSV files and start analysis
#==========================================================#
#for clim_var_csv_path in Clim_Var_CSVs:
clim_var_csv_path = Clim_Var_CSVs[0]
Filename = os.path.basename(os.path.normpath(clim_var_csv_path))
Clim_var_type = Filename.split('_', 1)[0]
print(clim_var_csv_path)
Full_df = pd.read_csv(clim_var_csv_path, index_col=0, parse_dates = True)
#pandas datestamp index to period (we don't need the 12 pm info in the index (daily periods are enough))
Full_df.index = Full_df.index.to_period('D')
Full_df = Full_df.drop(columns=['period'])
Ncols_df = len(Full_df)
#check data types of columns
#Full_df.dtypes
#==========================================================#
#==========================================================#
#substract a constant from all values to convert from kelvin to celcius (temp)
if Clim_var_type in ['tasmean','tasmax','sstmean']:
Full_df = Full_df.iloc[:,0:(Ncols_df-1)]-273.15
if Clim_var_type == 'evspsblmean' or Clim_var_type == 'potevpmean':
Full_df = Full_df.iloc[:,0:(Ncols_df-1)]*60*60*24
Fdf_1900_2080 = Full_df
if Clim_var_type in ['rsdsmean','rldsmean']:
Full_df = Full_df.iloc[:,0:(Ncols_df-1)]*60*60*24/1000000
#==========================================================#
#==========================================================#
#Aggregate daily df to annual time series
if Clim_var_type in ['pracc' ,'evspsblmean' ,'potevpmean' ,'pr1Hmaxtstep' ,
'wss1Hmaxtstep', 'rsdsmean', 'rldsmean']:
if(Stats == 'maxdaily'):
Fdf_1900_2080_annual = Fdf_1900_2080.resample('A').max()
Fdf_1900_2080_annual = Fdf_1900_2080_annual.replace(0, np.nan)
Fdf_Seas_means = Fdf_1900_2080.resample('Q-NOV').max() #seasonal means
Fdf_Seas_means = Fdf_Seas_means.replace(0, np.nan)
else:
Fdf_1900_2080_annual = Fdf_1900_2080.resample('A').sum()
Fdf_1900_2080_annual = Fdf_1900_2080_annual.replace(0, np.nan)
Fdf_Seas_means = Fdf_1900_2080.resample('Q-NOV').sum() #seasonal means
Fdf_Seas_means = Fdf_Seas_means.replace(0, np.nan)
else:
if(Stats == 'maxdaily'):
Fdf_1900_2080_annual = Fdf_1900_2080.resample('A').max()
Fdf_1900_2080_annual = Fdf_1900_2080_annual.replace(0, np.nan)
Fdf_Seas_means = Fdf_1900_2080.resample('Q-NOV').max() #seasonal means
Fdf_Seas_means = Fdf_Seas_means.replace(0, np.nan)
if(Stats[:4] =='days'):
Threshold = int(Stats[-2:])
#agg = ('abobe_27_count', lambda x: x.gt(27).sum()), ('average', 'mean')
agg = ('>'+ str(Threshold) + '_count', lambda x: x.gt(Threshold).sum()),
Fdf_1900_2080_annual = Fdf_1900_2080.resample('A').agg(agg)
#Fdf_1900_2080_annual = Fdf_1900_2080_annual.replace(0, np.nan)
Fdf_Seas_means = Fdf_1900_2080.resample('Q-NOV').agg(agg) #seasonal means
#Fdf_Seas_means = Fdf_Seas_means.replace(0, np.nan)
else:
Fdf_1900_2080_annual = Fdf_1900_2080.resample('A').mean()
Fdf_Seas_means = Fdf_1900_2080.resample('Q-NOV').mean() #seasonal means
Fdf_1900_2080_means = Fdf_1900_2080.mean()
#==========================================================#
#==========================================================#
#Select the 3 most representative models (min med and max difference betwen far future and present)
dfall, dfmin, dfmax, dfmed, Min_dif_mod_name, Med_dif_mod_name, Max_dif_mod_name = fct.select_min_med_max_dif_model(Fdf_1900_2080)
#==========================================================#
#==========================================================#
#create a dataframe that has 1 column for each of the three representative models
dfa = Fdf_1900_2080_annual.iloc[:,[0]]
dfa1 = Fdf_1900_2080_annual.iloc[:,[0,3,6]].loc[(Fdf_1900_2080_annual.index >= '1990') & (Fdf_1900_2080_annual.index <= '2009')]
dfa1.columns = [Min_dif_mod_name[:-5], Med_dif_mod_name[:-5], Max_dif_mod_name[:-5]]
dfa2 = Fdf_1900_2080_annual.iloc[:,[1,4,7]].loc[(Fdf_1900_2080_annual.index >= '2020') & (Fdf_1900_2080_annual.index <= '2039')]
dfa2.columns = [Min_dif_mod_name[:-5], Med_dif_mod_name[:-5], Max_dif_mod_name[:-5]]
dfa3 = Fdf_1900_2080_annual.iloc[:,[2,5,8]].loc[(Fdf_1900_2080_annual.index >= '2060') & (Fdf_1900_2080_annual.index <= '2079')]
dfa3.columns = [Min_dif_mod_name[:-5], Med_dif_mod_name[:-5], Max_dif_mod_name[:-5]]
dfall_annual = dfa1.append(dfa2).append(dfa3)
#==========================================================#
#==========================================================#
#Create Deltas of average change for annual and seasonal basis
#==========================================================#
delta_all_df = fct.calculate_deltas_NF_FF2(Fdf_1900_2080_annual, Fdf_Seas_means, Stats)
#==========================================================#
if delta_csv == 'yes':
out_file_name = Estuary + '_' + Clim_var_type + '_' + Stats + '_NARCliM_ensemble_changes.csv'
out_path = output_directory + '/' + out_file_name
delta_all_df.to_csv(out_path)
#==========================================================#
#==========================================================#
#create a dataframe that has a single column for present day, near and far future for the (3 selected models)
Full_current_df = Fdf_1900_2080.iloc[:,range(0,3)]
Full_current_df = Full_current_df.stack()
#nearfuture
Full_nearfuture_df = Fdf_1900_2080.iloc[:,range(3,6)]
Full_nearfuture_df = Full_nearfuture_df.stack()
#farfuture
Full_farfuture_df = Fdf_1900_2080.iloc[:,range(6,len(Fdf_1900_2080.columns))]
Full_farfuture_df = Full_farfuture_df.stack()
Summarized_df = pd.concat([Full_current_df, Full_nearfuture_df], axis=1, ignore_index=True)
Summarized_df = pd.concat([Summarized_df, Full_farfuture_df], axis=1, ignore_index=True)
Summarized_df.columns = ['present', 'near', 'far']
#==========================================================#
#==========================================================#
#output some summary plot into pdf
#==========================================================#
if plot_pdf == 'yes':
plotcolours36 = ['darkolivegreen','turquoise', 'lightgreen', 'darkgreen', 'lightpink','slateblue', 'slategray', 'orange', 'tomato', 'peru', 'navy', 'teal',
'darkolivegreen','turquoise', 'lightgreen', 'darkgreen', 'lightpink','slateblue', 'slategray', 'orange', 'tomato', 'peru', 'navy', 'teal',
'darkolivegreen','turquoise', 'lightgreen', 'darkgreen', 'lightpink','slateblue', 'slategray', 'orange', 'tomato', 'peru', 'navy', 'teal']
plotcolours36b = ['tomato', 'royalblue', 'mediumpurple' , 'tomato', 'royalblue', 'mediumpurple' , 'tomato', 'royalblue', 'mediumpurple' , 'tomato', 'royalblue', 'mediumpurple' ,
'tomato', 'royalblue', 'mediumpurple' , 'tomato', 'royalblue', 'mediumpurple' , 'tomato', 'royalblue', 'mediumpurple' , 'tomato', 'royalblue', 'mediumpurple' ,
'tomato', 'royalblue', 'mediumpurple' , 'tomato', 'royalblue', 'mediumpurple' , 'tomato', 'royalblue', 'mediumpurple' , 'tomato', 'royalblue', 'mediumpurple' ]
plotcolours12 = ['darkolivegreen','turquoise', 'lightgreen', 'darkgreen', 'lightpink','slateblue', 'slategray', 'orange', 'tomato', 'peru', 'navy', 'teal']
plotcolours15 = ['darkolivegreen','turquoise', 'lightgreen', 'darkgreen', 'lightpink','slateblue', 'slategray', 'orange', 'tomato', 'peru', 'navy', 'teal', 'lightgreen','lightpink','slateblue']
#plt.cm.Paired(np.arange(len(Fdf_1900_2080_means)))
#write the key plots to a single pdf document
pdf_out_file_name = Clim_var_type + '_' + Stats + '_start_' + Base_period_start + '_NARCliM_summary_' + Version + '.pdf'
pdf_out_path = output_directory +'/' + pdf_out_file_name
#open pdf and add the plots
with PdfPages(pdf_out_path) as pdf:
#barplot of model means
plt.title(Clim_var_type + ' - model means - full period')
ymin = min(Fdf_1900_2080_means)
ymax = max(Fdf_1900_2080_means) + 0.008 *min(Fdf_1900_2080_means)
Fdf_1900_2080_means.plot(kind='bar', ylim=(ymin,ymax), color=plotcolours36)
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
#
neardeltadf=delta_all_df['near']
ymin = min(neardeltadf) + 0.1 *min(neardeltadf)
ymax = max(neardeltadf) + 0.1 * max(neardeltadf)
neardeltadf=delta_all_df['far']
ymin2 = min(neardeltadf) + 0.1 *min(neardeltadf)
ymax2 = max(neardeltadf) + 0.1 * max(neardeltadf)
ymin = min(ymin, ymin2)
if (Clim_var_type == 'tasmax' or Clim_var_type == 'tasmean'):
ymin = 0
ymax = max(ymax, ymax2)
#
# delta barplot for report 1#################################
ax=plt.subplot(2,1,1)
plt.title(Clim_var_type + ' - model deltas - near-present')
neardeltadf=delta_all_df['near']
neardeltadf.plot(kind='bar', color=plotcolours15, ylim=(ymin,ymax), ax=ax)
plt.xticks([])
#ax.xaxis.set_ticklabels([])
#pdf.savefig(bbox_inches='tight', ylim=(ymin,ymax), pad_inches=0.4)
#plt.close()
#
ax=plt.subplot(2,1,2)
plt.title(Clim_var_type + ' - model deltas - far-present')
neardeltadf=delta_all_df['far']
neardeltadf.plot(kind='bar', color=plotcolours15, ylim=(ymin,ymax), ax=ax)
ax.xaxis.grid(False)
#fig.patch.set_alpha(0)
#plt.show()
pdf.savefig(bbox_inches='tight', ylim=(ymin,ymax), pad_inches=0.4)
plt.close()
# end delta barplot for report 1#################################
#
#full period density comparison
plt.title(Clim_var_type + ' - density comparison - full period - all models')
Summarized_df.plot.kde()
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
#full period density comparison
plt.title(Clim_var_type + ' - density comparison - full period - max delta model')
xmin = float(max(np.nanpercentile(Fdf_1900_2080.filter(regex= Max_dif_mod_name[:-5]),50) - 4 * np.std(Fdf_1900_2080.filter(regex= Max_dif_mod_name[:-5]))))
xmax = float(max(np.nanpercentile(Fdf_1900_2080.filter(regex= Max_dif_mod_name[:-5]),50) + 4 * np.std(Fdf_1900_2080.filter(regex= Max_dif_mod_name[:-5]))))
Fdf_1900_2080.filter(regex= Max_dif_mod_name[:-5]).plot.kde(xlim=(xmin,xmax))
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
#annual box
plt.title(Clim_var_type + ' - Annual means/sums for max diff model')
Fdf_1900_2080_annual.boxplot(rot=90)
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
#
#daily box
plt.title(Clim_var_type + ' - Daily means/sums')
Fdf_1900_2080.boxplot(rot=90)
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
# time series plot annual ALL models
plt.title(Clim_var_type + ' - Time series - all models')
Mod_order = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,19,20,21,16,17,18,22,23,24,31,32,33,25,26,27,28,29,30,34,35,36]
test = Fdf_1900_2080_annual
Mod_Names = test.columns
New_Mod_Name = []
for i in range(0,len(Mod_Names)):
New_Mod_Name.append(str(Mod_order[i]+10) + '_' + Mod_Names[i])
test.columns = New_Mod_Name
test_sorted = test.reindex_axis(sorted(test.columns), axis=1)
colnamest = test.columns
test_sorted.columns = [w[3:-5] for w in colnamest]
test_sorted.plot(legend=False, color = plotcolours36)
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
# time series plot annual ALL models
plt.title(Clim_var_type + ' - Time series - representative models')
dfall_annual.plot(legend=False)
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
# seasonal mean boxplots
ymin = min(Fdf_Seas_means[Fdf_Seas_means.index.quarter==1].mean())
ymax = max(Fdf_Seas_means[Fdf_Seas_means.index.quarter==1].mean())
plt.title(Clim_var_type + ' - DJF Summer means/sums')
pd.DataFrame(Fdf_Seas_means[Fdf_Seas_means.index.quarter==1].mean()).plot(kind='bar', ylim=(ymin,ymax))
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
plt.title(Clim_var_type + ' - DJF Summer means/sums')
Fdf_Seas_means[Fdf_Seas_means.index.quarter==1].boxplot(rot=90)
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
ymin = min(Fdf_Seas_means[Fdf_Seas_means.index.quarter==2].mean())
ymax = max(Fdf_Seas_means[Fdf_Seas_means.index.quarter==2].mean())
plt.title(Clim_var_type + ' - MAM Autumn means/sums')
pd.DataFrame(Fdf_Seas_means[Fdf_Seas_means.index.quarter==2].mean()).plot(kind='bar', ylim=(ymin,ymax))
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
plt.title(Clim_var_type + ' - MAM Autumn means/sums')
Fdf_Seas_means[Fdf_Seas_means.index.quarter==2].boxplot(rot=90)
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
ymin = min(Fdf_Seas_means[Fdf_Seas_means.index.quarter==3].mean())
ymax = max(Fdf_Seas_means[Fdf_Seas_means.index.quarter==3].mean())
plt.title(Clim_var_type + ' - JJA Winter means/sums')
pd.DataFrame(Fdf_Seas_means[Fdf_Seas_means.index.quarter==3].mean()).plot(kind='bar', ylim=(ymin,ymax))
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
plt.title(Clim_var_type + ' - JJA Winter means/sums')
Fdf_Seas_means[Fdf_Seas_means.index.quarter==3].boxplot(rot=90)
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
ymin = min(Fdf_Seas_means[Fdf_Seas_means.index.quarter==4].mean())
ymax = max(Fdf_Seas_means[Fdf_Seas_means.index.quarter==4].mean())
plt.title(Clim_var_type + ' - SON Spring means/sums')
pd.DataFrame(Fdf_Seas_means[Fdf_Seas_means.index.quarter==4].mean()).plot(kind='bar', ylim=(ymin,ymax))
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
plt.title(Clim_var_type + ' - SON Spring means/sums')
Fdf_Seas_means[Fdf_Seas_means.index.quarter==4].boxplot(rot=90)
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()

5
Analysis/Code/P1_SILO_Batch_Download.py
Unescape Escape View File

@ -38,12 +38,11 @@ os.chdir('C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdo
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', 'evap_syn']
Location = 'Catchment'
Location = 'Estuary' #'Catchment'
startdate= '19600101'
enddate= '20180101'
#==========================================================#
#==========================================================#
#set directory path for output files
output_directory = 'Data/SILO/' + Case_Study_Name + '/'
@ -85,7 +84,7 @@ for Estuary in mydict:
print('-------------------------------------------')
#==========================================================#
output_csv = output_directory_internal + 'SILO_climdata_' + Estuary +'_'+ Location +'_' + mydict[Estuary][0] + '_' + mydict[Estuary][1] + '2.csv'
output_csv = output_directory_internal + 'SILO_climdata_' + Estuary +'_'+ Location +'_' + mydict[Estuary][0] + '_' + mydict[Estuary][1] + '.csv'
silo_df = sil.pointdata(Silo_variables, 'Okl9EDxgS2uzjLWtVNIBM5YqwvVcCxOmpd3nCzJh',startdate, enddate,
None, mydict[Estuary][0], mydict[Estuary][1], True, output_csv)
#==========================================================#

98
Analysis/Code/climdata_fcts.py
Unescape Escape View File

@ -9,7 +9,7 @@ import matplotlib.pyplot as plt
from datetime import datetime, timedelta
import numpy as np
import pandas as pd
import glob
def compare_images(im1, im2):
@ -40,8 +40,6 @@ def datenum2datetime(datenum):
return time
def select_min_med_max_dif_model(NARCLIM_df):
#Select the 3 most representative models (min med and max difference betwen far future and present)
Fdf_1900_2080_sorted = NARCLIM_df.reindex_axis(sorted(NARCLIM_df.columns), axis=1)
@ -73,7 +71,7 @@ def select_min_med_max_dif_model(NARCLIM_df):
return dfall , dfmin, dfmed, dfmax, Min_dif_mod_name,Med_dif_mod_name, Max_dif_mod_name
def calculate_deltas_NF_FF2(Annual_df, Seasonal_df):
def calculate_deltas_NF_FF2(Annual_df, Seasonal_df, Stats):
"""calculates the "deltas" between nearfuture and present day for annual or seasonal climate data in pandas TS format"""
times = ['annual', 'DJF', 'MAM', 'JJA','SON']
@ -94,9 +92,11 @@ def calculate_deltas_NF_FF2(Annual_df, Seasonal_df):
if temp == 'SON':
Mean_df = Seasonal_df[Seasonal_df.index.quarter==4].mean()
Column_names = ['SON_near', 'SON_far']
models = list(Seasonal_df.mean().index)
if(Stats[:4] =='days'):
models = list(Seasonal_df.mean().index.get_level_values(0))
else:
models = list(Seasonal_df.mean().index)
newmodel = []
type(newmodel)
for each in models:
newmodel.append(each[:-5])
unique_models = set(newmodel)
@ -123,4 +123,88 @@ def calculate_deltas_NF_FF2(Annual_df, Seasonal_df):
delta_df.loc['90th'] = pd.Series({Column_names[0]:np.percentile(delta_df[Column_names[0]], 90), Column_names[1]:np.percentile(delta_df[Column_names[1]], 90)})
#append df to overall df
delta_all_df = pd.concat([delta_all_df, delta_df], axis=1)
return delta_all_df
if(Stats[:4] =='days'):
delta_all_df = delta_all_df .astype(int).fillna(0.0)
return delta_all_df
def import_present_day_climdata_csv(Estuary, Clim_var_type):
"""
this funciton imports the present day climate data used for
characterizing the present day climate varibility
If DataSource == 'Station', individual weather station data is used.
If DataSource == 'SILO' , SILO time series is used using the estuary centerpoint as reference locatoin for
selection of the grid cell
"""
#load present day climate data for the same variable
if Clim_var_type == 'evspsblmean': #ET time series that we have is not in the same format as the other variables, hence the different treatment
Present_day_Var_CSV = glob.glob('./Data/Wheather_Station_Data/**/' + Estuary + '_' + 'ET' + '*csv')
Present_day_df = pd.read_csv(Present_day_Var_CSV[0])
Dates = pd.to_datetime(Present_day_df.Date)
Present_day_df.index = Dates
Present_day_df = Present_day_df.iloc[:,1]
Present_day_df = Present_day_df.replace(r'\s+', np.nan, regex=True)
Present_day_df = pd.to_numeric(Present_day_df)
Present_day_df.index = Dates
[minplotDelta, maxplotDelta]=[50,50]
#for tasmean, observed min and max T need to be converted into mean T
elif Clim_var_type == 'tasmean':
Present_day_Var_CSV = glob.glob('./Data/Wheather_Station_Data/**/' + Estuary + '_MaxT*csv')
Present_day_df = pd.read_csv(Present_day_Var_CSV[0])
Dates = pd.to_datetime(Present_day_df.Year*10000+Present_day_df.Month*100+Present_day_df.Day,format='%Y%m%d')
Present_day_df.index = Dates
Present_day_MaxT_df = Present_day_df.iloc[:,5]
Present_day_Var_CSV = glob.glob('./Data/Wheather_Station_Data/**/' + Estuary + '_MinT*csv')
Present_day_df = pd.read_csv(Present_day_Var_CSV[0])
Dates = pd.to_datetime(Present_day_df.Year*10000+Present_day_df.Month*100+Present_day_df.Day,format='%Y%m%d')
Present_day_df.index = Dates
Present_day_MinT_df = Present_day_df.iloc[:,5]
Present_day_df = (Present_day_MaxT_df + Present_day_MinT_df)/2
[minplotDelta, maxplotDelta]=[1,2]
elif Clim_var_type == 'tasmax':
Present_day_Var_CSV = glob.glob('./Data/Wheather_Station_Data/**/' + Estuary + '_MaxT*csv')
Present_day_df = pd.read_csv(Present_day_Var_CSV[0])
Dates = pd.to_datetime(Present_day_df.Year*10000+Present_day_df.Month*100+Present_day_df.Day,format='%Y%m%d')
Present_day_df.index = Dates
Present_day_df = Present_day_df.iloc[:,5]
[minplotDelta, maxplotDelta]=[1,2]
elif Clim_var_type == 'wssmean' or Clim_var_type == 'wss1Hmaxtstep':
Present_day_Var_CSV = glob.glob('./Data/Wheather_Station_Data/**/Terrigal_Wind.csv')
Present_day_df = pd.read_csv(Present_day_Var_CSV[0])
Present_day_df.index = Present_day_df[['Year', 'Month', 'Day', 'Hour']].apply(lambda s : datetime(*s),axis = 1)
Present_day_df = Present_day_df.filter(regex= 'm/s')
Present_day_df = Present_day_df.replace(r'\s+', np.nan, regex=True)
Present_day_df['Wind speed measured in m/s'] = Present_day_df['Wind speed measured in m/s'].convert_objects(convert_numeric=True)
[minplotDelta, maxplotDelta]=[1, 1.5]
elif Clim_var_type == 'sstmean':
Estuary_Folder = glob.glob('./Data/NARCLIM_Site_CSVs/CASESTUDY2/' + Estuary + '*' )
Present_day_Var_CSV = glob.glob(Estuary_Folder[0] + '/' + Clim_var_type + '_NNRP*')
Present_day_df = pd.read_csv(Present_day_Var_CSV[0], parse_dates=True, index_col=0)
Present_day_df = Present_day_df.filter(regex= 'NNRP_R1_1950')
Present_day_df['NNRP_R1_1950'] = Present_day_df['NNRP_R1_1950'].convert_objects(convert_numeric=True)
[minplotDelta, maxplotDelta]=[1, 1]
else:
Present_day_Var_CSV = glob.glob('./Data/Wheather_Station_Data/**/' + Estuary + '_' + 'Rainfall' + '*csv')
Present_day_df = pd.read_csv(Present_day_Var_CSV[0])
Dates = pd.to_datetime(Present_day_df.Year*10000+Present_day_df.Month*100+Present_day_df.Day,format='%Y%m%d')
Present_day_df.index = Dates
Present_day_df = Present_day_df.iloc[:,5]
[minplotDelta, maxplotDelta]=[50,100]
return Present_day_df, minplotDelta, maxplotDelta

165
Analysis/Code/ggplot_time_series_with_trends_Waves_SST_OEH.R
Unescape Escape View File

@ -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()

45
Analysis/Code/ggplot_time_series_with_trends_sea_levle.R → Analysis/Code/ggplot_time_series_with_trends_sea_level.R
Unescape Escape View File

@ -19,41 +19,64 @@ setwd("C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdoc/P
#Set inputs
######################
Case.Study <- "CASESTUDY2"
Estuary <- "HUNTER"
Estuary <- "NADGEE"
Climvar <- 'tasmean'
ggplotGAM.k <- 7
######################
######################
#Set input file paths
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')
}
######################
AirT_CSV_Path <- "./Data/Ocean_Data/BOM_monthly_SL_Hunter_Newcastle.txt"
######################
#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,]
dat2 = dat[-c(1:3), ]
dat2 = dat2[-(720:726),]
dat2 = dat2[-(1:108),]
dat2 = dat2[,-1]
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
######################
@ -68,7 +91,7 @@ p1air <- ggplot(SeaLev.df, aes(y=MSL, x=Date)) + geom_line(alpha=0.5) +
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.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()
@ -84,7 +107,7 @@ p1air <- ggplot(SeaLev.df, aes(y=MSL, x=Date)) + geom_line(alpha=0.5) +
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.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()
Write Preview
Loading…
Cancel
Save