diff --git a/Analysis/Code/NARCLIM_Download/R_Preparing_BASH_script_for_NARCLIM_batch_download.R b/Analysis/Code/NARCLIM_Download/R_Preparing_BASH_script_for_NARCLIM_batch_download.R index 177ad13..12a905c 100644 --- a/Analysis/Code/NARCLIM_Download/R_Preparing_BASH_script_for_NARCLIM_batch_download.R +++ b/Analysis/Code/NARCLIM_Download/R_Preparing_BASH_script_for_NARCLIM_batch_download.R @@ -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")) diff --git a/Analysis/Code/P1_NARCliM_First_Pass_variab_deviation_plots.py b/Analysis/Code/P1_NARCliM_First_Pass_variab_deviation_plots.py index 6bb8770..6788090 100644 --- a/Analysis/Code/P1_NARCliM_First_Pass_variab_deviation_plots.py +++ b/Analysis/Code/P1_NARCliM_First_Pass_variab_deviation_plots.py @@ -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 + + diff --git a/Analysis/Code/P1_NARCliM_First_Pass_variab_deviation_plots_nonNARCLIM.py b/Analysis/Code/P1_NARCliM_First_Pass_variab_deviation_plots_nonNARCLIM.py new file mode 100644 index 0000000..7113eff --- /dev/null +++ b/Analysis/Code/P1_NARCliM_First_Pass_variab_deviation_plots_nonNARCLIM.py @@ -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) + + + + + + + + + + \ No newline at end of file diff --git a/Analysis/Code/P1_NARCliM_NC_to_CSV_CCRC_SS_BASH_script_readme.txt b/Analysis/Code/P1_NARCliM_NC_to_CSV_CCRC_SS_BASH_script_readme.txt index 6f9fb41..55dc031 100644 --- a/Analysis/Code/P1_NARCliM_NC_to_CSV_CCRC_SS_BASH_script_readme.txt +++ b/Analysis/Code/P1_NARCliM_NC_to_CSV_CCRC_SS_BASH_script_readme.txt @@ -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: diff --git a/Analysis/Code/P1_NARCliM_plots_Windows.py b/Analysis/Code/P1_NARCliM_plots_Windows.py index c1e0c69..175a12b 100644 --- a/Analysis/Code/P1_NARCliM_plots_Windows.py +++ b/Analysis/Code/P1_NARCliM_plots_Windows.py @@ -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() \ No newline at end of file +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() \ No newline at end of file diff --git a/Analysis/Code/P1_SILO_Batch_Download.py b/Analysis/Code/P1_SILO_Batch_Download.py index dfb5379..44c35f4 100644 --- a/Analysis/Code/P1_SILO_Batch_Download.py +++ b/Analysis/Code/P1_SILO_Batch_Download.py @@ -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) #==========================================================# diff --git a/Analysis/Code/climdata_fcts.py b/Analysis/Code/climdata_fcts.py index b81a2d4..a65ad5e 100644 --- a/Analysis/Code/climdata_fcts.py +++ b/Analysis/Code/climdata_fcts.py @@ -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 \ No newline at end of file + 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 + + + + + + + + + + + + + + + + diff --git a/Analysis/Code/ggplot_time_series_with_trends_sea_levle.R b/Analysis/Code/ggplot_time_series_with_trends_sea_levle.R deleted file mode 100644 index 58eee02..0000000 --- a/Analysis/Code/ggplot_time_series_with_trends_sea_levle.R +++ /dev/null @@ -1,90 +0,0 @@ -#R code for creating ggplots of time series with smooth (GAM) and linear term - - -###################### -#Import Libraries and set working directory -###################### -library(zoo) -library(hydroTSM) #you need to install these packages first before you can load them here -library(lubridate) -library(mgcv) -library(ggplot2) -library(gridExtra) -library(scales) -options(scipen=999) -setwd("C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdoc/Projects/Paper#1/") -###################### - -###################### -#Set inputs -###################### -Case.Study <- "CASESTUDY2" -Estuary <- "HUNTER" -Climvar <- 'tasmean' -ggplotGAM.k <- 7 -###################### - - -###################### -#Set input file paths -###################### -AirT_CSV_Path <- "./Data/Ocean_Data/BOM_monthly_SL_Hunter_Newcastle.txt" - - - - - - -dat = readLines(AirT_CSV_Path) -dat = as.data.frame(do.call(rbind, strsplit(dat, split=" {2,10}")), stringsAsFactors=FALSE) -colnames(dat) <-dat[3,] - -dat2 = dat[-c(1:3), ] -dat2 = dat2[-(720:726),] -dat2 = dat2[-(1:108),] -dat2 = dat2[,-1] - -dat2$Date <- as.yearmon(dat2[,1], "%m %Y") - -SeaLev.df <- dat2 -head(SeaLev.df) -SeaLev.df$MSL <- as.numeric(SeaLev.df$Mean) -SeaLev.df$Julday1 <- seq(1,length(SeaLev.df[,1]),1) -linear.trend.model_EC_all <- lm(MSL ~ Julday1, SeaLev.df) -SeaLev.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4] -SeaLev.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 12 - -###################### -#Plot -###################### - -##################################### Full Time Period -p1air <- ggplot(SeaLev.df, aes(y=MSL, x=Date)) + geom_line(alpha=0.5) + - ggtitle(paste(Estuary, " - Linear and smooth trend in monthly mean sea level (BOM Gauge) | lin trend was ", - round(100*SeaLev.lintrend,3), ' cm/year with p=', round(SeaLev.pvalNCV_ECall,10), sep=" ")) + - theme(plot.title=element_text(face="bold", size=9)) + - geom_smooth(method='lm',fill="green", formula=y~x, colour="darkgreen", size = 0.5) + - stat_smooth(method=gam, formula=y~s(x, k=4), se=T, size=0.5, col="red") + - ylab("Monthly Mean Sea Level [m]") + xlab("Time") - -#export to png -png.name <- paste('./Output/', Case.Study, '/', Estuary, '/Trends_MonthlyMeanSeaLevel_full_period_', Sys.Date(),".png", sep="") -png(file = png.name, width = 10.5, height = 7, units='in', res=500) -grid.arrange(p1air,ncol=1) -dev.off() - -#multiple smooths -p1air <- ggplot(SeaLev.df, aes(y=MSL, x=Date)) + geom_line(alpha=0.5) + - ggtitle(paste(Estuary, " - Linear and smooth trend in monthly mean sea level (BOM Gauge) | lin trend was ", - round(100*SeaLev.lintrend,3), ' cm/year with p=', round(SeaLev.pvalNCV_ECall,10), sep=" ")) + - theme(plot.title=element_text(face="bold", size=9)) + - stat_smooth(method=gam, formula=y~s(x, k=13), se=T, size=0.5, col="red") + - #stat_smooth(method=gam, formula=y~s(x, k=8), se=T, size=0.5, cor="blue") + - stat_smooth(method=gam, formula=y~s(x, k=5), se=T, size=0.5, col="green") + - ylab("Monthly Mean Sea Level [m]") + xlab("Time") - -#export to png -png.name <- paste('./Output/', Case.Study, '/', Estuary, '/Trends_MonthlyMeanSeaLevel_MultiGAM_full_period_', Sys.Date(),".png", sep="") -png(file = png.name, width = 10.5, height = 7, units='in', res=500) -grid.arrange(p1air,ncol=1) -dev.off()