# -*- coding: utf-8 -*- #==========================================================# #Last Updated - June 2018 #@author: z5025317 Valentin Heimhuber #code for creating climate prioritization plots for NARCLIM variables. #Inputs: Uses CSV files that contain all 12 NARCLIM model runs time series for 1 grid cell created with: P1_NARCliM_NC_to_CSV_CCRC_SS.py #==========================================================# #Load packages #==========================================================# import numpy as np import os import pandas as pd import glob import matplotlib import matplotlib.pyplot as plt from datetime import datetime from datetime import timedelta from matplotlib.backends.backend_pdf import PdfPages from ggplot import * matplotlib.style.use('ggplot') #==========================================================# #Input parameters #==========================================================# # Set working direcotry (where postprocessed NARClIM data is located) #set beginning and end years and corresponding scenario code fs=['Hwq003', 'Hwq005'] comparison_code = '003v005' startyear=1999 endyear=1999 year=range(startyear, endyear+1) #set directory path for in and output files output_directory = 'C:/Users/z5025317/OneDrive - UNSW/Hunter_CC_Modeling/08_Results/Output/Postprocessed' + comparison_code +'/' nodes_csv = 'H:/WRL_Projects/Hunter_CC_Modeling/Module_6/03_Results/Chainages/Hunter_nodes.csv' Flow_csv = 'H:/WRL_Projects/Hunter_CC_Modeling/Module_6/01_Input/BCGeneration/Scenarios/Calibration/Calibration_Greta.csv' png_output_directory = output_directory + '/Figures/' variables = ['Sal'] #Tem or Sal #read csv file with nodes to extract data from #node=[837, 666, 981, 59,527, 34, 1, 391] node = pd.read_csv(nodes_csv)['Hunter'].values chainages = pd.read_csv(nodes_csv)['x_km'].values #set plot parameters ALPHA_figs = 1 font = {'family' : 'sans-serif', 'weight' : 'normal', 'size' : 14} matplotlib.rc('font', **font) #==========================================================# #==========================================================# #100% automated part of the code doing the data extraction #========================================================== #output_directory = 'J:/Project wrl2016032/NARCLIM_Raw_Data/Extracted' if not os.path.exists(output_directory): os.makedirs(output_directory) #output_directory = 'J:/Project wrl2016032/NARCLIM_Raw_Data/Extracted' if not os.path.exists(png_output_directory): os.makedirs(png_output_directory) Summary_df = pd.DataFrame() for f in fs: for NODE in node: NODE = str(NODE) #set input and output directories input_directory = 'C:/Users/z5025317/OneDrive - UNSW/Hunter_CC_Modeling/08_Results/Output/'+ f #set directory path for individual png files png_output_directory = output_directory + '/Figures/' # Set working direcotry (where postprocessed NARClIM data is located) os.chdir(input_directory) #==========================================================# #Load data file Clim_Var_CSVs = glob.glob('*'+ NODE + '*WQ*') clim_var_csv_path = Clim_Var_CSVs[0] df = pd.read_csv(clim_var_csv_path, index_col=False, sep=' ') df.index = pd.to_datetime(df.Year, format = '%Y') + pd.to_timedelta(df.Hour, unit='h') df= df.drop(columns=['Year', 'Hour']) df.columns = [NODE+'_Sal', NODE+'_Tem'] #append to summary data frame if f in ['HCC010', 'HCC042']: scname = 'Present' if f == 'HCC011': scname = 'Near_Future' df.index = df.index - (datetime.strptime('2025 01 01', '%Y %m %d').date() - datetime.strptime('1995 01 01', '%Y %m %d').date()) if f == 'HCC012': scname = 'Far_Future' df.index = df.index - (datetime.strptime('2065 01 01', '%Y %m %d').date() - datetime.strptime('1995 01 01', '%Y %m %d').date()) if f in ['HCC040', 'HCC041']: scname = 'Present_Resto' df.columns = df.columns+'_'+ scname Summary_df = pd.concat([Summary_df, df], axis=1, join='outer') #cut down the summary df to common years Summary_df = Summary_df[datetime.strptime(str(startyear) + ' 01 05', '%Y %m %d').date():datetime.strptime(str(endyear) + ' 12 31', '%Y %m %d').date()] Summary_df.tail() for variable in variables: ii=1 for NODE in node: NODE = str(NODE) #=========# #comparison time series plot at node #=========# #out name png_out_file_name = comparison_code + '_' + variable + '_' + NODE + '_time_series.png' png_out_path = png_output_directory + '/' + png_out_file_name # plotin_df = Summary_df.filter(regex=variable) plotin_df = plotin_df.filter(regex=NODE) #prepare data frame for plot fig = plt.figure(figsize=(14,18)) ax=plt.subplot(4,1,1) #start and end dat of plots plt.title(variable + '_' + NODE + ' - time series') start_date = datetime(int('1996'), int('02'), int('01')) end_date = datetime(int('1996'), int('02'), int('5')) plotin_df1=plotin_df[start_date:end_date] ymin = min(plotin_df1.min()) - 0.01 *min(plotin_df1.min()) ymax = max(plotin_df1.max()) + 0.01 *max(plotin_df1.max()) plotin_df1.plot(ylim=(ymin,ymax) ,legend=False, ax=ax) #plt.legend( loc=1) plt.xticks([]) ax=plt.subplot(4,1,2) start_date = datetime(int('1996'), int('02'), int('01')) end_date = datetime(int('1996'), int('02'), int('3')) plotin_df2=plotin_df[start_date:end_date] ymin = min(plotin_df2.min()) - 0.1 *min(plotin_df2.min()) ymax = max(plotin_df2.max()) + 0.1 *max(plotin_df2.max()) plotin_df2.plot(ylim=(ymin,ymax),legend=False, ax=ax) #plt.legend( loc=1) ax.xaxis.grid(False) ax=plt.subplot(4,1,3) plt.title(variable + '_' + NODE + ' - monthly maximum') Monthly_max = plotin_df.resample('M').max() Monthly_max.plot(ax=ax,legend=False) plt.legend( loc=1) ax.xaxis.grid(False) ax=plt.subplot(4,1,4) plt.title(variable + '_' + NODE + ' - monthly mean') Monthly_mean = plotin_df.resample('M').mean() Monthly_mean.plot(ax=ax) #plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, # ncol=2, mode="expand", borderaxespad=0.) #plt.legend( loc=1) ax.xaxis.grid(False) #ax.patch.set_alpha(ALPHA_figs) fig.patch.set_alpha(ALPHA_figs) fig.tight_layout() fig.savefig(png_out_path) plt.close() #=========# png_out_file_name = comparison_code + '_' + variable + '_' + NODE + '_exceedance.png' png_out_path = png_output_directory + '/' + png_out_file_name fig = plt.figure(figsize=(10,5)) ax=plt.subplot(1,1,1) start_date = datetime(int('1995'), int('02'), int('01')) end_date = datetime(int('1998'), int('12'), int('31')) plotin_df3=plotin_df[start_date:end_date] Qhist = plotin_df3.filter(regex='Present') col = Qhist.columns[0] Qhist = Qhist.sort_values(by=[col],ascending=False).values Qnearf = plotin_df3.filter(regex='Near_Future') col = Qnearf.columns[0] Qnearf = Qnearf.sort_values(by=[col],ascending=False).values Qmodern = plotin_df3.filter(regex='Far_Future') col = Qmodern.columns[0] Qmodern = Qmodern.sort_values(by=[col],ascending=False).values Nh = len(Qhist) # number of data points Nm = len(Qmodern) # here's one place types matter: must divide by float not int plt.plot(np.arange(Nh)/float(Nh)*100, Qhist) plt.plot(np.arange(Nm)/float(Nm)*100, Qnearf) plt.plot(np.arange(Nm)/float(Nm)*100, Qmodern) #plt.plot(np.sort(Qmodern)[::-1]) plt.xlabel('% Exceedance') if variable == 'Sal': plt.ylabel('Salinity [PSU]') if variable == 'Tem': plt.ylabel('Temperature [C]') plt.legend(['Present', 'Near Future', 'Far Future']) ax.xaxis.grid(False) #ax.patch.set_alpha(ALPHA_figs) fig.patch.set_alpha(ALPHA_figs) fig.tight_layout() fig.savefig(png_out_path) plt.close()