#added new code for also deriving quantile scaling factors from the NARCLIM data

Analysis/Code/P1_NARCliM_plots_Windows.py

@@ -50,7 +50,7 @@ for Est in Estuaries:
     Estuary = Est # 'Belongil'
     print Estuary
     #Clim_var_type  = 'potevpmean'    #  '*' will create pdf for all variables in folder  "pracc*|tasmax*"   
-    Clim_var_types  = ['pracc']
+    Clim_var_types  = ['evspsblmean']
     for climvar in Clim_var_types:
         Clim_var_type  = climvar 
         
@@ -91,6 +91,8 @@ for Est in Estuaries:
         if Location == 'Catchment':
             Estuary_Folder = glob.glob('./Data/NARCLIM_Site_CSVs/'+ Case_Study_Name + '/' + Estuary + '*catch' )
         Clim_Var_CSVs = glob.glob(Estuary_Folder[0] + '/' + Clim_var_type + '*')
+        Coordinates_foldername = os.path.basename(os.path.normpath(Estuary_Folder[0])) 
+        Coordinates = Coordinates_foldername.split('_', 3)[1] + '_' +Coordinates_foldername.split('_', 3)[2]
         #==========================================================#        
         
         #==========================================================#
@@ -104,7 +106,7 @@ for Est in Estuaries:
         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'])
+        Full_df = Full_df.drop(columns=['period'])
         Ncols_df = len(Full_df)
         #check data types of columns
         #Full_df.dtypes
@@ -213,11 +215,11 @@ for Est in Estuaries:
             
         #==========================================================#
         if delta_csv == 'yes':
-            out_file_name = Estuary  + '_' + Clim_var_type + '_' + Stats + '_NARCliM_ensemble_changes.csv'
+            out_file_name = Estuary  + '_'  + Clim_var_type + '_' + Stats + '_NARCliM_ensemble_changes' + Perc_Vs_Abs +  '_' + Location + '_' + Coordinates + '.csv'
             out_path = output_directory + '/' + out_file_name
             delta_all_df.to_csv(out_path)
         if delta_csv == 'yes' and DoMonthly ==True:
-            out_file_name = Estuary  + '_' + Clim_var_type + '_' + Stats + '_NARCliM_ensemble_changes_monthly.csv'
+            out_file_name = Estuary  + '_' + Clim_var_type + '_' + Stats + '_NARCliM_ensemble_changes_monthly' + Perc_Vs_Abs + '_' + Location + '_' + Coordinates + '.csv'
             out_path = output_directory + '/NARCLIM_Changes_Hunter/' + out_file_name
             delta_all_df_monthly.to_csv(out_path)
         #==========================================================#

Analysis/Code/P1_NARCliM_plots_Windows_quantile_scaling_factors.py

@@ -0,0 +1,520 @@
+# -*- coding: utf-8 -*-
+#==========================================================#
+#Last Updated - June 2018
+#@author: z5025317 Valentin Heimhuber
+#code for creating climate prioritization plots for NARCLIM variables. 
+#Inputs: Uses CSV files that contain all 12 NARCLIM model runs time series for 1 grid cell created with: P1_NARCliM_NC_to_CSV_CCRC_SS.py
+
+#==========================================================#
+#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 *
+matplotlib.style.use('ggplot')
+
+# import own modules
+# Set working direcotry (where postprocessed NARClIM data is located)
+os.chdir('C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdoc/Projects/Paper#1/Analysis/Code')
+import climdata_fcts as fct
+import silo as sil
+
+ALPHA_figs = 0 
+font = {'family' : 'sans-serif',
+        'weight' : 'normal',
+        'size'   : 14}
+matplotlib.rc('font', **font)
+#==========================================================#
+# 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'
+Version = 'V5' #V4 is latest for all model runs #V5 is just for the Hunter catchment climatology
+
+Estuaries = ['HUNTER', 'RICHMOND', 'NADGEE', 'SHOALHAVEN', 'GEORGES','CATHIE']
+Estuaries = ['HUNTER'] #,'HUNTER','RICHMOND']
+#Estuaries = ['SHOALHAVEN']
+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  = ['pracc']
+    for climvar in Clim_var_types:
+        Clim_var_type  = climvar 
+        
+        #==========================================================#
+        #set input preferences
+        #==========================================================#
+        plot_pdf = 'no'
+        plot_pngs = 'no'
+        delta_csv = 'yes'
+        Stats = 'mean'   # 'maxdaily', 'mean', 'days_h_25'
+        DoMonthly = True
+        Perc_Vs_Abs = 'percent' #'percent' vs. 'absolute' deltas for near and far future
+        Location = 'Catchment'
+        #==========================================================# 
+        
+        #==========================================================#
+        #set directory path for output files
+        output_directory = 'Output/' + Case_Study_Name + '_' + Version + '/' + 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('-------------------------------------------')
+            
+        #set directory path for individual png files
+        png_output_directory = 'Output/' + Case_Study_Name +  '_' + Version + '/' + Estuary + '/NARCLIM_Key_Figs'
+        #output_directory = 'J:/Project wrl2016032/NARCLIM_Raw_Data/Extracted'
+        if not os.path.exists(png_output_directory):
+            os.makedirs(png_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 + '*' )
+        if Location == 'Catchment':
+            Estuary_Folder = glob.glob('./Data/NARCLIM_Site_CSVs/'+ Case_Study_Name + '/' + Estuary + '*catch' )
+        Clim_Var_CSVs = glob.glob(Estuary_Folder[0] + '/' + Clim_var_type + '*')
+        Coordinates_foldername = os.path.basename(os.path.normpath(Estuary_Folder[0])) 
+        Coordinates = Coordinates_foldername.split('_', 3)[1] + '_' +Coordinates_foldername.split('_', 3)[2]
+        #==========================================================#        
+        
+        #==========================================================#
+        #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']) #not part of all input csv files
+        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
+        if Clim_var_type in ['pr30maxtstep','pr10maxtstep','pr20maxtstep']:
+            Full_df = Full_df.iloc[:,0:(Ncols_df-1)]*60*30 #int(Clim_var_type[2:4])
+        if Clim_var_type in ['pr1Hmaxtstep']:
+            Full_df = Full_df.iloc[:,0:(Ncols_df-1)]*60*60           
+        if Clim_var_type in ['rsdsmean','rldsmean']:
+            Full_df = Full_df.iloc[:,0:(Ncols_df-1)]*60*60*24/1000000  #convert to unit used in SILO
+        Fdf_1900_2080 = Full_df
+        #==========================================================#
+        
+        
+        
+        
+        
+        #==========================================================#
+        #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)
+                Fdf_Monthly_means = Fdf_1900_2080.resample('M').max() #seasonal means
+                Fdf_Monthly_means = Fdf_Monthly_means.replace(0, np.nan)
+            elif(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)
+                Fdf_Monthly_means = Fdf_1900_2080.resample('M').agg(agg)
+            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)
+                Fdf_Monthly_means = Fdf_1900_2080.resample('M').sum() #Monthlyonal means
+                Fdf_Monthly_means = Fdf_Monthly_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)
+                Fdf_Monthly_means = Fdf_1900_2080.resample('M').max() #Monthlyonal means
+                Fdf_Monthly_means = Fdf_Monthly_means.replace(0, np.nan)
+            elif(Stats == 'maxmonthly'):
+                Fdf_1900_2080_annual = Fdf_1900_2080.resample('M').max().resample('A').mean()
+                Fdf_1900_2080_annual = Fdf_1900_2080_annual.replace(0, np.nan)
+                Fdf_Seas_means = Fdf_1900_2080.resample('M').max().resample('Q-NOV').mean() #seasonal means
+                Fdf_Seas_means = Fdf_Seas_means.replace(0, np.nan)
+                Fdf_Monthly_means = Fdf_1900_2080.resample('M').max().resample('Q-NOV').mean() #Monthlyonal means
+                Fdf_Monthly_means = Fdf_Monthly_means.replace(0, np.nan)
+            elif(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)
+                Fdf_Monthly_means = Fdf_1900_2080.resample('M').agg(agg)
+            else:
+                Fdf_1900_2080_annual = Fdf_1900_2080.resample('A').mean()
+                Fdf_Seas_means = Fdf_1900_2080.resample('Q-NOV').mean() #seasonal means
+                Fdf_Monthly_means = Fdf_1900_2080.resample('M').mean() #seasonal means
+        Fdf_1900_2080_means = Fdf_1900_2080.mean()
+        #==========================================================#
+        
+        
+        #==========================================================#
+        # construction of quantile scaling method
+        #==========================================================#
+        quantiles = [0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,0.91,0.92,0.93,0.94,0.95,0.96,0.97,0.98,0.99,1]
+        Quantile_change_df = quant_quant_scaling(Full_df, quantiles, False)
+        
+        png_out_file_name = Clim_var_type + '_' + Stats + '_Deltas_Barplot_' + Version + '.png'
+        png_out_path = png_output_directory + '/' + png_out_file_name
+            
+        #==========================================================#
+        #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, Perc_Vs_Abs)
+        if DoMonthly ==True:
+            delta_all_df_monthly = fct.calculate_deltas_monthly(Fdf_Monthly_means, Stats, Perc_Vs_Abs)  
+            delta_all_df_monthly = pd.concat([delta_all_df_monthly.filter(regex= 'near'),delta_all_df_monthly.filter(regex= 'far')], axis=1)[-3:]  
+            
+        #==========================================================#
+        if delta_csv == 'yes':
+            out_file_name = Estuary  + '_'  + Clim_var_type + '_' + Stats + '_NARCliM_ensemble_changes_' + Perc_Vs_Abs +  '_' + Location + '_' + Coordinates + '.csv'
+            delta_all_df.to_csv(output_directory + '/' + out_file_name)
+            #quantile scaling csv
+            out_file_name = Estuary  + '_'  + Clim_var_type + '_' + Stats + '_NARCliM_quant_quant_' + Coordinates + '.csv'
+            Quantile_change_df.to_csv(output_directory + '/NARCLIM_Changes_Hunter/' + out_file_name)
+        if delta_csv == 'yes' and DoMonthly ==True:
+            out_file_name = Estuary  + '_' + Clim_var_type + '_' + Stats + '_NARCliM_ensemble_changes_monthly_' + Perc_Vs_Abs + '_' + Location + '_' + Coordinates + '.csv'
+            delta_all_df_monthly.to_csv(output_directory + '/NARCLIM_Changes_Hunter/' + out_file_name)
+        #==========================================================#
+        
+        
+        #==========================================================#       
+        pd.qcut(range(5), 4)
+        #==========================================================#        
+        
+        #==========================================================#
+        #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']
+        #==========================================================#
+        
+
+        #==========================================================#
+        #generate colour schemes for plotting
+        #==========================================================#
+        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']
+        color_dict = {'CCCMA3.1_R1': 'darkolivegreen', 'CCCMA3.1_R2': 'turquoise', 'CCCMA3.1_R3': 'lightgreen', 'CSIRO-MK3.0_R2': 'darkgreen',
+                  'CSIRO-MK3.0_R1':'lightpink', 'CSIRO-MK3.0_R3':'slateblue','ECHAM5_R1':'slategray', 'ECHAM5_R2': 'orange', 'ECHAM5_R3':  'tomato',
+                  'MIROC3.2_R1': 'peru', 'MIROC3.2_R2': 'navy' ,'MIROC3.2_R3': 'teal', '10th':'lightgreen','median':'lightpink','90th':'slateblue'}
+        #==========================================================#
+         
+         
+        #==========================================================#
+        #output crucial summary plots into individual png files
+        #==========================================================#
+        if plot_pngs == 'yes':
+            #=========#
+            #Barplot of near and far future deltas
+            #=========#
+            #out name
+            png_out_file_name = Clim_var_type + '_' + Stats + '_Deltas_Barplot_' + Version + '.png'
+            png_out_path = png_output_directory + '/' + png_out_file_name
+            #prepare data frame for plot
+            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)
+            #
+            fig = plt.figure(figsize=(5,6))
+            ax=plt.subplot(2,1,1) 
+            plt.title(Clim_var_type + ' - ' + Stats + ' - deltas - near')
+            neardeltadf=delta_all_df['near']
+            neardeltadf.plot(kind='bar', color=plotcolours15, ylim=(ymin,ymax), ax=ax)  
+            plt.xticks([])
+            ax=plt.subplot(2,1,2) 
+            plt.title(Clim_var_type + ' - ' + Stats + ' - deltas - far')
+            neardeltadf=delta_all_df['far']
+            neardeltadf.plot(kind='bar', color=plotcolours15, ylim=(ymin,ymax), ax=ax)
+            ax.xaxis.grid(False)
+            #ax.patch.set_alpha(ALPHA_figs)
+            fig.patch.set_alpha(ALPHA_figs)
+            fig.tight_layout()
+            fig.savefig(png_out_path)
+            plt.close()
+            #=========#
+            
+            #=========#
+            #full period density comparison
+            #=========#
+            #out name
+            png_out_file_name = Clim_var_type + '_' + Stats + '_MaxDeltaMod_Histogram_' + Version + '.png'
+            png_out_path = png_output_directory + '/' + png_out_file_name
+            plt.title(Clim_var_type + ' - ' + Stats + ' -  hist - full period - max delta model')
+            #prepare data
+            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]))))
+            fig = plt.figure(figsize=(5,5))
+            ax=plt.subplot(2,1,1) 
+            Fdf_1900_2080.filter(regex= Max_dif_mod_name[:-5]).plot.kde(xlim=(xmin,xmax), ax=ax)
+            plt.legend(loc=9, bbox_to_anchor=(0.5, -0.3))
+            #ax.xaxis.grid(False)
+            ax.yaxis.grid(False)
+            fig.patch.set_alpha(ALPHA_figs)
+            fig.tight_layout()
+            fig.savefig(png_out_path)
+            plt.close()
+            #=========#
+            
+            #=========#
+            # time series plot annual ALL models
+            #=========#
+            png_out_file_name = Clim_var_type + '_' + Stats + '_TimeSeries_AllMods_' + Version + '3.png'
+            png_out_path = png_output_directory + '/' + png_out_file_name
+            plt.title(Clim_var_type + ' - ' + Stats + ' -  Time series - representative models')
+            #prepare data
+#            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)):
+#                if(Stats[:4] =='days'):
+#                    New_Mod_Name.append(str(Mod_order[i]+10) + '_' +  Mod_Names[i][0])
+#                else:
+#                    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]
+
+            #plot
+            fig = plt.figure(figsize=(8,7))
+            ax=plt.subplot(2,1,1)
+            #test_sorted.plot(ax=ax) #,color =  plotcolours36)
+            Fdf_1900_2080_annual_b = Fdf_1900_2080_annual
+            New_Mod_Name = []
+            if(Stats[:4] =='days'):
+                for each in Fdf_1900_2080_annual_b.columns:
+                    New_Mod_Name.append(each[0][:-5])
+            else:
+                for each in Fdf_1900_2080_annual_b.columns:
+                    New_Mod_Name.append(each[:-5])
+            Fdf_1900_2080_annual_b.columns = New_Mod_Name
+            Fdf_1900_2080_annual_b.plot(ax=ax, color=[color_dict[x] for x in Fdf_1900_2080_annual_b.columns])
+            plt.legend(loc=9, bbox_to_anchor=(0.5, -0.2))
+            ax.xaxis.grid(False)
+            fig.patch.set_alpha(ALPHA_figs)
+            fig.tight_layout()
+            fig.savefig(png_out_path)
+            plt.close()
+            #=========#
+        
+        #==========================================================#
+        #output some summary plot into pdf
+        #==========================================================#
+        if plot_pdf == 'yes':
+            #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 + '_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)  
+                ax.patch.set_alpha(ALPHA_figs)
+                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']
+                fig = neardeltadf.plot(kind='bar', color=plotcolours15, ylim=(ymin,ymax), ax=ax)
+                ax.xaxis.grid(False)
+                ax.patch.set_alpha(ALPHA_figs)
+                #fig.patch.set_alpha(ALPHA_figs)
+                #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)
+                fig.patch.set_alpha(ALPHA_figs)
+                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]
+                fig = test_sorted.plot(legend=False, color =  plotcolours36)
+                pdf.savefig(bbox_inches='tight', pad_inches=0.4)
+                fig.patch.set_alpha(ALPHA_figs)
+                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()

Analysis/Code/climdata_fcts.py

@@ -245,6 +245,60 @@ def import_present_day_climdata_csv(Estuary, Clim_var_type):
     return Present_day_df, minplotDelta, maxplotDelta
 
 
+def quant_quant_scaling(Full_df, quantiles, Plotbool):
+    """
+    calculates the % "deltas" for each quantile in line with Climate Change in Australia recommendations provided here:
+    https://www.climatechangeinaustralia.gov.au/en/support-and-guidance/using-climate-projections/application-ready-data/scaling-methods/    
+    """
+    Periods = ['1990', '2020', '2060']
+    quantiles_srings = [str(x) for x in quantiles][:-1]
+    models = list(Full_df.resample('A').max().mean().index)
+    newmodel = []
+    for each in models:
+        newmodel.append(each[:-5])
+    unique_models = list(set(newmodel))
+    #create empty df and loop through models and periods to derive the % change factors for all quantiles and models
+    Quant_diff_df_outer = pd.DataFrame()
+    for unique_model in unique_models:
+        dfdiff = Full_df.filter(regex= unique_model)
+        Quant_diff_df = pd.DataFrame() 
+        for period in Periods:
+            x=dfdiff.filter(regex= period).dropna().values
+            x.sort(axis=0)
+            df=pd.DataFrame(x)
+            df.columns = ['A']
+            cut_df = pd.DataFrame(df.groupby(pd.qcut(df.rank(method='first').A, quantiles))['A'].mean().values)
+            cut_df.columns = [period]
+            Quant_diff_df = pd.concat([Quant_diff_df, cut_df], axis=1, join='outer')  
+        if Plotbool:
+            Quant_diff_df.plot(x=Quant_diff_df.index, y=Periods, kind="bar", title = unique_model)
+        Quant_diff_df['NF_%change_'+ unique_model] =  (Quant_diff_df['2020'] - Quant_diff_df['1990'])/Quant_diff_df['1990']*100
+        Quant_diff_df['FF_%change_'+ unique_model] =  (Quant_diff_df['2060'] - Quant_diff_df['1990'])/Quant_diff_df['1990']*100
+        Quant_diff_df = Quant_diff_df.replace([np.inf, -np.inf], np.nan)
+        Quant_diff_df = Quant_diff_df.iloc[:,[3,4]].fillna(0)
+        Quant_diff_df_outer = pd.concat([Quant_diff_df_outer, Quant_diff_df], axis=1, join='outer')
+    Quant_diff_df_outer.index = quantiles_srings 
+    if Plotbool:
+        Quant_diff_df_outer.plot(x=Quant_diff_df_outer.index, y=Quant_diff_df_outer.columns, kind="bar", legend = False, title='% intensity change per quantile')
+    #add new cols and rows with summary statistics
+    Quantile_change_NF_df = Quant_diff_df_outer.filter(regex= 'NF')
+    Quantile_change_FF_df = Quant_diff_df_outer.filter(regex= 'FF')
+    Quantile_change_stats_df = pd.DataFrame()
+    for loop_df in [Quantile_change_NF_df, Quantile_change_FF_df]:
+        Sum95_df = pd.DataFrame(loop_df.iloc[-5:,:].sum()).transpose()
+        Sum95_df.index = ['Sum>95perc']
+        Sum_df = pd.DataFrame(loop_df.sum()).transpose()
+        Sum_df.index = ['Sum']
+        Med_df = pd.DataFrame(loop_df.median()).transpose()
+        Med_df.index = ['Median']
+        loop_df = loop_df.append(Sum_df)
+        loop_df = loop_df.append(Sum95_df)
+        loop_df = loop_df.append(Med_df)
+        loop_df['Median'] = loop_df.median(axis=1)
+        loop_df['Maximum'] = loop_df.max(axis=1)
+        loop_df['Minimum'] = loop_df.min(axis=1)
+        Quantile_change_stats_df = pd.concat([Quantile_change_stats_df, loop_df], axis=1, join='outer')
+    return Quantile_change_stats_df
         
 
 

Analysis/Code/ggplot_time_series_with_trends.R

@@ -15,6 +15,26 @@ options(scipen=999)
 setwd("C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdoc/Projects/Paper#1/")
 ######################
 
+
+
+GRACE.monthly.av.df <- data.frame(read.csv("./Data/GRACE/Zunyi/CSV/With_border_pixels/MDB.monthly.TWSA.Mean.2002-2017 with border pixels.csv"))
+GRACE.monthly.std.df <- data.frame(read.csv("./Data/GRACE/Zunyi/CSV/With_border_pixels/MDB.monthly.TWSA.Stdv.2002-2017 with border pixels.csv"))
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 ######################
 #Set inputs
 ######################
@@ -38,7 +58,7 @@ if (file.exists(Output.Directory)){
   dir.create(file.path(Output.Directory))
   print('output folder did not exist and was created')
   }
-Output.Directory <- paste('./Output/CASESTUDY2_V4/', Estuary,'/Recent_Trends/Riv_', Riv.Gauge.loc,'_Est_',Est.Gauge.loc,'_GAMk', ggplotGAM.k, 'b/', sep="")
+Output.Directory <- paste('./Output/CASESTUDY2_V4/', Estuary,'/Recent_Trends/Riv_', Riv.Gauge.loc,'_Est_',Est.Gauge.loc,'_GAMk', ggplotGAM.k, 'd/', sep="")
 if (file.exists(Output.Directory)){
   print('output folder already existed and was not created again')
 } else {
@@ -82,12 +102,12 @@ colnames(AirT.full.df) <- 'tasmean'
 AirT.df$Julday1 <- seq(1,length(AirT.df[,1]),1)
 linear.trend.model_EC_all <- lm(tasmean ~ Julday1, AirT.df)
 AirT.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4] 
-AirT.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 356
+AirT.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 365
 ############
 AirT.full.df$Julday1 <- seq(1,length(AirT.full.df[,1]),1)
 linear.trend.model_EC_all <- lm(tasmean ~ Julday1, AirT.full.df)
 AirT.full.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4] 
-AirT.full.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 356
+AirT.full.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 365
 ############tasmean
 
 
@@ -116,11 +136,11 @@ colnames(RivT.full.df) <- 'rivTmean'
 RivT.df$Julday1 <- seq(1,length(RivT.df[,1]),1)
 linear.trend.model_EC_all <- lm(rivTmean ~ Julday1, RivT.df)
 RivT.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4] 
-RivT.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 356
+RivT.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 365
 RivT.full.df$Julday1 <- seq(1,length(RivT.full.df[,1]),1)
 linear.trend.model_EC_all <- lm(rivTmean ~ Julday1, RivT.full.df)
 RivT.full.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4] 
-RivT.full.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 356
+RivT.full.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 365
 ############River temp
 
 #export interpolated data as csv
@@ -141,10 +161,9 @@ colnames(dat) <- gsub(x=colnames(dat), pattern="Discharge (ML/d)",replacement="F
 RivQ.df <- dat
 rm(dat,char.df)
 #RivQ.full.TS <- zoo(log10((as.numeric(as.character(RivQ.df$Flow))) +1) , order.by= as.Date(RivQ.df[,"Date"], format = "%d/%m/%Y")) #=daily time series of rainfall for creation of clean, daily TS of ET and Q
-RivQ.full.TS <- zoo(as.numeric(as.character(RivQ.df$Flow)), order.by= as.Date(RivQ.df[,"Date"], format = "%d/%m/%Y")) #=daily time series of rainfall for creation of clean, daily TS of ET and Q
-
-#RivQ.TS <- window(RivQ.full.TS, start=as.Date("1990-01-01"), end=as.Date("2018-01-01"))
-RivQ.full.TS <- window(RivQ.full.TS, start=as.Date("2013-06-01"), end=as.Date("2018-06-01"))
+RivQ.full.TS <- zoo((as.numeric(as.character(RivQ.df$Flow)))/86.4, order.by= as.Date(RivQ.df[,"Date"], format = "%d/%m/%Y")) #=daily time series of rainfall for creation of clean, daily TS of ET and Q
+RivQ.TS <- window(RivQ.full.TS, start=as.Date("1990-01-01"), end=as.Date("2018-01-01"))
+#RivQ.full.TS <- window(RivQ.full.TS, start=as.Date("2013-06-01"), end=as.Date("2018-06-01"))
 RivQ.full.df <- data.frame(RivQ.full.TS)   ### This is only done because 
 RivQ.df <- data.frame(RivQ.TS)
 colnames(RivQ.df) <- 'RivQmean'
@@ -153,11 +172,11 @@ colnames(RivQ.full.df) <- 'RivQmean'
 RivQ.df$Julday1 <- seq(1,length(RivQ.df[,1]),1)
 linear.trend.model_EC_all <- lm(RivQmean ~ Julday1, RivQ.df)
 RivQ.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4] 
-RivQ.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 356
+RivQ.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 365
 RivQ.full.df$Julday1 <- seq(1,length(RivQ.full.df[,1]),1)
 linear.trend.model_EC_all <- lm(RivQmean ~ Julday1, RivQ.full.df)
 RivQ.full.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4] 
-RivQ.full.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 356
+RivQ.full.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 365
 ############River Flow
 
 
@@ -165,24 +184,23 @@ RivQ.full.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 35
 #Load a daily (no gaps) time series as a time serie baseline for other time series used here
 SST.df <- data.frame(read.csv(SST_CSV_Path))
 SST.full.TS <- zoo(SST.df$NNRP_R1_1950-273.15, order.by= as.Date(SST.df[,"X"], format = "%Y-%m-%d")) #=daily time series of rainfall for creation of clean, daily TS of ET and Q
-#SST.TS <- window(SST.full.TS, start=as.Date("1990-01-01"), end=as.Date("2018-01-01"))
+SST.TS <- window(SST.full.TS, start=as.Date("1990-01-01"), end=as.Date("2018-01-01"))
 
-SST.full.TS <- window(SST.full.TS, start=as.Date("2013-06-01"), end=as.Date("2018-06-01"))
+#SST.full.TS <- window(SST.full.TS, start=as.Date("2013-06-01"), end=as.Date("2018-06-01"))
 
 SST.full.df <- data.frame(SST.full.TS)
 SST.df <- data.frame(SST.TS)
-str(SST.df)
 colnames(SST.df) <- 'SSTmean'
 colnames(SST.full.df) <- 'SSTmean'
 ############
 SST.full.df$Julday1 <- seq(1,length(SST.full.df[,1]),1)
 linear.trend.model_EC_all <- lm(SSTmean ~ Julday1, SST.full.df)
 SST.full.pvalNCV_ECall <- summary(linear.trend.model_EC_all)$coefficients[2,4] 
-SST.full.lintrend <- summary(linear.trend.model_EC_all)$coefficients[2,1] * 356
+SST.full.lintrend <- summary(linear.trend.model_EC_all)$coefficients[2,1] * 365
 SST.df$Julday1 <- seq(1,length(SST.df[,1]),1)
 linear.trend.model_EC_all2 <- lm(SSTmean ~ Julday1, SST.df)
 SST.pvalNCV_ECall <- summary(linear.trend.model_EC_all2)$coefficients[2,4] 
-SST.lintrend <- summary(linear.trend.model_EC_all2)$coefficients[2,1] * 356
+SST.lintrend <- summary(linear.trend.model_EC_all2)$coefficients[2,1] * 365
 ############ SST
 
 
@@ -214,11 +232,11 @@ colnames(EstT.full.df) <- 'EstTmean'
 EstT.df$Julday1 <- seq(1,length(EstT.df[,1]),1)
 linear.trend.model_EC_all <- lm(EstTmean ~ Julday1, EstT.df)
 EstT.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4] 
-EstT.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 356
+EstT.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 365
 EstT.full.df$Julday1 <- seq(1,length(EstT.full.df[,1]),1)
 linear.trend.model_EC_all <- lm(EstTmean ~ Julday1, EstT.full.df)
 EstT.full.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4] 
-EstT.full.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 356
+EstT.full.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 365
 ############Est temp
 
 
@@ -353,45 +371,45 @@ colnames(combined.df) <- c('tasmean','rivTmean', 'SSTmean', 'EstTmean', 'rivQmea
 
 #Air temp 
 p2air <- ggplot(combined.df, aes(y=tasmean, x=index(combined.TS))) + geom_line(alpha=0.5) + 
-  ggtitle(paste(Estuary," Catchment centroid - Linear and smooth trend in catchment airT (SILO) lin trend was ",
+  ggtitle(paste(Estuary," Catchment centroid - Linear and smooth trend in catchment airT (SILO) linear trend was ",
                 round(AirT.lintrend,3), ' C°/year with p=', round(AirT.pvalNCV_ECall,10), sep=" ")) +
   theme(plot.title=element_text(face="bold", size=9)) +
   geom_smooth(method='lm',fill="green", formula=y~x, colour="darkgreen", size = 0.5) +
   stat_smooth(method=gam, formula=y~s(x, k=ggplotGAM.k), se=T, size=0.5) +
-  ylab("Air Temperature [C°]") + xlab("Time")+
+  ylab("Air Temperature [C°]") + xlab(NULL)+
   theme(axis.text=element_text(size=Fontsize)) + 
   theme(panel.grid.major.x = element_blank() ,panel.grid.minor.x = element_blank(), panel.grid.major.y = element_line( size=.1, color="white" ))
 
 #Riv temp 
 p2riv <- ggplot(combined.df, aes(y=rivTmean, x=index(combined.TS))) + geom_line(alpha=0.5) + 
-  ggtitle(paste(Estuary,'@', Riv.Gauge.loc, " - Linear and smooth trend in river temperature (GAUGE) lin trend was ",
+  ggtitle(paste(Estuary,'@', Riv.Gauge.loc, " - Linear and smooth trend in river temperature (GAUGE) linear trend was ",
                 round(RivT.lintrend,3), ' C°/year with p=', round(RivT.pvalNCV_ECall,10), sep=" ")) +
   theme(plot.title=element_text(face="bold", size=9)) +
   geom_smooth(method='lm',fill="green", formula=y~x, colour="darkgreen", size = 0.5) +
   stat_smooth(method=gam, formula=y~s(x, k=ggplotGAM.k), se=T, size=0.5) +
-  ylab("River Temperature [C°]") + xlab("Time")+
+  ylab("River Temperature [C°]") + xlab(NULL) +  
   theme(axis.text=element_text(size=Fontsize)) + 
   theme(panel.grid.major.x = element_blank() ,panel.grid.minor.x = element_blank(), panel.grid.major.y = element_line( size=.1, color="white" ))
 
 #Riv flow
 if(logtransformFlow ==TRUE){
   p2rivQ <- ggplot(combined.df, aes(y=log10(rivQmean+2), x=index(combined.TS))) + geom_line(alpha=0.5) + 
-    ggtitle(paste(Estuary,'@', Riv.Gauge.loc, " - Linear and smooth trend in river flow (GAUGE) lin trend was ",
-                  round(RivQ.full.lintrend,3), ' ML/day /year with p=', round(RivQ.full.pvalNCV_ECall,10), sep=" ")) +
+    ggtitle(paste(Estuary,'@', Riv.Gauge.loc, " - Linear and smooth trend in river flow (GAUGE) linear trend was ",
+                  round(RivQ.full.lintrend,3),  'cubic-meter/year with p=', round(RivQ.full.pvalNCV_ECall,10), sep=" ")) +
     theme(plot.title=element_text(face="bold", size=9)) +
     geom_smooth(method='lm',fill="green", formula=y~x, colour="darkgreen", size = 0.5) +
     stat_smooth(method=gam, formula=y~s(x, k=ggplotGAM.k), se=T, size=0.5) +
-    ylab(expression(paste("log10(River Flow [ML/day] + 2)", sep=""))) + xlab("Time")+
+    ylab(expression(paste("log10(River flow [m"^"3" *"/s] + 2)", sep=""))) + xlab(NULL)+
     theme(axis.text=element_text(size=Fontsize)) + 
     theme(panel.grid.major.x = element_blank() ,panel.grid.minor.x = element_blank(), panel.grid.major.y = element_line( size=.1, color="white" ))
   } else {
   p2rivQ <- ggplot(combined.df, aes(y=rivQmean, x=index(combined.TS))) + geom_line(alpha=0.5) + 
-    ggtitle(paste(Estuary,'@', Riv.Gauge.loc, " - Linear and smooth trend in river flow (GAUGE) lin trend was ",
-                  round(RivT.lintrend,3), ' C°/year with p=', round(RivT.pvalNCV_ECall,10), sep=" ")) +
+    ggtitle(paste(Estuary,'@', Riv.Gauge.loc, " - Linear and smooth trend in river flow (GAUGE) linear trend was ",
+                  round(RivT.lintrend,3), 'cubic-meter/year with p=', round(RivT.pvalNCV_ECall,10), sep=" ")) +
     theme(plot.title=element_text(face="bold", size=9)) +
     geom_smooth(method='lm',fill="green", formula=y~x, colour="darkgreen", size = 0.5) +
     stat_smooth(method=gam, formula=y~s(x, k=ggplotGAM.k), se=T, size=0.5) +
-    ylab(expression(paste("River Flow [ML/day]", sep=""))) + xlab("Time")+
+    ylab(expression(paste("River flow [m"^"3" *"/s]", sep=""))) + xlab(NULL)+
     theme(axis.text=element_text(size=Fontsize)) + 
     theme(panel.grid.major.x = element_blank() ,panel.grid.minor.x = element_blank(), panel.grid.major.y = element_line( size=.1, color="white" ))
     }
@@ -409,7 +427,7 @@ p2sst <- ggplot(combined.df, aes(y=SSTmean, x=index(combined.TS))) + geom_line(a
 
 #sst temp 
 p2Est <- ggplot(combined.df, aes(y=EstTmean, x=index(combined.TS))) + geom_line(alpha=0.5) + 
-  ggtitle(paste(Estuary,'@', Est.Gauge.loc, " - Linear and smooth trend in Estuary temperature (GAUGE) lin trend was ",
+  ggtitle(paste(Estuary,'@', Est.Gauge.loc, " - Linear and smooth trend in Estuary temperature (GAUGE) linear trend was ",
                 round(EstT.lintrend,3), ' C°/year with p=', round(EstT.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) +
@@ -457,6 +475,16 @@ png(file = png.name, width = 10.5, height = 7, units='in', res=500)
 grid.arrange(gB,gC, ncol=1)
 dev.off()
 
+png.name <- paste(Output.Directory, Estuary, '@', Riv.Gauge.loc, '_Trends_AirT_RivT_1990-present_', Sys.Date(),".png", sep="")
+png(file = png.name, width = 10.5, height = 7, units='in', res=500)
+grid.arrange(gA,gB, ncol=1)
+dev.off()
+
+png.name <- paste(Output.Directory, Estuary, '@', Riv.Gauge.loc, '_Trends_RivQ_RivT_1990-present_', Sys.Date(),".png", sep="")
+png(file = png.name, width = 10.5, height = 7, units='in', res=500)
+grid.arrange(gB ,gE, ncol=1)
+dev.off()
+
 png.name <- paste(Output.Directory, Estuary, '@', Riv.Gauge.loc, '_Trends_SST_RivT_EstT_1990-present_', Sys.Date(),".png", sep="")
 png(file = png.name, width = 10.5, height = 7, units='in', res=500)
 grid.arrange(gB,gD,gC,ncol=1)
@@ -467,9 +495,9 @@ png(file = png.name, width = 10.5, height = 7, units='in', res=500)
 grid.arrange(gE,gB,gD,ncol=1)
 dev.off()
 
-png.name <- paste(Output.Directory, Estuary, '@', Riv.Gauge.loc, '_Trends_AirT_RivT_RivQ_1990-present_', Sys.Date(),".png", sep="")
-png(file = png.name, width = 10.5, height = 7, units='in', res=500)
-grid.arrange(gA,gB,gE,ncol=1)
+png.name <- paste(Output.Directory, Estuary, '@', Riv.Gauge.loc, '_Trends_AirT_RivT_RivQ_1990-present_', Sys.Date(),"b.png", sep="")
+png(file = png.name, width = 10.5, height = 11, units='in', res=500)
+grid.arrange(gB, gA,gE,ncol=1)
 dev.off()
 
 png.name <- paste(Output.Directory, Estuary, '@', Riv.Gauge.loc, '_Trends_AirT_RivT_1990-present_', Sys.Date(),".png", sep="")
@@ -482,8 +510,45 @@ png(file = png.name, width = 10.5, height = 7, units='in', res=500)
 grid.arrange(gB,gE,ncol=1)
 dev.off()
 
+lm_eqn <- function(df){
+  m <- lm(y ~ x, df);
+  eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2, 
+                   list(a = format(coef(m)[1], digits = 2), 
+                        b = format(coef(m)[2], digits = 2), 
+                        r2 = format(summary(m)$r.squared, digits = 3)))
+  as.character(as.expression(eq));  
+}
+
+Modeling <- 'yes'
+if (Modeling == 'yes'){
+  source("C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdoc/Projects/Paper#1/Analysis/Code/FUN_do_ggplot_Scatterplot_2vars.R")
+  AirT.full.TS.Monthly <- aggregate(AirT.full.TS, as.yearmon(time(AirT.full.TS)),FUN = mean)
+  RivT.full.TS.Monthly <- aggregate(RivT.full.TS, as.yearmon(time(RivT.full.TS)),FUN = mean)
+  RivQ.full.TS.Monthly <- aggregate(RivQ.full.TS, as.yearmon(time(RivQ.full.TS)),FUN = mean)
+  
+  if(logtransformFlow ==TRUE){
+    RivQ.full.TS.Monthly <- log10(RivQ.full.TS.Monthly +2) 
+  }
+  All.TS.zoo <- merge(AirT.full.TS.Monthly, RivT.full.TS.Monthly ,RivQ.full.TS.Monthly , all=F)
+  AirvsRivT <- Generate.ggplot.scatterplot(All.TS.zoo[,2], All.TS.zoo[,1],"River temp. [°C]" , "Air temp. [°C]")
+  if(logtransformFlow ==TRUE){
+    QvsRivT <- Generate.ggplot.scatterplot(All.TS.zoo[,2], All.TS.zoo[,3], "River temp. [°C]",  expression(paste("log10(River flow [m"^"3" *"/s] + 2)", sep="")))
+  } else {
+    QvsRivT <- Generate.ggplot.scatterplot(All.TS.zoo[,2], All.TS.zoo[,3], "River temp. [°C]",  expression(paste("River flow [m"^"3" *"/s]", sep="")))
+    
+  }
+  #export an overview plot as png
+  gA <- ggplotGrob(AirvsRivT)
+  gB <- ggplotGrob(QvsRivT)
+  gB$widths <- gA$widths
 
   
+  png.name <- paste(Output.Directory, Estuary, '@', Riv.Gauge.loc, '_Scatterplots_of_RivT_vs_predictors_1990-present_', Sys.Date(),"f.png", sep="")
+  png(file = png.name, width = 8, height = 4, units='in', res=500)
+  grid.arrange(gA,gB , ncol=2)
+  dev.off()
+  
+}