#several major imporvements on all fronts

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"))

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,232 +19,340 @@ 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':
+#==========================================================#
+     
+#loop through case study estuaries    
+Estuaries = ['HUNTER', 'RICHMOND', 'NADGEE', 'SHOALHAVEN', 'GEORGES','CATHIE']
+Estuaries = ['HUNTER']
+for Est in Estuaries:
+    
+    #==========================================================#
+    #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
+    
+    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)
+    #==========================================================#
+    
+    
+    
+    #==========================================================#
+    #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)
-    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 = Present_day_df.resample('Q-NOV').max() #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 )
-    
-    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
+    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:
-        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)
-    
-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)
+        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:
-        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)
-    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_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_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)
         
-# use transparent=True if you want the whole figure with a transparent background
+        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
         
         
     

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)
+    
+
+    
+    
+    
+    
+    
+    
+    
+    

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:

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('-------------------------------------------')
-#==========================================================#    
+Estuaries = ['HUNTER', 'RICHMOND', 'NADGEE', 'SHOALHAVEN', 'GEORGES','CATHIE']
+Estuaries = ['HUNTER']
 
-#==========================================================#
-Estuary_Folder = glob.glob('./Data/NARCLIM_Site_CSVs/'+ Case_Study_Name + '/' + Estuary + '*' )
-Clim_Var_CSVs = glob.glob(Estuary_Folder[0] + '/' + Clim_var_type + '*')
-#==========================================================#
+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'
+        #==========================================================#
         
-#==========================================================#
-#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
-#==========================================================#
+        #==========================================================#
+        #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('-------------------------------------------')
+        #==========================================================#    
         
-#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()
-#==========================================================#
+        #==========================================================#
+        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
+        #==========================================================#
         
-#==========================================================#
-#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)
-#==========================================================#
+        #==========================================================#
+        #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
+        #==========================================================#
         
-#==========================================================#
-#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)
         
+        #==========================================================#
+        #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()
+        #==========================================================#
         
-#==========================================================#
-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']
-#==========================================================#
+        #==========================================================#
+        #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)
         
-#==========================================================#
-#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()
+        
+        #==========================================================#
+        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()

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)
 #==========================================================# 

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
+        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
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

Analysis/Code/ggplot_time_series_with_trends_Waves_SST_OEH.R

@@ -0,0 +1,165 @@
+#R code for creating ggplots of time series with smooth (GAM) and linear term for the OEH buoy data on wave heigth and SST
+
+
+######################
+#Import Libraries and set working directory
+######################
+library(zoo)
+library(hydroTSM) #you need to install these packages first before you can load them here
+library(lubridate)
+library(mgcv)
+library(ggplot2)
+library(gridExtra)
+library(scales)
+options(scipen=999)
+setwd("C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdoc/Projects/Paper#1/")
+######################
+
+######################
+#Set inputs
+######################
+Case.Study <- "CASESTUDY2"
+Estuary <- "RICHMOND"
+Climvar <- 'tasmean'
+ggplotGAM.k <- 7
+######################
+
+
+######################
+Output.Directory <- paste('./Output/', Case.Study, '/', Estuary,'/Recent_Trends/', sep="")
+if (file.exists(Output.Directory)){
+  print('output folder already existed and was not created again')
+} else {
+  dir.create(file.path(Output.Directory))
+  print('output folder did not exist and was created')
+}
+######################
+
+
+######################
+#Set input file paths
+######################
+pattern = paste(Estuary, '.*.csv', sep="")
+Buoy_CSV_Path <- list.files("./Data/Ocean_Data/Waves_SST/", pattern, full.names=T)
+Buoy.name <- substr(list.files("./Data/Ocean_Data/Waves_SST/", pattern, full.names=F), nchar(Estuary)+2, nchar(Estuary)+7)
+
+SST.df <- data.frame(read.csv(AirT_CSV_Path, colClasses = "character"))
+colnames(SST.df) <- as.character(SST.df[9,])
+SST.df = SST.df[-c(1:9), ] 
+
+colnames(SST.df) <- gsub(x=colnames(SST.df), pattern="Date/Time",replacement="Date",fixed=T)
+colnames(SST.df) <- gsub(x=colnames(SST.df), pattern=" ",replacement="",fixed=T)
+colnames(SST.df) <- gsub(x=colnames(SST.df), pattern="Hsig(m)",replacement="Hsig",fixed=T)
+colnames(SST.df) <- gsub(x=colnames(SST.df), pattern="SeaTemp(C)",replacement="SST",fixed=T)
+SST.df$Date <- substr(SST.df$Date, 1, 10)
+
+
+######################
+#Hsig
+######################
+Hsig.full.TS <- zoo(as.numeric(SST.df$Hsig), order.by= as.Date(SST.df$Date, format = "%d/%m/%Y")) #=daily time series of rainfall for creation of clean, daily TS of ET and Q
+Hsig.full.TS <- aggregate(Hsig.full.TS, index(Hsig.full.TS) , FUN=mean)
+Hsig.full.TS <-  aggregate(Hsig.full.TS , as.yearmon(time(Hsig.full.TS )),FUN = mean)
+Hsig.full.TS <- Hsig.full.TS [1:length(Hsig.full.TS)-1]
+Hsig.TS <- window(Hsig.full.TS, start=as.Date("1990-01-01"), end=as.Date("2018-01-01"))
+Hsig.full.df <- data.frame(Hsig.full.TS)
+Hsig.df <- data.frame(Hsig.TS)
+str(Hsig.df)
+colnames(Hsig.df) <- 'Hsigmean'
+colnames(Hsig.full.df) <- 'Hsigmean'
+#rid of outliers
+#Hsig.full.df$MSL <- replace(Hsig.full.df$MSL, which(Hsig.full.df$MSL > 4), NA)
+Hsig.full.df$Julday1 <- seq(1,length(Hsig.full.df[,1]),1)
+linear.trend.model_EC_all <- lm(Hsigmean ~ Julday1, Hsig.full.df)
+Hsig.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4] 
+Hsig.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 12
+######################
+
+######################
+#SST
+######################
+SST.full.TS <- zoo(as.numeric(SST.df$SST), order.by= as.Date(SST.df$Date, format = "%d/%m/%Y")) #=daily time series of rainfall for creation of clean, daily TS of ET and Q
+SST.full.TS  <- na.omit(SST.full.TS)
+SST.full.TS <- aggregate(SST.full.TS, index(SST.full.TS) , FUN=mean)
+SST.full.TS <- SST.full.TS[1:length(SST.full.TS)-1]
+SST.TS <- window(SST.full.TS, start=as.Date("1990-01-01"), end=as.Date("2018-01-01"))
+SST.full.df <- data.frame(SST.full.TS)
+SST.df <- data.frame(SST.TS)
+str(SST.df)
+tail(SST.full.TS)
+colnames(SST.df) <- 'SSTmean'
+colnames(SST.full.df) <- 'SSTmean'
+#rid of outliers
+#SST.full.df$MSL <- replace(SST.full.df$MSL, which(SST.full.df$MSL > 4), NA)
+SST.full.df$Julday1 <- seq(1,length(SST.full.df[,1]),1)
+linear.trend.model_EC_all <- lm(SSTmean ~ Julday1, SST.full.df)
+SST.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4] 
+SST.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 365
+######################
+
+
+######################
+#Plot
+######################
+
+##################################### Full Time Period
+p1air <- ggplot(Hsig.full.df, aes(y=Hsigmean, x=index(Hsig.full.TS))) + geom_line(alpha=0.5) + 
+  ggtitle(paste(Estuary, " - Linear and smooth trend in significant wave height (OEH Buoy:",Buoy.name,  "| lin trend was ",
+                round(Hsig.lintrend,3), ' m/year with p=', round(Hsig.pvalNCV_ECall,10), sep=" ")) +
+  theme(plot.title=element_text(face="bold", size=9)) +
+  geom_smooth(method='lm',fill="green", formula=y~x, colour="darkgreen", size = 0.5) +
+  stat_smooth(method=gam, formula=y~s(x,  k=4), se=T, size=0.5, col="red") +
+  ylab("Significant Wave Height [m]") + xlab("Time")
+
+#export to png
+png.name <- paste(Output.Directory, '/Trends_Significant_Wave_Height_',Buoy.name,  '_full_period_', Sys.Date(),".png", sep="")
+png(file = png.name, width = 10.5, height = 7, units='in', res=500)
+grid.arrange(p1air,ncol=1)
+dev.off()
+
+#multiple smooths
+p1air <- ggplot(Hsig.full.df, aes(y=Hsigmean, x=index(Hsig.full.TS))) + geom_line(alpha=0.5) + 
+  ggtitle(paste(Estuary, " - Linear and smooth trend in monthly mean sea level (OEH Buoy:",Buoy.name,  "| lin trend was ",
+                round(Hsig.lintrend,3), ' m/year with p=', round(Hsig.pvalNCV_ECall,10), sep=" ")) +
+  theme(plot.title=element_text(face="bold", size=9)) +
+  stat_smooth(method=gam, formula=y~s(x,  k=13), se=T, size=0.5, col="red") +
+  #stat_smooth(method=gam, formula=y~s(x,  k=8), se=T, size=0.5, cor="blue") +
+  stat_smooth(method=gam, formula=y~s(x,  k=5), se=T, size=0.5, col="green") +
+  ylab("Significant Wave Height [m]") + xlab("Time")
+
+#export to png
+png.name <- paste(Output.Directory, '/Trends_Significant_Wave_Height_',Buoy.name,  '_MultiGAM_full_period_', Sys.Date(),".png", sep="")
+png(file = png.name, width = 10.5, height = 7, units='in', res=500)
+grid.arrange(p1air,ncol=1)
+dev.off()
+
+##################################### Full Time Period SST
+p1air <- ggplot(SST.full.df, aes(y=SSTmean, x=index(SST.full.TS))) + geom_line(alpha=0.5) + 
+  ggtitle(paste(Estuary, " - Linear and smooth trend in SST (OEH Buoy:",Buoy.name,  "| lin trend was ",
+                round(SST.lintrend,3), ' deg C/year with p=', round(SST.pvalNCV_ECall,10), sep=" ")) +
+  theme(plot.title=element_text(face="bold", size=9)) +
+  geom_smooth(method='lm',fill="green", formula=y~x, colour="darkgreen", size = 0.5) +
+  stat_smooth(method=gam, formula=y~s(x,  k=4), se=T, size=0.5, col="red") +
+  ylab("Sea Surface Temperature [C°]") + xlab("Time")
+
+#export to png
+png.name <- paste(Output.Directory, '/Trends_SST_',Buoy.name,  '_full_period_',  Sys.Date(),".png", sep="")
+png(file = png.name, width = 10.5, height = 7, units='in', res=500)
+grid.arrange(p1air,ncol=1)
+dev.off()
+
+#multiple smooths
+p1air <- ggplot(SST.full.df, aes(y=SSTmean, x=index(SST.full.TS))) + geom_line(alpha=0.5) + 
+  ggtitle(paste(Estuary, " - Linear and smooth trend in monthly mean sea level (OEH Buoy:",Buoy.name,  "| lin trend was ",
+                round(SST.lintrend,3), ' deg C/year with p=', round(SST.pvalNCV_ECall,10), sep=" ")) +
+  theme(plot.title=element_text(face="bold", size=9)) +
+  stat_smooth(method=gam, formula=y~s(x,  k=13), se=T, size=0.5, col="red") +
+  #stat_smooth(method=gam, formula=y~s(x,  k=8), se=T, size=0.5, cor="blue") +
+  stat_smooth(method=gam, formula=y~s(x,  k=5), se=T, size=0.5, col="green") +
+  ylab("Sea Surface Temperature [C°]") + xlab("Time")
+
+#export to png
+png.name <- paste(Output.Directory, '/Trends_SST_',Buoy.name, '_MultiGAM_full_period_', Sys.Date(),".png", sep="")
+png(file = png.name, width = 10.5, height = 7, units='in', res=500)
+grid.arrange(p1air,ncol=1)
+dev.off()

Analysis/Code/ggplot_time_series_with_trends_sea_level.R

@@ -19,41 +19,64 @@ setwd("C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdoc/P
 #Set inputs
 ######################
 Case.Study <- "CASESTUDY2"
-Estuary <- "HUNTER"
+Estuary <- "NADGEE"
 Climvar <- 'tasmean'
 ggplotGAM.k <- 7
 ######################
 
 
 ######################
-#Set input file paths
+Output.Directory <- paste('./Output/', Case.Study, '/', Estuary,'/Recent_Trends/', sep="")
+if (file.exists(Output.Directory)){
+  print('output folder already existed and was not created again')
+} else {
+  dir.create(file.path(Output.Directory))
+  print('output folder did not exist and was created')
+}
 ######################
-AirT_CSV_Path <- "./Data/Ocean_Data/BOM_monthly_SL_Hunter_Newcastle.txt"
-
-
 
 
+######################
+#Set input file paths
+######################
+AirT_CSV_Path <- paste("./Data/Ocean_Data/BOM_monthly_SL_",Estuary,  "_Eden.txt", sep="")
 
 
 dat = readLines(AirT_CSV_Path)
 dat = as.data.frame(do.call(rbind, strsplit(dat, split=" {2,10}")), stringsAsFactors=FALSE)
 colnames(dat) <-dat[3,]
 
-dat2 = dat[-c(1:3), ] 
-dat2 = dat2[-(720:726),]
-dat2 = dat2[-(1:108),]
-dat2 = dat2[,-1]
+if (Estuary=='NADGEE'){
+  dat2 = dat[-c(1:11), ] 
+  dat2 = dat2[-(365:381),]
+  dat2 = dat2[,-1]
+  }
+
+if (Estuary=='HUNTER'){
+  dat2 = dat[-c(1:3), ]  
+  dat2 = dat2[-(720:726),]
+  dat2 = dat2[-(1:108),]
+  dat2 = dat2[,-1]
+  }
+
+
 
 dat2$Date <-  as.yearmon(dat2[,1], "%m %Y")
 
 SeaLev.df <- dat2
 head(SeaLev.df)
 SeaLev.df$MSL <- as.numeric(SeaLev.df$Mean)
+
+#rid of outliers
+SeaLev.df$MSL <- replace(SeaLev.df$MSL, which(SeaLev.df$MSL > 4), NA)
+  
 SeaLev.df$Julday1 <- seq(1,length(SeaLev.df[,1]),1)
 linear.trend.model_EC_all <- lm(MSL ~ Julday1, SeaLev.df)
 SeaLev.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4] 
 SeaLev.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 12
 
+
+
 ######################
 #Plot
 ######################
@@ -68,7 +91,7 @@ p1air <- ggplot(SeaLev.df, aes(y=MSL, x=Date)) + geom_line(alpha=0.5) +
   ylab("Monthly Mean Sea Level [m]") + xlab("Time")
 
 #export to png
-png.name <- paste('./Output/', Case.Study, '/', Estuary, '/Trends_MonthlyMeanSeaLevel_full_period_', Sys.Date(),".png", sep="")
+png.name <- paste(Output.Directory, '/Trends_MonthlyMeanSeaLevel_full_period_', Sys.Date(),".png", sep="")
 png(file = png.name, width = 10.5, height = 7, units='in', res=500)
 grid.arrange(p1air,ncol=1)
 dev.off()
@@ -84,7 +107,7 @@ p1air <- ggplot(SeaLev.df, aes(y=MSL, x=Date)) + geom_line(alpha=0.5) +
   ylab("Monthly Mean Sea Level [m]") + xlab("Time")
 
 #export to png
-png.name <- paste('./Output/', Case.Study, '/', Estuary, '/Trends_MonthlyMeanSeaLevel_MultiGAM_full_period_', Sys.Date(),".png", sep="")
+png.name <- paste(Output.Directory, '/Trends_MonthlyMeanSeaLevel_MultiGAM_full_period_', Sys.Date(),".png", sep="")
 png(file = png.name, width = 10.5, height = 7, units='in', res=500)
 grid.arrange(p1air,ncol=1)
 dev.off()