You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
215 lines
8.8 KiB
Python
215 lines
8.8 KiB
Python
# -*- coding: utf-8 -*-
|
|
"""
|
|
Created on Mon Nov 25 10:56:21 2019
|
|
|
|
@author: z5025317
|
|
"""
|
|
|
|
|
|
#==========================================================#
|
|
#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
|
|
import seaborn
|
|
matplotlib.style.use('ggplot')
|
|
#========================================================#
|
|
|
|
|
|
#========================================================#
|
|
#set input parameters and paths
|
|
#========================================================#
|
|
Gauges = [ 'Stockton', 'Raymond', 'Oakhampton', 'Paterson'] #name of all gauges or IDs in a given river
|
|
River = 'Clarence'
|
|
|
|
#name of plot scenario/experiment . This is a unique identifier for this particular scenario comparison
|
|
Plot_scenario = 'Clarence_gauge_analysis_V1'
|
|
|
|
#time frames for analysis
|
|
tideanalysis_day = '2006 06 07' # the time series plots will be centred around this day.
|
|
year = 2006
|
|
|
|
#set directories
|
|
Folder = 'C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdoc/Projects/Paper#2_HD_modeling/01_Data/Tides/CSV/' #needs to be adjusted - should be a folder with all the MHL data in individual CSV files
|
|
Out_Folder = 'C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdoc/Projects/Paper#2_HD_modeling/02_Analysis/Tides/' #needs to be adjusted
|
|
#========================================================#
|
|
|
|
|
|
|
|
#========================================================#
|
|
#laod MHL data and store in a single dataframe
|
|
#========================================================#
|
|
WL_df = pd.DataFrame()
|
|
for gauge in Gauges:
|
|
# gauge = 'Oyster'
|
|
print 'loading data for ' + gauge
|
|
CSVpath = glob.glob(Folder + '*' + gauge + '*' ) # will find the csv file corresponding to the gauge name or ID
|
|
print(CSVpath)
|
|
Full_df = pd.read_csv(CSVpath ,header= 0, skiprows= 30, parse_dates = True) #needs to be adjusted to fit the MHL format
|
|
#Full_df.index = pd.to_datetime(Full_df.Date, format = '%d/%m/%Y') + pd.to_timedelta(Full_df.Time, unit='h') #time series index may have to be adjusted, depending on the outcome of parse_dates
|
|
Full_df.index = pd.to_datetime(Full_df.Date, format = '%d/%m/%Y') + pd.to_timedelta(Full_df.Time, unit='h')
|
|
Full_df = Full_df.iloc[:,[-2]]
|
|
Full_df.columns = ['WL[AHD]']
|
|
Full_df = Full_df.loc[(Full_df.index >= year + '-01-01 00:00:00') & (Full_df.index <= str(np.int(year)+1) + '-01-01 00:00:00')]
|
|
Full_df = Full_df.loc[~Full_df.index.duplicated(keep='first')] #get rid of possible duplicate datetimes
|
|
WL_df = pd.concat([WL_df, Full_df], axis=1)
|
|
WL_df.columns = Gauges
|
|
#subset the data in time if desired (e.g. to a given year)
|
|
WL_df_1 = WL_df.loc[(WL_df.index >= '2017-01-10 00:00:00') & (WL_df.index <= '2017-01-18 00:00:00')]
|
|
#better
|
|
#WL_df_1 = WL_df.loc[(WL_df.index >= datetime.strptime(startdate , '%Y %m %d').date()) & (WL_df.index <= datetime.strptime(enddate , '%Y %m %d').date())] #may need to be adjusted
|
|
|
|
#export dataframe to CSV if desired
|
|
#export the pandas data frame as a CSV file within the output directory is desired
|
|
#WL_df.to_csv(Folder + 'data_combined.csv')
|
|
#========================================================#
|
|
|
|
|
|
|
|
|
|
#========================================================#
|
|
# Plot water level histogram at all gauges on a singel plot
|
|
#========================================================#
|
|
#full period density comparison all four gauges
|
|
png_out_name = Out_Folder + River + '_water level histograms.pdf'
|
|
fig = plt.figure(figsize=(40,18))
|
|
ax=plt.subplot(1,1,1)
|
|
plt.title('Hunter water level density comparison - Jan 2016- Feb 2020')
|
|
WL_df.iloc[:,[0,1,2,3]].plot.kde(xlim=(-1.5,1.5), ax=ax)
|
|
fig.tight_layout()
|
|
fig.savefig(png_out_name)
|
|
plt.close()
|
|
fig.tight_layout()
|
|
fig.savefig(png_out_name)
|
|
plt.close()
|
|
#========================================================#
|
|
|
|
|
|
|
|
|
|
#========================================================#
|
|
# Plot water level histogram at all gauges
|
|
#========================================================#
|
|
ylims = [-0.9,2.2]
|
|
#set plot parameters
|
|
ALPHA_figs = 1
|
|
font = {'family' : 'sans-serif',
|
|
'weight' : 'normal',
|
|
'size' : 6}
|
|
matplotlib.rc('font', **font)
|
|
|
|
#plot
|
|
i=1
|
|
png_out_name = Out_Folder + '/' + Plot_scenario + '/' + River + 'B_water level histogram and exceedance planes '+ startdate + '_' + enddate + '3.pdf'
|
|
fig = plt.figure(figsize=(len(Gauges)*4,8))
|
|
for Gauge in Gauges:
|
|
#Gauge = 'Yamba'
|
|
xlocation = 0.6
|
|
#full period density comparison for individual gauges overlaid by tidal planes
|
|
Wollow_13ExPr = np.nanpercentile(WL_df[Gauge],87)
|
|
Wollow_6p4ExPr = np.nanpercentile(WL_df[Gauge],93.6)
|
|
Wollow_1_3ExPr = np.nanpercentile(WL_df[Gauge],98.7)
|
|
Wollow_0p01ExPr = np.nanpercentile(WL_df[Gauge],99.99)
|
|
Wollow_50ExPr = np.nanpercentile(WL_df[Gauge],50)
|
|
ax=plt.subplot(4,1,i)
|
|
plt.title( Gauge + ' water level density comparison - Jan 2016- Feb 2020')
|
|
#WL_df[Gauge].plot.kde(xlim=(-0.5,1),vertical=True, lw=3, ax=ax)
|
|
seaborn.kdeplot(WL_df[Gauge], vertical=True, color='royalblue', lw=3, ax=ax)
|
|
|
|
plt.axhline(y=Wollow_13ExPr, color='firebrick', linestyle='--', lw=3, alpha=1)
|
|
plt.text( xlocation , (Wollow_13ExPr), '13% exceedance', ha='left', va='bottom')
|
|
|
|
plt.axhline(y=Wollow_6p4ExPr, color='orange', linestyle='--', lw=3, alpha=1)
|
|
plt.text(xlocation , (Wollow_6p4ExPr), '6.4% exceedance', ha='left', va='bottom')
|
|
|
|
plt.axhline(y=Wollow_1_3ExPr, color='yellow', linestyle='--', lw=3, alpha=1)
|
|
plt.text( xlocation ,Wollow_1_3ExPr , '1.3% exceedance', ha='left', va='bottom')
|
|
|
|
plt.axhline(y=Wollow_0p01ExPr, color='limegreen', linestyle='--', lw=3, alpha=1)
|
|
plt.text( xlocation , Wollow_0p01ExPr , '0.01% exceedance', ha='left', va='bottom')
|
|
|
|
plt.axhline(y=Wollow_50ExPr, color='cyan', linestyle='--', lw=3, alpha=1)
|
|
plt.text( xlocation , Wollow_50ExPr , '50% exceedance', ha='left', va='bottom')
|
|
|
|
plt.ylim(-1,1.7)
|
|
|
|
plt.xlabel("Probability/density")
|
|
plt.ylabel("Water level [m AHD]")
|
|
|
|
fig.tight_layout()
|
|
fig.savefig(png_out_name)
|
|
plt.close()
|
|
fig.tight_layout()
|
|
fig.savefig(png_out_name)
|
|
i=i+1
|
|
plt.close()
|
|
#========================================================#
|
|
|
|
|
|
|
|
#========================================================#
|
|
# Save exceedance planes to csv
|
|
#========================================================#
|
|
Summary={}
|
|
Summary['13% exceedance'] = Wollow_13ExPr
|
|
Summary['6.4% exceedance'] = Wollow_6p4ExPr
|
|
Summary['1.3% exceedance'] = Wollow_1_3ExPr
|
|
Summary['0.01% exceedance'] = Wollow_0p01ExPr
|
|
Summary['50% exceedance'] = Wollow_50ExPr
|
|
pdf1 = pd.DataFrame(Summary, index=[0])
|
|
pdf1.to_csv(os.path.join( Out_Folder + Gauge + '_entrance_stats.csv'))
|
|
#========================================================#
|
|
|
|
|
|
|
|
#========================================================#
|
|
# Plot water level time series at control points
|
|
#========================================================#
|
|
#set plot parameters
|
|
ALPHA_figs = 1
|
|
font = {'family' : 'sans-serif',
|
|
'weight' : 'normal',
|
|
'size' : 5}
|
|
matplotlib.rc('font', **font)
|
|
|
|
halfwindows = [1,2,5] #Days to add before and after central date chosen for analysis
|
|
for halfwindow in halfwindows:
|
|
png_out_name = Out_Folder + '/' + Plot_scenario + '/' +River+ 'water level time series '+ startdate + '_' + enddate + ' _window_of_' + str(halfwindow*2) +'_days.pdf'
|
|
fig = plt.figure(figsize=(10,len(Gauges)*3))
|
|
i=1
|
|
for Gauge in Gauges:
|
|
Gauge_only_df = WL_df.filter(regex=Gauge)
|
|
#Gauge_only_df_1 = Gauge_only_df.loc[(Gauge_only_df.index >= '2008-05-09 12:00:00') & (Gauge_only_df.index <= '2008-05-20 12:00:00')]
|
|
Gauge_only_df_1 = Gauge_only_df[datetime.strptime(tideanalysis_day , '%Y %m %d').date()- timedelta(days=halfwindow) : datetime.strptime(tideanalysis_day , '%Y %m %d').date() + timedelta(days=halfwindow) ]
|
|
ax=plt.subplot(len(Gauges),1,i)
|
|
plt.title('Tidal water levels at 4 locations in the ' + River +' Estuary')
|
|
ax.set_prop_cycle(color=['steelblue', 'coral', 'steelblue', 'deeppink'],linestyle=[ '-', '-','--','-'], lw=[1,1,1,2])
|
|
Gauge_only_df_1.iloc[:,[0,1,2]].plot( legend=False, ax=ax)
|
|
plt.axhline(y=0, color='slateblue', linestyle='--', lw=1, alpha=0.5)
|
|
#plt.axhline(y=0.9, color='tomato', linestyle='--', lw=1, alpha=0.5)
|
|
plt.ylim(ylims)
|
|
plt.ylabel("Water level [m AHD]")
|
|
plt.xlabel("Time")
|
|
ax.legend()
|
|
fig.tight_layout()
|
|
i=i+1
|
|
fig.savefig(png_out_name)
|
|
plt.close()
|
|
#========================================================#
|
|
# Plot water level histograms at control points
|
|
#========================================================#
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|