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

Compare commits

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

No commits in common. '9eda19058ec68183bb67cd0f92826be077638274' and 'c5f2f4f743a590b155e2aaebd3de661556b454c7' have entirely different histories.
9eda19058e ... c5f2f4f743

9 changed files with 517 additions and 1142 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 <- 'rldsmean' Clim_Var <- 'sstmean'
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/CASESTUDY2/' Directory <- 'C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdoc/Projects/Paper#1/Data/NARCLIM_Site_CSVs/'
Filename <- 'CASESTDUY2_NARCLIM_Point_Sites.csv' Filename <- '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"))

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

@ -1,14 +1,13 @@
# -*- 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
@ -19,340 +18,232 @@ 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')
# Load my own functions #plt.rcParams.update(plt.rcParamsDefault)
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/')
#==========================================================# #
#####################################----------------------------------
#loop through case study estuaries #set input parameters
Estuaries = ['HUNTER', 'RICHMOND', 'NADGEE', 'SHOALHAVEN', 'GEORGES','CATHIE'] Base_period_start = '1986-01-01'
Estuaries = ['HUNTER'] Base_period_end = '2005-01-01' #use last day that's not included in period as < is used for subsetting
for Est in Estuaries: Estuary = 'Nadgee' # 'Belongil'
Clim_var_type = "*" # '*' will create pdf for all variables in folder
#==========================================================# Clim_var_type = "pracc*" # '*' will create pdf for all variables in folder
#set input parameters Present_Day_Clim_Var = 'Rainfall' #MaxT, MinT, Rainfall, ET Wind
#==========================================================# present_day_plot = 'yes'
Case_Study_Name = 'CASESTUDY2' Version = "V2"
Casestudy2_csv_path = "Data/NARCLIM_Site_CSVs/CASESTUDY2/CASESTDUY2_NARCLIM_Point_Sites.csv" Stats = 'dailymax' #maximum takes the annual max Precipitation instead of the sum
#Estuary = 'HUNTER' # 'Belongil' #####################################----------------------------------
Estuary = Est # 'Belongil'
print Estuary #set directory path for output files
Base_period_start = '1970-01-01' #Start of interval for base period of climate variability output_directory = 'Output/Clim_Deviation_Plots/'+ Estuary
Base_period_DELTA_start = '1990-01-01' #Start of interval used for calculating mean and SD for base period (Delta base period) #output_directory = 'J:/Project wrl2016032/NARCLIM_Raw_Data/Extracted'
Base_period_end = '2009-01-01' #use last day that's not included in period as < is used for subsetting if not os.path.exists(output_directory):
Clim_var_type = 'tasmax' #Name of climate variable in NARCLIM models '*' will create pdf for all variables in folder os.makedirs(output_directory)
print('-------------------------------------------')
PD_Datasource = 'SILO' #Source for present day climate data (historic time series) can be either: 'Station' or 'SILO' print("output directory folder didn't exist and was generated")
SILO_Clim_var = ['max_temp'] #select the SILO clim variable to be used for base period. - see silo.py for detailed descriptions print('-------------------------------------------')
Location = 'Estuary' # pick locaiton for extracting the SILO data can be: Estuary, Catchment, or Ocean print('-------------------')
Clim_Var_CSVs = glob.glob('./Output/' + Estuary + '/' + Estuary + '_' + Clim_var_type[:-1] + '_' + Stats + '*')
presentdaybar = False #include a bar for present day variability in the plots? #read CSV file
present_day_plot = 'no' #save a time series of present day clim_var_csv_path = Clim_Var_CSVs[0]
Version = "V1" Filename = os.path.basename(os.path.normpath(clim_var_csv_path))
Stats = 'days_h_35' #'maxdaily' #maximum takes the annual max Precipitation instead of the sum Clim_var_type = Filename.split('_', 2)[1]
ALPHA_figs = 1 #Set alpha of figure background (0 makes everything around the plot panel transparent) 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
#set directory path for output files Present_day_Var_CSV = glob.glob('./Data/Wheather_Station_Data/**/' + Estuary + '_' + Present_Day_Clim_Var + '*csv')
output_directory = 'Output/' + Case_Study_Name + '/' + Estuary + '/' + '/Clim_Deviation_Plots/' Present_day_df = pd.read_csv(Present_day_Var_CSV[0])
#output_directory = 'J:/Project wrl2016032/NARCLIM_Raw_Data/Extracted' Dates = pd.to_datetime(Present_day_df.Date)
if not os.path.exists(output_directory): Present_day_df.index = Dates
os.makedirs(output_directory) Present_day_df = Present_day_df.iloc[:,1]
print('-------------------------------------------') Present_day_df = Present_day_df.replace(r'\s+', np.nan, regex=True)
print("output directory folder didn't exist and was generated") Present_day_df = pd.to_numeric(Present_day_df)
print('-------------------------------------------') Present_day_df.index = Dates
print('-------------------') [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])
#read CSV file 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
Clim_Var_CSVs = glob.glob('./Output/' + Case_Study_Name + '/' + Estuary + '/' + Estuary + '_' + Clim_var_type + '_' + Stats + '*') Present_day_MaxT_df = Present_day_df.iloc[:,5]
clim_var_csv_path = Clim_Var_CSVs[0] Present_day_Var_CSV = glob.glob('./Data/Wheather_Station_Data/**/' + Estuary + '_MinT*csv')
Filename = os.path.basename(os.path.normpath(clim_var_csv_path)) Present_day_df = pd.read_csv(Present_day_Var_CSV[0])
Clim_var_type = Filename.split('_', 2)[1] Dates = pd.to_datetime(Present_day_df.Year*10000+Present_day_df.Month*100+Present_day_df.Day,format='%Y%m%d')
print(Clim_var_type) Present_day_df.index = Dates
Ensemble_Delta_full_df = pd.read_csv(clim_var_csv_path, index_col=0, parse_dates = True) Present_day_MinT_df = Present_day_df.iloc[:,5]
Present_day_df = (Present_day_MaxT_df + Present_day_MinT_df)/2
if Clim_var_type == 'rsdsmean': [minplotDelta, maxplotDelta]=[1,2]
Clim_Var_CSVs = glob.glob('./Output/' + Case_Study_Name + '/' + Estuary + '/' + Estuary + '_rsdsmean_' + Stats + '*') elif Clim_var_type == 'tasmax':
clim_var_csv_path = Clim_Var_CSVs[0] [minplotDelta, maxplotDelta]=[1,2]
Filename = os.path.basename(os.path.normpath(clim_var_csv_path)) elif Clim_var_type == 'wssmean' or Clim_var_type == 'wss1Hmaxtstep':
Clim_var_type = Filename.split('_', 2)[1] Present_day_Var_CSV = glob.glob('./Data/Wheather_Station_Data/**/' + Estuary + '_' + Present_Day_Clim_Var + '*csv')
print(Clim_var_type) Present_day_df = pd.read_csv(Present_day_Var_CSV[0])
Ensemble_Delta_full_df1 = pd.read_csv(clim_var_csv_path, index_col=0, parse_dates = True) Present_day_df.index = Present_day_df[['Year', 'Month', 'Day', 'Hour']].apply(lambda s : datetime(*s),axis = 1)
Clim_Var_CSVs = glob.glob('./Output/' + Case_Study_Name + '/' + Estuary + '/' + Estuary + '_rldsmean_' + Stats + '*') Present_day_df = Present_day_df.filter(regex= 'm/s')
clim_var_csv_path = Clim_Var_CSVs[0] Present_day_df = Present_day_df.replace(r'\s+', np.nan, regex=True)
Filename = os.path.basename(os.path.normpath(clim_var_csv_path)) Present_day_df['Wind speed measured in m/s'] = Present_day_df['Wind speed measured in m/s'].convert_objects(convert_numeric=True)
Clim_var_type = Filename.split('_', 2)[1] [minplotDelta, maxplotDelta]=[1, 1.5]
print(Clim_var_type) else:
Ensemble_Delta_full_df2 = pd.read_csv(clim_var_csv_path, index_col=0, parse_dates = True) Present_day_Var_CSV = glob.glob('./Data/Wheather_Station_Data/**/' + Estuary + '_' + Present_Day_Clim_Var + '*csv')
Ensemble_Delta_full_df = Ensemble_Delta_full_df2 + Ensemble_Delta_full_df1 Present_day_df = pd.read_csv(Present_day_Var_CSV[0])
#Ensemble_Delta_full_df = pd.to_numeric(Ensemble_Delta_full_df) 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.
#load present day climate variable 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 PD_Datasource == 'Station': if Stats == 'dailymax':
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.resample('A').max()
Present_day_df_annual = Present_day_df_annual.replace(0, np.nan) 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: else:
if Clim_var_type in ['pracc' ,'evspsblmean' ,'potevpmean' , Present_day_df_annual = Present_day_df.resample('A').sum()
'pr1Hmaxtstep' ,'wss1Hmaxtstep', Present_day_df_annual = Present_day_df_annual.replace(0, np.nan)
'rsdsmean', 'rldsmean']: Fdf_Seas_means = Present_day_df.resample('Q-NOV').sum() #seasonal means
Present_day_df_annual = Present_day_df.resample('A').sum() Fdf_Seas_means = Fdf_Seas_means.replace(0, np.nan)
Present_day_df_annual = Present_day_df_annual.replace(0, np.nan) else:
Fdf_Seas_means = Present_day_df.resample('Q-NOV').sum() #seasonal means Present_day_df_annual = Present_day_df.resample('A').mean()
Fdf_Seas_means = Fdf_Seas_means.replace(0, np.nan) Present_day_df_monthly = Present_day_df.resample('M').mean()
elif Clim_var_type in ['tasmean' ,'tasmax' ,'sstmean']: Present_day_df_weekly = Present_day_df.resample('W').mean()
Present_day_df_annual = Present_day_df.resample('A').mean() Fdf_Seas_means = Present_day_df.resample('Q-NOV').mean() #seasonal means
Present_day_df_monthly = Present_day_df.resample('M').mean()
Present_day_df_weekly = Present_day_df.resample('W').mean() #Loop through annual and seasons and create a deviation plot for each.
Fdf_Seas_means = Present_day_df.resample('Q-NOV').mean() #seasonal means times = ['annual', 'DJF', 'MAM', 'JJA','SON']
#==========================================================# fig = plt.figure(figsize=(14,8))
delta_all_df = pd.DataFrame()
i=1
#==========================================================# for temp in times:
#Loop through annual and seasons and create a deviation plot for each. temp = 'annual'
#==========================================================# #subset the ensemble dataframe for the period used:
times = ['annual', 'DJF', 'MAM', 'JJA','SON'] if temp == 'annual':
fig = plt.figure(figsize=(14,8)) Ensemble_Delta_df = Ensemble_Delta_full_df.iloc[:,range(0,2)]
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: Present_Day_ref_df = Present_day_df_annual
index=['-2std', '-1std', 'Med', '1std', '2std'] else:
columns = ['near future', 'far future'] Ensemble_Delta_df = Ensemble_Delta_full_df.filter(regex= temp)
Plot_in_df2 = pd.DataFrame(index=index, columns =columns ) Ensemble_Delta_df.columns = ['near', 'far']
#Plot_in_df2['present'] = [float(0),float(0), float(0), if temp == 'DJF':
# float(0), float(0)] Mean_df = Fdf_Seas_means[Fdf_Seas_means.index.quarter==1]
else: Column_names = ['DJF_near', 'DJF_far']
Plot_in_df2['present'] = [float(Present_Day_Mean-2*Present_Day_SD),float(Present_Day_SD), float(Present_Day_SD), if temp == 'MAM':
float(Present_Day_SD), float(Present_Day_SD)] Mean_df = Fdf_Seas_means[Fdf_Seas_means.index.quarter==2]
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]), Column_names = ['MAM_near', 'MAM_far']
np.NaN, float(Ensemble_Delta_df.near[-3:][2]-Ensemble_Delta_df.near[-3:][1])] if temp == 'JJA':
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]), Mean_df = Fdf_Seas_means[Fdf_Seas_means.index.quarter==3]
np.NaN, float(Ensemble_Delta_df.far[-3:][2]-Ensemble_Delta_df.far[-3:][1])] Column_names = ['JJA_near', 'JJA_far']
#transpose the data frame if temp == 'SON':
Plot_in_df_tp = Plot_in_df2.transpose() Mean_df = Fdf_Seas_means[Fdf_Seas_means.index.quarter==4]
#do the individual plots Present_Day_ref_df = Mean_df
if temp == 'annual': print(Ensemble_Delta_df.columns.values)
xmin = int(min(Plot_in_df.min(axis=1))-minplotDelta) #Subset to present day variability period
xmax = int(max(Plot_in_df.max(axis=1))+maxplotDelta) 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)])
else: Present_Day_Mean = np.percentile(Present_Day_ref_df, 50)
xmin = int(min(Plot_in_df.min(axis=1))-minplotDelta) Present_Day_SD = np.std(Present_Day_ref_df)
xmax = int(max(Plot_in_df.max(axis=1))+maxplotDelta) #create data frame for floating stacked barplots
#define colour scheme index=['-2std', '-1std', 'Med', '1std', '2std']
#likert_colors = ['none', 'firebrick','firebrick','lightcoral','lightcoral'] columns = ['present','near future', 'far future']
likert_colors = ['none', 'darkblue', 'darkblue','cornflowerblue','cornflowerblue'] Plot_in_df = pd.DataFrame(index=index, columns =columns)
uni_colors = ['none', 'cornflowerblue', 'cornflowerblue','cornflowerblue','cornflowerblue'] #
#plot the stacked barplot Plot_in_df['present'] = [float(Present_Day_Mean-2*Present_Day_SD),float(Present_Day_Mean-Present_Day_SD), float(Present_Day_Mean),
if temp == 'annual': float(Present_Day_Mean+Present_Day_SD), float(Present_Day_Mean+2*Present_Day_SD)]
ax=plt.subplot(2,4,3) 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]),
else: np.NaN, float(Present_Day_Mean + Ensemble_Delta_df.near[-3:][2])]
ax=plt.subplot(2,4,i) 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]),
Plot_in_df_tp.plot.bar(stacked=True, color=uni_colors, edgecolor='none', legend=False, ax=ax, width=0.5) np.NaN, float(Present_Day_Mean + Ensemble_Delta_df.far[-3:][2])]
df = pd.DataFrame(Plot_in_df.iloc[2,[1,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
index2=['Med'] Plot_in_df2 = pd.DataFrame(index=index, columns =columns )
columns2 = ['near future', 'far future']
Plot_in_df3 = pd.DataFrame(index=index2, columns =columns2) Plot_in_df2['present'] = [float(Present_Day_Mean-2*Present_Day_SD),float(Present_Day_SD), float(Present_Day_SD),
Plot_in_df3['near future'] = df.iloc[1,0] float(Present_Day_SD), float(Present_Day_SD)]
Plot_in_df3['far future'] = df.iloc[0,0] 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]),
Plot_in_df3 = Plot_in_df3.transpose() np.NaN, float(Ensemble_Delta_df.near[-3:][2]-Ensemble_Delta_df.near[-3:][1])]
plt.plot(Plot_in_df3['Med'], linestyle="", markersize=52, 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]),
marker="_", color='darkblue', label="Median") np.NaN, float(Ensemble_Delta_df.far[-3:][2]-Ensemble_Delta_df.far[-3:][1])]
z = plt.axhline(float(Present_Day_Mean-2*Present_Day_SD), linestyle='-', color='black', alpha=.5) #transpose the data frame
z.set_zorder(-1) Plot_in_df_tp = Plot_in_df2.transpose()
z = plt.axhline(float(Present_Day_Mean+2*Present_Day_SD), linestyle='-', color='black', alpha=.5) #do the individual plots
z.set_zorder(-1) if temp == 'annual':
z = plt.axhline(float(Present_Day_Mean-Present_Day_SD), linestyle='--', color='black', alpha=.5) xmin = int(min(Plot_in_df.min(axis=1))-minplotDelta)
z.set_zorder(-1) xmax = int(max(Plot_in_df.max(axis=1))+maxplotDelta)
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: else:
out_file_name = Estuary + '_' + Clim_var_type + '_' + Stats + '_' + PD_Datasource + '_' + SILO_Clim_var[0] + Version + '_' + '.png' xmin = int(min(Plot_in_df.min(axis=1))-minplotDelta)
xmax = int(max(Plot_in_df.max(axis=1))+maxplotDelta)
#define colour scheme
#likert_colors = ['none', 'firebrick','firebrick','lightcoral','lightcoral']
likert_colors = ['none', 'darkblue', 'darkblue','cornflowerblue','cornflowerblue']
#plot the stacked barplot
ax=plt.subplot(2,3,i)
Plot_in_df_tp.plot.bar(stacked=True, color=likert_colors, edgecolor='none', legend=False, ax=ax)
z = plt.axhline(float(Present_Day_Mean-2*Present_Day_SD), linestyle='-', color='black', alpha=.5)
z.set_zorder(-1)
z = plt.axhline(float(Present_Day_Mean+2*Present_Day_SD), linestyle='-', color='black', alpha=.5)
z.set_zorder(-1)
z = plt.axhline(float(Present_Day_Mean-Present_Day_SD), linestyle='--', color='black', alpha=.5)
z.set_zorder(-1)
z = plt.axhline(float(Present_Day_Mean+Present_Day_SD), linestyle='--', color='black', alpha=.5)
z.set_zorder(-1)
z = plt.axhline(float(Present_Day_Mean), linestyle='--', color='red', alpha=.5)
z.set_zorder(-1)
plt.ylim(xmin, xmax)
plt.title(Clim_var_type + ' ' + temp )
ax.grid(False)
for tick in ax.get_xticklabels():
tick.set_rotation(0)
fig.tight_layout()
fig.patch.set_alpha(0)
#reset i to i+1 for next step
if temp == 'MAM':
i=i+2
else:
i=i+1
print(i)
plt.show()
out_file_name = Estuary + '_' + Clim_var_type + '_' + Stats + '_' + '_CC_prio_plot' + Version + '.png'
out_path = output_directory + '/' + out_file_name
fig.savefig(out_path)
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)
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 out_path = output_directory + '/' + out_file_name
fig.savefig(out_path) fig.savefig(out_file_name)
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) # use transparent=True if you want the whole figure with a transparent background
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

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

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

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

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

@ -38,11 +38,12 @@ 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 = 'Estuary' #'Catchment' Location = '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 + '/'
@ -84,7 +85,7 @@ for Estuary in mydict:
print('-------------------------------------------') print('-------------------------------------------')
#==========================================================# #==========================================================#
output_csv = output_directory_internal + 'SILO_climdata_' + Estuary +'_'+ Location +'_' + mydict[Estuary][0] + '_' + mydict[Estuary][1] + '.csv' output_csv = output_directory_internal + 'SILO_climdata_' + Estuary +'_'+ Location +'_' + mydict[Estuary][0] + '_' + mydict[Estuary][1] + '2.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,6 +40,8 @@ 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)
@ -71,7 +73,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, Stats): def calculate_deltas_NF_FF2(Annual_df, Seasonal_df):
"""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']
@ -92,11 +94,9 @@ def calculate_deltas_NF_FF2(Annual_df, Seasonal_df, Stats):
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']
if(Stats[:4] =='days'): models = list(Seasonal_df.mean().index)
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,88 +123,4 @@ def calculate_deltas_NF_FF2(Annual_df, Seasonal_df, Stats):
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)
if(Stats[:4] =='days'): return delta_all_df
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

@ -1,165 +0,0 @@
#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_level.R → Analysis/Code/ggplot_time_series_with_trends_sea_levle.R
Unescape Escape View File

@ -19,64 +19,41 @@ setwd("C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdoc/P
#Set inputs #Set inputs
###################### ######################
Case.Study <- "CASESTUDY2" Case.Study <- "CASESTUDY2"
Estuary <- "NADGEE" Estuary <- "HUNTER"
Climvar <- 'tasmean' Climvar <- 'tasmean'
ggplotGAM.k <- 7 ggplotGAM.k <- 7
###################### ######################
###################### ######################
Output.Directory <- paste('./Output/', Case.Study, '/', Estuary,'/Recent_Trends/', sep="") #Set input file paths
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,]
if (Estuary=='NADGEE'){ dat2 = dat[-c(1:3), ]
dat2 = dat[-c(1:11), ] dat2 = dat2[-(720:726),]
dat2 = dat2[-(365:381),] dat2 = dat2[-(1:108),]
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
###################### ######################
@ -91,7 +68,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.Directory, '/Trends_MonthlyMeanSeaLevel_full_period_', Sys.Date(),".png", sep="") png.name <- paste('./Output/', Case.Study, '/', Estuary, '/Trends_MonthlyMeanSeaLevel_full_period_', Sys.Date(),".png", sep="")
png(file = png.name, width = 10.5, height = 7, units='in', res=500) 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()
@ -107,7 +84,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.Directory, '/Trends_MonthlyMeanSeaLevel_MultiGAM_full_period_', Sys.Date(),".png", sep="") png.name <- paste('./Output/', Case.Study, '/', Estuary, '/Trends_MonthlyMeanSeaLevel_MultiGAM_full_period_', Sys.Date(),".png", sep="")
png(file = png.name, width = 10.5, height = 7, units='in', res=500) 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