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

Compare commits

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

2 Commits
c5f2f4f743 ... 9eda19058e

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

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

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

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

@ -1,13 +1,14 @@
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
#####################################---------------------------------- #==========================================================#
#Last Updated - March 2018 #Last Updated - March 2018
#@author: z5025317 Valentin Heimhuber #@author: z5025317 Valentin Heimhuber
#code for creating future climate variability deviation plots for NARCLIM variables. #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 #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 #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 #Load packages
#####################################---------------------------------- #==========================================================#
import numpy as np import numpy as np
import os import os
import pandas as pd import pandas as pd
@ -18,234 +19,342 @@ from datetime import datetime
from datetime import timedelta from datetime import timedelta
from matplotlib.backends.backend_pdf import PdfPages from matplotlib.backends.backend_pdf import PdfPages
from ggplot import * from ggplot import *
import csv
matplotlib.style.use('ggplot') 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) # 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/') os.chdir('C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdoc/Projects/Paper#1/')
# #==========================================================#
#####################################----------------------------------
#set input parameters #loop through case study estuaries
Base_period_start = '1986-01-01' Estuaries = ['HUNTER', 'RICHMOND', 'NADGEE', 'SHOALHAVEN', 'GEORGES','CATHIE']
Base_period_end = '2005-01-01' #use last day that's not included in period as < is used for subsetting Estuaries = ['HUNTER']
Estuary = 'Nadgee' # 'Belongil' for Est in Estuaries:
Clim_var_type = "*" # '*' will create pdf for all variables in folder
Clim_var_type = "pracc*" # '*' will create pdf for all variables in folder
Present_Day_Clim_Var = 'Rainfall' #MaxT, MinT, Rainfall, ET Wind
present_day_plot = 'yes'
Version = "V2"
Stats = 'dailymax' #maximum takes the annual max Precipitation instead of the sum
#####################################----------------------------------
#set directory path for output files
output_directory = 'Output/Clim_Deviation_Plots/'+ Estuary
#output_directory = 'J:/Project wrl2016032/NARCLIM_Raw_Data/Extracted'
if not os.path.exists(output_directory):
os.makedirs(output_directory)
print('-------------------------------------------')
print("output directory folder didn't exist and was generated")
print('-------------------------------------------')
print('-------------------')
Clim_Var_CSVs = glob.glob('./Output/' + Estuary + '/' + Estuary + '_' + Clim_var_type[:-1] + '_' + Stats + '*')
#read CSV file
clim_var_csv_path = Clim_Var_CSVs[0]
Filename = os.path.basename(os.path.normpath(clim_var_csv_path))
Clim_var_type = Filename.split('_', 2)[1]
print(Clim_var_type)
Ensemble_Delta_full_df = pd.read_csv(clim_var_csv_path, index_col=0, parse_dates = True)
#Ensemble_Delta_full_df = pd.to_numeric(Ensemble_Delta_full_df)
#
#load present day climate data for the same variable
if Clim_var_type == 'evspsblmean': #ET time series that we have is not in the same format as the other variables, hence the different treatment
Present_day_Var_CSV = glob.glob('./Data/Wheather_Station_Data/**/' + Estuary + '_' + Present_Day_Clim_Var + '*csv')
Present_day_df = pd.read_csv(Present_day_Var_CSV[0])
Dates = pd.to_datetime(Present_day_df.Date)
Present_day_df.index = Dates
Present_day_df = Present_day_df.iloc[:,1]
Present_day_df = Present_day_df.replace(r'\s+', np.nan, regex=True)
Present_day_df = pd.to_numeric(Present_day_df)
Present_day_df.index = Dates
[minplotDelta, maxplotDelta]=[50,50]
#for tasmean, observed min and max T need to be converted into mean T
elif Clim_var_type == 'tasmean':
Present_day_Var_CSV = glob.glob('./Data/Wheather_Station_Data/**/' + Estuary + '_MaxT*csv')
Present_day_df = pd.read_csv(Present_day_Var_CSV[0])
Dates = pd.to_datetime(Present_day_df.Year*10000+Present_day_df.Month*100+Present_day_df.Day,format='%Y%m%d')
Present_day_df.index = Dates
Present_day_MaxT_df = Present_day_df.iloc[:,5]
Present_day_Var_CSV = glob.glob('./Data/Wheather_Station_Data/**/' + Estuary + '_MinT*csv')
Present_day_df = pd.read_csv(Present_day_Var_CSV[0])
Dates = pd.to_datetime(Present_day_df.Year*10000+Present_day_df.Month*100+Present_day_df.Day,format='%Y%m%d')
Present_day_df.index = Dates
Present_day_MinT_df = Present_day_df.iloc[:,5]
Present_day_df = (Present_day_MaxT_df + Present_day_MinT_df)/2
[minplotDelta, maxplotDelta]=[1,2]
elif Clim_var_type == 'tasmax':
[minplotDelta, maxplotDelta]=[1,2]
elif Clim_var_type == 'wssmean' or Clim_var_type == 'wss1Hmaxtstep':
Present_day_Var_CSV = glob.glob('./Data/Wheather_Station_Data/**/' + Estuary + '_' + Present_Day_Clim_Var + '*csv')
Present_day_df = pd.read_csv(Present_day_Var_CSV[0])
Present_day_df.index = Present_day_df[['Year', 'Month', 'Day', 'Hour']].apply(lambda s : datetime(*s),axis = 1)
Present_day_df = Present_day_df.filter(regex= 'm/s')
Present_day_df = Present_day_df.replace(r'\s+', np.nan, regex=True)
Present_day_df['Wind speed measured in m/s'] = Present_day_df['Wind speed measured in m/s'].convert_objects(convert_numeric=True)
[minplotDelta, maxplotDelta]=[1, 1.5]
else:
Present_day_Var_CSV = glob.glob('./Data/Wheather_Station_Data/**/' + Estuary + '_' + Present_Day_Clim_Var + '*csv')
Present_day_df = pd.read_csv(Present_day_Var_CSV[0])
Dates = pd.to_datetime(Present_day_df.Year*10000+Present_day_df.Month*100+Present_day_df.Day,format='%Y%m%d')
Present_day_df.index = Dates
Present_day_df = Present_day_df.iloc[:,5]
[minplotDelta, maxplotDelta]=[50,50]
#create seasonal sums etc.
if (Clim_var_type == 'pracc' or Clim_var_type == 'evspsblmean' or Clim_var_type == 'potevpmean'
or Clim_var_type == 'pr1Hmaxtstep' or Clim_var_type == 'wss1Hmaxtstep'):
if Stats == 'dailymax':
Present_day_df_annual = Present_day_df.resample('A').max()
Present_day_df_annual = Present_day_df_annual.replace(0, np.nan)
else:
Present_day_df_annual = Present_day_df.resample('A').sum()
Present_day_df_annual = Present_day_df_annual.replace(0, np.nan)
Fdf_Seas_means = Present_day_df.resample('Q-NOV').sum() #seasonal means
Fdf_Seas_means = Fdf_Seas_means.replace(0, np.nan)
else:
Present_day_df_annual = Present_day_df.resample('A').mean()
Present_day_df_monthly = Present_day_df.resample('M').mean()
Present_day_df_weekly = Present_day_df.resample('W').mean()
Fdf_Seas_means = Present_day_df.resample('Q-NOV').mean() #seasonal means
#Loop through annual and seasons and create a deviation plot for each.
times = ['annual', 'DJF', 'MAM', 'JJA','SON']
fig = plt.figure(figsize=(14,8))
delta_all_df = pd.DataFrame()
i=1
for temp in times:
temp = 'annual'
#subset the ensemble dataframe for the period used:
if temp == 'annual':
Ensemble_Delta_df = Ensemble_Delta_full_df.iloc[:,range(0,2)]
#
Present_Day_ref_df = Present_day_df_annual
else:
Ensemble_Delta_df = Ensemble_Delta_full_df.filter(regex= temp)
Ensemble_Delta_df.columns = ['near', 'far']
if temp == 'DJF':
Mean_df = Fdf_Seas_means[Fdf_Seas_means.index.quarter==1]
Column_names = ['DJF_near', 'DJF_far']
if temp == 'MAM':
Mean_df = Fdf_Seas_means[Fdf_Seas_means.index.quarter==2]
Column_names = ['MAM_near', 'MAM_far']
if temp == 'JJA':
Mean_df = Fdf_Seas_means[Fdf_Seas_means.index.quarter==3]
Column_names = ['JJA_near', 'JJA_far']
if temp == 'SON':
Mean_df = Fdf_Seas_means[Fdf_Seas_means.index.quarter==4]
Present_Day_ref_df = Mean_df
print(Ensemble_Delta_df.columns.values)
#Subset to present day variability period
Present_Day_ref_df = pd.DataFrame(Present_Day_ref_df.loc[(Present_Day_ref_df.index >= Base_period_start) & (Present_Day_ref_df.index <= Base_period_end)])
Present_Day_Mean = np.percentile(Present_Day_ref_df, 50)
Present_Day_SD = np.std(Present_Day_ref_df)
#create data frame for floating stacked barplots
index=['-2std', '-1std', 'Med', '1std', '2std']
columns = ['present','near future', 'far future']
Plot_in_df = pd.DataFrame(index=index, columns =columns)
#
Plot_in_df['present'] = [float(Present_Day_Mean-2*Present_Day_SD),float(Present_Day_Mean-Present_Day_SD), float(Present_Day_Mean),
float(Present_Day_Mean+Present_Day_SD), float(Present_Day_Mean+2*Present_Day_SD)]
Plot_in_df['near future'] = [float(Present_Day_Mean + Ensemble_Delta_df.near[-3:][0]),np.NaN, float(Present_Day_Mean + Ensemble_Delta_df.near[-3:][1]),
np.NaN, float(Present_Day_Mean + Ensemble_Delta_df.near[-3:][2])]
Plot_in_df['far future'] = [float(Present_Day_Mean + Ensemble_Delta_df.far[-3:][0]),np.NaN, float(Present_Day_Mean + Ensemble_Delta_df.far[-3:][1]),
np.NaN, float(Present_Day_Mean + Ensemble_Delta_df.far[-3:][2])]
#Create a second data frame that has the values in a way that they can be stacked up in bars with the correct absolute values
Plot_in_df2 = pd.DataFrame(index=index, columns =columns )
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)] #set input parameters
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])] Case_Study_Name = 'CASESTUDY2'
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]), Casestudy2_csv_path = "Data/NARCLIM_Site_CSVs/CASESTUDY2/CASESTDUY2_NARCLIM_Point_Sites.csv"
np.NaN, float(Ensemble_Delta_df.far[-3:][2]-Ensemble_Delta_df.far[-3:][1])] #Estuary = 'HUNTER' # 'Belongil'
#transpose the data frame Estuary = Est # 'Belongil'
Plot_in_df_tp = Plot_in_df2.transpose() print Estuary
#do the individual plots Base_period_start = '1970-01-01' #Start of interval for base period of climate variability
if temp == 'annual': Base_period_DELTA_start = '1990-01-01' #Start of interval used for calculating mean and SD for base period (Delta base period)
xmin = int(min(Plot_in_df.min(axis=1))-minplotDelta) Base_period_end = '2009-01-01' #use last day that's not included in period as < is used for subsetting
xmax = int(max(Plot_in_df.max(axis=1))+maxplotDelta) Clim_var_type = 'tasmax' #Name of climate variable in NARCLIM models '*' will create pdf for all variables in folder
else:
xmin = int(min(Plot_in_df.min(axis=1))-minplotDelta) PD_Datasource = 'SILO' #Source for present day climate data (historic time series) can be either: 'Station' or 'SILO'
xmax = int(max(Plot_in_df.max(axis=1))+maxplotDelta) SILO_Clim_var = ['max_temp'] #select the SILO clim variable to be used for base period. - see silo.py for detailed descriptions
#define colour scheme Location = 'Estuary' # pick locaiton for extracting the SILO data can be: Estuary, Catchment, or Ocean
#likert_colors = ['none', 'firebrick','firebrick','lightcoral','lightcoral']
likert_colors = ['none', 'darkblue', 'darkblue','cornflowerblue','cornflowerblue']
#plot the stacked barplot
ax=plt.subplot(2,3,i)
Plot_in_df_tp.plot.bar(stacked=True, color=likert_colors, edgecolor='none', legend=False, ax=ax)
z = plt.axhline(float(Present_Day_Mean-2*Present_Day_SD), linestyle='-', color='black', alpha=.5)
z.set_zorder(-1)
z = plt.axhline(float(Present_Day_Mean+2*Present_Day_SD), linestyle='-', color='black', alpha=.5)
z.set_zorder(-1)
z = plt.axhline(float(Present_Day_Mean-Present_Day_SD), linestyle='--', color='black', alpha=.5)
z.set_zorder(-1)
z = plt.axhline(float(Present_Day_Mean+Present_Day_SD), linestyle='--', color='black', alpha=.5)
z.set_zorder(-1)
z = plt.axhline(float(Present_Day_Mean), linestyle='--', color='red', alpha=.5)
z.set_zorder(-1)
plt.ylim(xmin, xmax)
plt.title(Clim_var_type + ' ' + temp )
ax.grid(False)
for tick in ax.get_xticklabels():
tick.set_rotation(0)
fig.tight_layout()
fig.patch.set_alpha(0)
#reset i to i+1 for next step
if temp == 'MAM':
i=i+2
else:
i=i+1
print(i)
plt.show()
out_file_name = Estuary + '_' + Clim_var_type + '_' + Stats + '_' + '_CC_prio_plot' + Version + '.png'
out_path = output_directory + '/' + out_file_name
fig.savefig(out_path)
if present_day_plot == 'yes': presentdaybar = False #include a bar for present day variability in the plots?
#print present day climate data present_day_plot = 'no' #save a time series of present day
fig = plt.figure(figsize=(5,4)) Version = "V1"
ax = fig.add_subplot(1, 1, 1) Stats = 'days_h_35' #'maxdaily' #maximum takes the annual max Precipitation instead of the sum
if temp == 'annual': ALPHA_figs = 1 #Set alpha of figure background (0 makes everything around the plot panel transparent)
xmin = int(min(Plot_in_df.min(axis=1))-minplotDelta) #==========================================================#
xmax = int(max(Plot_in_df.max(axis=1))+maxplotDelta)
else:
xmin = int(min(Plot_in_df.min(axis=1))-minplotDelta)
xmax = int(max(Plot_in_df.max(axis=1))+maxplotDelta)
Present_Day_ref_df.plot(legend=False, ax=ax)
z = plt.axhline(float(Present_Day_Mean-2*Present_Day_SD), linestyle='-', color='black', alpha=.5)
z.set_zorder(-1)
z = plt.axhline(float(Present_Day_Mean+2*Present_Day_SD), linestyle='-', color='black', alpha=.5)
z.set_zorder(-1)
z = plt.axhline(float(Present_Day_Mean-Present_Day_SD), linestyle='--', color='black', alpha=.5)
z.set_zorder(-1)
z = plt.axhline(float(Present_Day_Mean+Present_Day_SD), linestyle='--', color='black', alpha=.5)
z.set_zorder(-1)
z = plt.axhline(float(Present_Day_Mean), linestyle='--', color='red', alpha=.5)
z.set_zorder(-1)
#fig.patch.set_facecolor('deepskyblue')
fig.patch.set_alpha(0)
plt.ylim(13, xmax)
plt.show()
out_file_name = 'C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdoc/Projects/OEH_Coastal_Node_Deliverables/Technical_Report_1/Figures/tasmean_present_day_manual_backgroundTP.png'
out_path = output_directory + '/' + out_file_name
fig.savefig(out_file_name)
# use transparent=True if you want the whole figure with a transparent background #==========================================================#
#set directory path for output files
output_directory = 'Output/' + Case_Study_Name + '/' + Estuary + '/' + '/Clim_Deviation_Plots/'
#output_directory = 'J:/Project wrl2016032/NARCLIM_Raw_Data/Extracted'
if not os.path.exists(output_directory):
os.makedirs(output_directory)
print('-------------------------------------------')
print("output directory folder didn't exist and was generated")
print('-------------------------------------------')
print('-------------------')
#==========================================================#
#==========================================================#
#read CSV file
#==========================================================#
Clim_Var_CSVs = glob.glob('./Output/' + Case_Study_Name + '/' + Estuary + '/' + Estuary + '_' + Clim_var_type + '_' + Stats + '*')
clim_var_csv_path = Clim_Var_CSVs[0]
Filename = os.path.basename(os.path.normpath(clim_var_csv_path))
Clim_var_type = Filename.split('_', 2)[1]
print(Clim_var_type)
Ensemble_Delta_full_df = pd.read_csv(clim_var_csv_path, index_col=0, parse_dates = True)
if Clim_var_type == 'rsdsmean':
Clim_Var_CSVs = glob.glob('./Output/' + Case_Study_Name + '/' + Estuary + '/' + Estuary + '_rsdsmean_' + Stats + '*')
clim_var_csv_path = Clim_Var_CSVs[0]
Filename = os.path.basename(os.path.normpath(clim_var_csv_path))
Clim_var_type = Filename.split('_', 2)[1]
print(Clim_var_type)
Ensemble_Delta_full_df1 = pd.read_csv(clim_var_csv_path, index_col=0, parse_dates = True)
Clim_Var_CSVs = glob.glob('./Output/' + Case_Study_Name + '/' + Estuary + '/' + Estuary + '_rldsmean_' + Stats + '*')
clim_var_csv_path = Clim_Var_CSVs[0]
Filename = os.path.basename(os.path.normpath(clim_var_csv_path))
Clim_var_type = Filename.split('_', 2)[1]
print(Clim_var_type)
Ensemble_Delta_full_df2 = pd.read_csv(clim_var_csv_path, index_col=0, parse_dates = True)
Ensemble_Delta_full_df = Ensemble_Delta_full_df2 + Ensemble_Delta_full_df1
#Ensemble_Delta_full_df = pd.to_numeric(Ensemble_Delta_full_df)
#==========================================================#
#==========================================================#
#load present day climate variable time series
#==========================================================#
if PD_Datasource == 'Station':
Present_day_df,minplotDelta, maxplotDelta = fct.import_present_day_climdata_csv(Estuary, Clim_var_type)
if PD_Datasource == 'SILO':
#read the CSV to extract the lat long and case stuy sites
with open(Casestudy2_csv_path, mode='r') as infile:
reader = csv.reader(infile)
next(reader, None)
with open('coors_new.csv', mode='w') as outfile:
writer = csv.writer(outfile)
if Location == 'Estuary':
mydict = dict((rows[0],[rows[1],rows[2]]) for rows in reader)
if Location == 'Ocean':
mydict = dict((rows[0],[rows[3],rows[4]]) for rows in reader)
if Location == 'Catchment':
mydict = dict((rows[0],[rows[5],rows[6]]) for rows in reader)
if Clim_var_type == 'tasmean':
silo_df = sil.pointdata(["max_temp", "min_temp"], 'Okl9EDxgS2uzjLWtVNIBM5YqwvVcCxOmpd3nCzJh', Base_period_start.replace("-", ""), Base_period_end.replace("-", ""),
None, mydict[Estuary][0], mydict[Estuary][1], False, None)
#take mean of daily min and max temp
Present_day_df = silo_df.iloc[:,[2,4]].mean(axis=1)
else:
silo_df = sil.pointdata(SILO_Clim_var, 'Okl9EDxgS2uzjLWtVNIBM5YqwvVcCxOmpd3nCzJh', Base_period_start.replace("-", ""), Base_period_end.replace("-", ""),
None, mydict[Estuary][0], mydict[Estuary][1], False, None)
Present_day_df = silo_df.iloc[:,[2]]
#set the x and y limit deltas - for plotting only
if Clim_var_type in ['evspsblmean', 'potevpmean']: #ET time series that we have is not in the same format as the other variables, hence the different treatment
[minplotDelta, maxplotDelta]=[50,50]
#for tasmean, observed min and max T need to be converted into mean T
elif Clim_var_type == 'tasmean':
[minplotDelta, maxplotDelta]=[0.2,1]
elif Clim_var_type == 'tasmax':
[minplotDelta, maxplotDelta]=[1,2]
elif Clim_var_type == 'wssmean' or Clim_var_type == 'wss1Hmaxtstep':
[minplotDelta, maxplotDelta]=[1, 1.5]
elif Clim_var_type == 'pracc':
[minplotDelta, maxplotDelta]=[50,100]
elif Clim_var_type in ['rsdsmean', 'rldsmean']:
[minplotDelta, maxplotDelta]=[100,100]
#[minplotDelta, maxplotDelta]=[1,1]
#==========================================================#
#substract a constant from all values to convert from kelvin to celcius (temp)
if Clim_var_type == 'sstmean':
Present_day_df = Present_day_df.iloc[:,0:(len(Present_day_df)-1)]-273.15
Present_day_df.index = Present_day_df.index.to_period('D')
#==========================================================#
#create seasonal sums etc.
#==========================================================#
if Stats == 'maxdaily':
Present_day_df_annual = Present_day_df.resample('A').max()
Present_day_df_annual = Present_day_df_annual.replace(0, np.nan)
Fdf_Seas_means = Present_day_df.resample('Q-NOV').max() #seasonal means
Fdf_Seas_means = Fdf_Seas_means.replace(0, np.nan)
if(Stats[:4] =='days'):
Threshold = int(Stats[-2:])
agg = ('>'+ str(Threshold) + '_count', lambda x: x.gt(Threshold).sum()),
Present_day_df_annual = Present_day_df.resample('A').agg(agg)
Fdf_Seas_means = Present_day_df.resample('Q-NOV').agg(agg) #seasonal means
else:
if Clim_var_type in ['pracc' ,'evspsblmean' ,'potevpmean' ,
'pr1Hmaxtstep' ,'wss1Hmaxtstep',
'rsdsmean', 'rldsmean']:
Present_day_df_annual = Present_day_df.resample('A').sum()
Present_day_df_annual = Present_day_df_annual.replace(0, np.nan)
Fdf_Seas_means = Present_day_df.resample('Q-NOV').sum() #seasonal means
Fdf_Seas_means = Fdf_Seas_means.replace(0, np.nan)
elif Clim_var_type in ['tasmean' ,'tasmax' ,'sstmean']:
Present_day_df_annual = Present_day_df.resample('A').mean()
Present_day_df_monthly = Present_day_df.resample('M').mean()
Present_day_df_weekly = Present_day_df.resample('W').mean()
Fdf_Seas_means = Present_day_df.resample('Q-NOV').mean() #seasonal means
#==========================================================#
#==========================================================#
#Loop through annual and seasons and create a deviation plot for each.
#==========================================================#
times = ['annual', 'DJF', 'MAM', 'JJA','SON']
fig = plt.figure(figsize=(14,8))
delta_all_df = pd.DataFrame()
i=1
for temp in times:
#temp = 'annual'
#subset the ensemble dataframe for the period used:
if temp == 'annual':
Ensemble_Delta_df = Ensemble_Delta_full_df.iloc[:,range(0,2)]
Present_Day_ref_df = Present_day_df_annual
else:
Ensemble_Delta_df = Ensemble_Delta_full_df.filter(regex= temp)
Ensemble_Delta_df.columns = ['near', 'far']
if temp == 'DJF':
Mean_df = Fdf_Seas_means[Fdf_Seas_means.index.quarter==1]
Column_names = ['DJF_near', 'DJF_far']
if temp == 'MAM':
Mean_df = Fdf_Seas_means[Fdf_Seas_means.index.quarter==2]
Column_names = ['MAM_near', 'MAM_far']
if temp == 'JJA':
Mean_df = Fdf_Seas_means[Fdf_Seas_means.index.quarter==3]
Column_names = ['JJA_near', 'JJA_far']
if temp == 'SON':
Mean_df = Fdf_Seas_means[Fdf_Seas_means.index.quarter==4]
Present_Day_ref_df = Mean_df
#Subset to present day variability period
Present_Day_ref_df = pd.DataFrame(Present_Day_ref_df.loc[(Present_Day_ref_df.index >= Base_period_start) & (Present_Day_ref_df.index <= Base_period_end)])
#Subset to present day variability delta base period (this is the statistical baseline used for present day conditions)
Present_Day_Delta_ref_df = pd.DataFrame(Present_Day_ref_df.loc[(Present_Day_ref_df.index >= Base_period_DELTA_start) & (Present_Day_ref_df.index <= Base_period_end)])
Present_Day_Mean = np.percentile(Present_Day_Delta_ref_df, 50)
Present_Day_SD = np.std(Present_Day_Delta_ref_df)
#create data frame for floating stacked barplots
index=['-2std', '-1std', 'Med', '1std', '2std']
columns = ['present','near future', 'far future']
Plot_in_df = pd.DataFrame(index=index, columns =columns)
#
Plot_in_df['present'] = [float(Present_Day_Mean-2*Present_Day_SD),float(Present_Day_Mean-Present_Day_SD), float(Present_Day_Mean),
float(Present_Day_Mean+Present_Day_SD), float(Present_Day_Mean+2*Present_Day_SD)]
Plot_in_df['near future'] = [float(Present_Day_Mean + Ensemble_Delta_df.near[-3:][0]),np.NaN, float(Present_Day_Mean + Ensemble_Delta_df.near[-3:][1]),
np.NaN, float(Present_Day_Mean + Ensemble_Delta_df.near[-3:][2])]
Plot_in_df['far future'] = [float(Present_Day_Mean + Ensemble_Delta_df.far[-3:][0]),np.NaN, float(Present_Day_Mean + Ensemble_Delta_df.far[-3:][1]),
np.NaN, float(Present_Day_Mean + Ensemble_Delta_df.far[-3:][2])]
#Create a second data frame that has the values in a way that they can be stacked up in bars with the correct absolute values
Plot_in_df2 = pd.DataFrame(index=index, columns =columns )
#
if presentdaybar == False:
index=['-2std', '-1std', 'Med', '1std', '2std']
columns = ['near future', 'far future']
Plot_in_df2 = pd.DataFrame(index=index, columns =columns )
#Plot_in_df2['present'] = [float(0),float(0), float(0),
# float(0), float(0)]
else:
Plot_in_df2['present'] = [float(Present_Day_Mean-2*Present_Day_SD),float(Present_Day_SD), float(Present_Day_SD),
float(Present_Day_SD), float(Present_Day_SD)]
Plot_in_df2['near future'] = [float(Present_Day_Mean + Ensemble_Delta_df.near[-3:][0]),np.NaN, float(Ensemble_Delta_df.near[-3:][1]-Ensemble_Delta_df.near[-3:][0]),
np.NaN, float(Ensemble_Delta_df.near[-3:][2]-Ensemble_Delta_df.near[-3:][1])]
Plot_in_df2['far future'] = [float(Present_Day_Mean + Ensemble_Delta_df.far[-3:][0]),np.NaN, float(Ensemble_Delta_df.far[-3:][1]-Ensemble_Delta_df.far[-3:][0]),
np.NaN, float(Ensemble_Delta_df.far[-3:][2]-Ensemble_Delta_df.far[-3:][1])]
#transpose the data frame
Plot_in_df_tp = Plot_in_df2.transpose()
#do the individual plots
if temp == 'annual':
xmin = int(min(Plot_in_df.min(axis=1))-minplotDelta)
xmax = int(max(Plot_in_df.max(axis=1))+maxplotDelta)
else:
xmin = int(min(Plot_in_df.min(axis=1))-minplotDelta)
xmax = int(max(Plot_in_df.max(axis=1))+maxplotDelta)
#define colour scheme
#likert_colors = ['none', 'firebrick','firebrick','lightcoral','lightcoral']
likert_colors = ['none', 'darkblue', 'darkblue','cornflowerblue','cornflowerblue']
uni_colors = ['none', 'cornflowerblue', 'cornflowerblue','cornflowerblue','cornflowerblue']
#plot the stacked barplot
if temp == 'annual':
ax=plt.subplot(2,4,3)
else:
ax=plt.subplot(2,4,i)
Plot_in_df_tp.plot.bar(stacked=True, color=uni_colors, edgecolor='none', legend=False, ax=ax, width=0.5)
df = pd.DataFrame(Plot_in_df.iloc[2,[1,2]])
index2=['Med']
columns2 = ['near future', 'far future']
Plot_in_df3 = pd.DataFrame(index=index2, columns =columns2)
Plot_in_df3['near future'] = df.iloc[1,0]
Plot_in_df3['far future'] = df.iloc[0,0]
Plot_in_df3 = Plot_in_df3.transpose()
plt.plot(Plot_in_df3['Med'], linestyle="", markersize=52,
marker="_", color='darkblue', label="Median")
z = plt.axhline(float(Present_Day_Mean-2*Present_Day_SD), linestyle='-', color='black', alpha=.5)
z.set_zorder(-1)
z = plt.axhline(float(Present_Day_Mean+2*Present_Day_SD), linestyle='-', color='black', alpha=.5)
z.set_zorder(-1)
z = plt.axhline(float(Present_Day_Mean-Present_Day_SD), linestyle='--', color='black', alpha=.5)
z.set_zorder(-1)
z = plt.axhline(float(Present_Day_Mean+Present_Day_SD), linestyle='--', color='black', alpha=.5)
z.set_zorder(-1)
z = plt.axhline(float(Present_Day_Mean), linestyle='--', color='red', alpha=.5)
z.set_zorder(-1)
plt.ylim(xmin, xmax)
plt.title(Clim_var_type + ' ' + temp )
ax.grid(False)
for tick in ax.get_xticklabels():
tick.set_rotation(0)
fig.tight_layout()
fig.patch.set_alpha(ALPHA_figs)
if temp == 'annual':
ax=plt.subplot(2,2,1)
Present_Day_ref_df.plot(legend=False, ax=ax)
z = plt.axhline(float(Present_Day_Mean-2*Present_Day_SD), linestyle='-', color='black', alpha=.5)
z.set_zorder(-1)
z = plt.axhline(float(Present_Day_Mean+2*Present_Day_SD), linestyle='-', color='black', alpha=.5)
z.set_zorder(-1)
z = plt.axhline(float(Present_Day_Mean-Present_Day_SD), linestyle='--', color='black', alpha=.5)
z.set_zorder(-1)
z = plt.axhline(float(Present_Day_Mean+Present_Day_SD), linestyle='--', color='black', alpha=.5)
z.set_zorder(-1)
z = plt.axhline(float(Present_Day_Mean), linestyle='--', color='red', alpha=.5)
z.set_zorder(-1)
#fig.patch.set_facecolor('deepskyblue')
fig.tight_layout()
fig.patch.set_alpha(ALPHA_figs)
plt.title(Clim_var_type + ' ' + Stats + ' ' + temp +' present day')
plt.ylim(xmin, xmax)
ax.grid(False)
#if temp == 'MAM':
i=i+4
else:
i=i+1
#plt.show()
if presentdaybar == False:
out_file_name = Estuary + '_' + Clim_var_type + '_' + Stats + '_' + PD_Datasource + '_' + SILO_Clim_var[0] + Version + '_' + '_NPDB.png'
else:
out_file_name = Estuary + '_' + Clim_var_type + '_' + Stats + '_' + PD_Datasource + '_' + SILO_Clim_var[0] + Version + '_' + '.png'
out_path = output_directory + '/' + out_file_name
fig.savefig(out_path)
if present_day_plot == 'yes':
#print present day climate data
fig = plt.figure(figsize=(5,4))
ax = fig.add_subplot(1, 1, 1)
if temp == 'annual':
xmin = int(min(Plot_in_df.min(axis=1))-minplotDelta)
xmax = int(max(Plot_in_df.max(axis=1))+maxplotDelta)
else:
xmin = int(min(Plot_in_df.min(axis=1))-minplotDelta)
xmax = int(max(Plot_in_df.max(axis=1))+maxplotDelta)
Present_Day_ref_df.plot(legend=False, ax=ax)
z = plt.axhline(float(Present_Day_Mean-2*Present_Day_SD), linestyle='-', color='black', alpha=.5)
z.set_zorder(-1)
z = plt.axhline(float(Present_Day_Mean+2*Present_Day_SD), linestyle='-', color='black', alpha=.5)
z.set_zorder(-1)
z = plt.axhline(float(Present_Day_Mean-Present_Day_SD), linestyle='--', color='black', alpha=.5)
z.set_zorder(-1)
z = plt.axhline(float(Present_Day_Mean+Present_Day_SD), linestyle='--', color='black', alpha=.5)
z.set_zorder(-1)
z = plt.axhline(float(Present_Day_Mean), linestyle='--', color='red', alpha=.5)
z.set_zorder(-1)
#fig.patch.set_facecolor('deepskyblue')
fig.patch.set_alpha(0)
plt.ylim(13, xmax)
#export plot to png
out_file_name = Estuary + '_' + Clim_var_type + '_' + Stats + '_' + Base_period_start + '_' + Base_period_end + Version + 'Present_Day_Period.png'
out_path = output_directory + '/' + out_file_name
fig.savefig(out_path)
# use transparent=True if you want the whole figure with a transparent background

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

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

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

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

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

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

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

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

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

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

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

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

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

@ -19,41 +19,64 @@ setwd("C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdoc/P
#Set inputs #Set inputs
###################### ######################
Case.Study <- "CASESTUDY2" Case.Study <- "CASESTUDY2"
Estuary <- "HUNTER" Estuary <- "NADGEE"
Climvar <- 'tasmean' Climvar <- 'tasmean'
ggplotGAM.k <- 7 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 = readLines(AirT_CSV_Path)
dat = as.data.frame(do.call(rbind, strsplit(dat, split=" {2,10}")), stringsAsFactors=FALSE) dat = as.data.frame(do.call(rbind, strsplit(dat, split=" {2,10}")), stringsAsFactors=FALSE)
colnames(dat) <-dat[3,] colnames(dat) <-dat[3,]
dat2 = dat[-c(1:3), ] if (Estuary=='NADGEE'){
dat2 = dat2[-(720:726),] dat2 = dat[-c(1:11), ]
dat2 = dat2[-(1:108),] dat2 = dat2[-(365:381),]
dat2 = dat2[,-1] 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") dat2$Date <- as.yearmon(dat2[,1], "%m %Y")
SeaLev.df <- dat2 SeaLev.df <- dat2
head(SeaLev.df) head(SeaLev.df)
SeaLev.df$MSL <- as.numeric(SeaLev.df$Mean) 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) SeaLev.df$Julday1 <- seq(1,length(SeaLev.df[,1]),1)
linear.trend.model_EC_all <- lm(MSL ~ Julday1, SeaLev.df) linear.trend.model_EC_all <- lm(MSL ~ Julday1, SeaLev.df)
SeaLev.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4] SeaLev.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4]
SeaLev.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 12 SeaLev.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 12
###################### ######################
#Plot #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") ylab("Monthly Mean Sea Level [m]") + xlab("Time")
#export to png #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) png(file = png.name, width = 10.5, height = 7, units='in', res=500)
grid.arrange(p1air,ncol=1) grid.arrange(p1air,ncol=1)
dev.off() 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") ylab("Monthly Mean Sea Level [m]") + xlab("Time")
#export to png #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) png(file = png.name, width = 10.5, height = 7, units='in', res=500)
grid.arrange(p1air,ncol=1) grid.arrange(p1air,ncol=1)
dev.off() dev.off()
Write Preview
Loading…
Cancel
Save