#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Standardised graphics generating functions
Created on Wed Aug 17 15:54:43 2022
@author: chris
"""
# ----------------------------------------------------------------------------
# IMPORTS
# ----------------------------------------------------------------------------
# Standard Python Dependencies
import copy
# Non-Standard Python Dependencies
import numpy as np
from matplotlib import tri
import matplotlib.pyplot as plt
# Local Module Dependencies
# Other Dependencies
# ----------------------------------------------------------------------------
# GLOBAL VARIABLES
# ----------------------------------------------------------------------------
# ----------------------------------------------------------------------------
# CLASS DEFINITIONS
# ----------------------------------------------------------------------------
#--------------------------------------------------------------------------
# FUNCTION DEFINITIONS
#--------------------------------------------------------------------------
[docs]def plotMesh(x,y,elems,figsize=(9,6),UTM=True):
"""
Plot an unstrcutured triangular mesh from node coordinates and elements list.
Parameters
----------
x : numpy.array, float
x-coordinate of mesh nodes.
y : numpy.array, float
y-coordinate of mesh nodes.
elems : numpy.array, integer
triangular element corner nodes.
figsize : tuple, float, optional
Figure size in inches. The default is (9,6).
UTM : bool, optional
Flag for CRS of node coordinates (True=UTM or False=WGS84) . The default is True.
Returns
-------
figHnd : obj
Handle to the figure object created.
"""
msh = tri.Triangulation(x,y,elems)
figHnd = plt.figure(figsize=figsize)
plt.triplot(msh)
ax = plt.gca()
ax.set_aspect(1)
if UTM:
plt.xlabel('X [ m ]')
plt.ylabel('Y [ m ]')
else:
plt.xlabel('Lon. [ '+r'$^{\circ}$'+' E ]')
plt.ylabel('Lat. [ '+r'$^{\circ}$'+' N ]')
return figHnd
[docs]def plot_validation_stats(obs,mod,metrics,paramStr,unitStr,min_val,max_val,info):
# Plot results
figHnd = plt.figure(figsize=[10,6])
plt.plot(obs,mod,'k.',markersize=1)
plt.xlim([min_val,max_val])
plt.ylim([min_val,max_val])
plt.plot((min_val,max_val),(min_val,max_val),'--',color='gray',linewidth=1)
plt.title('Modelled vs Observed '+paramStr)
plt.xlabel('Observed [ '+unitStr+' ]')
plt.ylabel('Modelled [ '+unitStr+' ]')
ax = plt.gca()
ax.set_aspect('equal')
pos = ax.get_position()
ax.set_position([0.1,0.08,pos.x1-pos.x0,pos.y1])
plt.grid('on')
dr = (max_val-min_val)
x = max_val + 0.1*dr
hdfontsize = 12
txfontsize = 8
plt.text(x,max_val-0.00*dr,'Validation Statistics',font='DejaVu Sans Mono',fontsize=hdfontsize)
plt.text(x,max_val-0.05*dr,'Correlation = '+"{0:+.3f}".format(metrics['R']),font='DejaVu Sans Mono',fontsize=txfontsize)
plt.text(x,max_val-0.10*dr,'Goodness-of-fit = '+"{0:+.3f}".format(metrics['R2']),font='DejaVu Sans Mono',fontsize=txfontsize)
plt.text(x,max_val-0.15*dr,'Mean Bias = '+"{0:+.3f}".format(metrics['MB'])+unitStr,font='DejaVu Sans Mono',fontsize=txfontsize)
plt.text(x,max_val-0.20*dr,'Normalised MB = '+"{0:+.2f}".format(metrics['NMB'])+'%',font='DejaVu Sans Mono',fontsize=txfontsize)
plt.text(x,max_val-0.25*dr,'Mean Absolute Error = '+"{0:+.3f}".format(metrics['MAE'])+unitStr,font='DejaVu Sans Mono',fontsize=txfontsize)
plt.text(x,max_val-0.30*dr,'Normalised MAE = '+"{0:+.2f}".format(metrics['NMAE'])+'%',font='DejaVu Sans Mono',fontsize=txfontsize)
plt.text(x,max_val-0.35*dr,'RMS Error = '+"{0:+.3f}".format(metrics['RMSE'])+unitStr,font='DejaVu Sans Mono',fontsize=txfontsize)
plt.text(x,max_val-0.40*dr,'Normalised RMSE = '+"{0:+.2f}".format(metrics['NRMSE'])+'%',font='DejaVu Sans Mono',fontsize=txfontsize)
plt.text(x,max_val-0.45*dr,'Scatter Index = '+"{0:+.2f}".format(metrics['SI'])+'%',font='DejaVu Sans Mono',fontsize=txfontsize)
plt.text(x,max_val-0.50*dr,'Number of Samples = '+str(np.int32(metrics['NSMPL'])),font='DejaVu Sans Mono',fontsize=txfontsize)
plt.text(x,max_val-0.60*dr,'Data Description',font='DejaVu Sans Mono',fontsize=hdfontsize)
plt.text(x,max_val-0.65*dr,'Observation Data Source: '+info['Obs_Src'],font='DejaVu Sans Mono',fontsize=txfontsize)
plt.text(x,max_val-0.70*dr,'Deployment ID: '+info['Obs_Deploy'],font='DejaVu Sans Mono',fontsize=txfontsize)
plt.text(x,max_val-0.75*dr,'Model Data Source: '+info['Mod_Src'],font='DejaVu Sans Mono',fontsize=txfontsize)
plt.text(x,max_val-0.80*dr,'Data Geo Location: '+str(info['DataLoc']),font='DejaVu Sans Mono',fontsize=txfontsize)
plt.text(x,max_val-0.85*dr,'Data CRS: '+info['CRS'],font='DejaVu Sans Mono',fontsize=txfontsize)
if info['hab']:
plt.text(x,max_val-0.90*dr,'Data Height Above Bed: '+str(info['zloc'])+' m',font='DejaVu Sans Mono',fontsize=txfontsize)
else:
plt.text(x,max_val-0.90*dr,'Depth Below Surface: '+str(info['zloc'])+' m',font='DejaVu Sans Mono',fontsize=txfontsize)
plt.text(x,max_val-0.95*dr,'Coverage: '+str(info['NumDays'])+' days from '+info['StartDate'],font='DejaVu Sans Mono',fontsize=txfontsize)
plt.text(x,max_val-1.00*dr,'Tidal State: '+info['tideState'],font='DejaVu Sans Mono',fontsize=txfontsize)
return figHnd
[docs]def flow_class_diagram(data,tec=None,velmax=None):
# Extract Information
dataset = data['dataset']
datasrc = data['datasrc']
deploystr = data['deploystr']
hab = data['hab']
zloc = data['zloc']
thta = data['thta'].copy()
magn = data['magn'].copy()
res_thta = data['res_thta']
res_magn = data['res_magn']
pkdir = data['peak_dir'].copy()
peak_id = data['peak_id'].copy()
if tec is not None:
tec_name = tec.model
cut_in_speed = tec.cutInSpeed
rated_speed = tec.ratedSpeed
# IEC Power Curve Limits
IEC_PC_lower = 0.5 * cut_in_speed
IEC_PC_upper = 1.2 * rated_speed
# Set default velmax value if not given
if velmax is None:
velmax = 4.5
# Get Flood and Ebb peak angles
fld_id = peak_id['flood']
ebb_id = peak_id['ebb']
# Magnitude vs Direction Diagram
figHnd = plt.figure(figsize=(11.5,6.5))
plt.xlim([-180.0,180.0])
plt.ylim([0.0,velmax])
plt.xticks(np.arange(-180.0, 180.0, 30.0))
plt.grid(True)
ax = plt.gca()
ax.set_position([0.08,0.12,0.86,0.75])
xl = plt.xlim()
yl = plt.ylim()
# Overlay Flood Peak Information
fld_dir = pkdir['dir'][fld_id]
fld_bearing = pkdir['bearing'][fld_id]
plt.plot(np.asarray([1.0,1.0])*fld_dir,yl,'k-')
plt.text(fld_dir-2.0,velmax-0.2,'FLOOD',fontsize=8, \
horizontalalignment='right',verticalalignment='center')
plt.text(fld_dir+2.0,velmax-0.18,'[ '+"{0:+.1f}".format(fld_dir)+r'$^{\circ}$'+' ]',fontsize=8, \
horizontalalignment='left',verticalalignment='center')
plt.text(fld_dir-2.0,velmax-0.4,'( Bearing',fontsize=8, \
color=(0.5,0.5,0.5), \
horizontalalignment='right',verticalalignment='center')
plt.text(fld_dir+2.0,velmax-0.38,"{0:+.1f}".format(fld_bearing)+r'$^{\circ}$'+' )',fontsize=8, \
color=(0.5,0.5,0.5), \
horizontalalignment='left',verticalalignment='center')
#Overlay Ebb Peak Information
ebb_dir = pkdir['dir'][ebb_id]
ebb_bearing = pkdir['bearing'][ebb_id]
plt.plot(np.asarray([1.0,1.0])*ebb_dir,yl,'k-')
plt.text(ebb_dir-2.0,velmax-0.2,'EBB',fontsize=8, \
horizontalalignment='right',verticalalignment='center')
plt.text(ebb_dir+2.0,velmax-0.18,'[ '+"{0:+.1f}".format(ebb_dir)+r'$^{\circ}$'+' ]',fontsize=8, \
horizontalalignment='left',verticalalignment='center')
plt.text(ebb_dir-2.0,velmax-0.4,'( Bearing',fontsize=8, \
color=(0.5,0.5,0.5), \
horizontalalignment='right',verticalalignment='center')
plt.text(ebb_dir+2.0,velmax-0.38,"{0:+.1f}".format(ebb_bearing)+r'$^{\circ}$'+' )',fontsize=8, \
color=(0.5,0.5,0.5), \
horizontalalignment='left',verticalalignment='center')
# Plot heat map of magnitude vs direction
#dens = (240,300)
dens = (160,200)
#dens = (120,150)
dtheta = 360/dens[0]
dvmagn = velmax/dens[1]
h, XBinEdges, YBinEdges = np.histogram2d(thta,magn,dens,([-180.0,180.0],[0.0,velmax]),density=True)
pdf2d = h
pdf2d[pdf2d==0.0] = np.nan
hcmap = copy.copy(plt.cm.get_cmap('viridis'))
hcmap.set_bad(alpha=0)
heatmap = plt.imshow(pdf2d.transpose(),interpolation='nearest',origin='lower', \
extent=[XBinEdges[0], XBinEdges[-1], YBinEdges[0], YBinEdges[-1]], \
aspect='auto',cmap=hcmap)
cb = plt.colorbar(heatmap, fraction=0.05, aspect=30, pad=0.1, shrink=1.0)
cb.ax.tick_params(labelsize=8)
cb.ax.tick_params(axis='y', labelrotation=90)
cblabel = 'Probability Density [ '+'$\Delta$'+'$\\theta$'+' = '+ \
"{0:+.2f}".format(dtheta)+r'$^{\circ}$'+' , '+'$\Delta\it{v}$'+ \
' = '+"{0:+.2f}".format(dvmagn)+' m s'+r'$^{-1}$'+']'
cb.set_label(cblabel,fontsize=8,rotation=90)
cb.ax.set_position([0.92125, 0.12, 0.9364166666666666, 0.75])
ax.set_position([0.08,0.12,0.82,0.75])
textx = pkdir['mean_dir']
# Overlay Cut-in speed, rated speed, IEC Power Curve upper and lower thresholds
if tec is not None:
plt.plot(xl,np.asarray([1.0,1.0])*cut_in_speed,color=(0,0.5,0),linestyle='--')
plt.text(textx,cut_in_speed+0.05,'CUT IN SPEED ['+tec_name+']',color=(0,0.5,0),fontsize=8, \
horizontalalignment='center',verticalalignment='bottom')
plt.plot(xl,np.asarray([1.0,1.0])*rated_speed,color=(0,0.5,0),linestyle='--')
plt.text(textx,rated_speed+0.05,'RATED SPEED ['+tec_name+']',color=(0,0.5,0),fontsize=8, \
horizontalalignment='center',verticalalignment='bottom')
plt.plot(xl,np.asarray([1.0,1.0])*IEC_PC_lower,color=(0.3,0,0),linestyle='-')
plt.text(textx,IEC_PC_lower+0.05,'IEC PC Lower Limit',color=(0.3,0,0),fontsize=8, \
horizontalalignment='center',verticalalignment='bottom')
plt.plot(xl,np.asarray([1.0,1.0])*IEC_PC_upper,color=(0.3,0,0),linestyle='-')
plt.text(textx,IEC_PC_upper+0.05,'IEC PC Upper Limit',color=(0.3,0,0),fontsize=8, \
horizontalalignment='center',verticalalignment='bottom')
# Add further information
deltheta = pkdir['del_dir']
plt.text(textx,2.0,'$\Delta$ $\\theta$ = '+"{0:+.1f}".format(deltheta)+r'$^{\circ}$',fontsize=8, \
horizontalalignment='center',verticalalignment='top')
'''
plt.plot(res_thta, res_magn,color=(0.6,0.0,0.0),marker='o',markersize=3,markerfacecolor=(0.6,0.0,0.0))
plt.text(pkdir['mean_dir'], 2.2, \
'Mean Residual Flow [ '+"{0:+.1f}".format(res_thta)+r'$^{\circ}$'+' , '+"{0:+.1f}".format(res_magn)+' '+r'$m s^{-1}$'+' ]', \
color=(0.6,0.0,0.0),fontsize=8, \
horizontalalignment='center',verticalalignment='bottom')
'''
plt.text(0,velmax-2.5,'EAST',color=(0.5,0.5,0.5),fontsize=8,rotation=90, \
horizontalalignment='center',verticalalignment='center')
plt.text(90,velmax-2.5,'NORTH',color=(0.5,0.5,0.5),fontsize=8,rotation=90, \
horizontalalignment='center',verticalalignment='center')
plt.text(180,velmax-2.5,'WEST',color=(0.5,0.5,0.5),fontsize=8,rotation=90, \
horizontalalignment='center',verticalalignment='center')
plt.text(-180,velmax-2.5,'WEST',color=(0.5,0.5,0.5),fontsize=8,rotation=90, \
horizontalalignment='center',verticalalignment='center')
plt.text(-90,velmax-2.5,'SOUTH',color=(0.5,0.5,0.5),fontsize=8,rotation=90, \
horizontalalignment='center',verticalalignment='center')
plt.xlabel('Flow Direction [$\circ$, cartesian coordinates]',fontsize=9)
plt.ylabel('Flow Speed [m/s]',fontsize=9)
if hab:
plt.title('Flow Classification - '+dataset+' [ '+datasrc+' ] at HAB = '+"{0:+.1f}".format(zloc)+'m\n' + deploystr + '\n')
else:
plt.title('Flow Classification - '+dataset+' [ '+datasrc+' ] at DBS = '+"{0:+.1f}".format(-zloc)+'m\n' + deploystr + '\n')
return figHnd
[docs]def speed_bin_histograms(data,tec=None,binwidth=0.2,maxfrac=0.2,maxspd=4.5):
dspd = binwidth
fraclim = maxfrac
txtx = fraclim * 0.875
# Extract Information
dataset = data['dataset']
magn = data['magn'].copy()
fldindices = data['flood_indices']
ebbindices = data['ebb_indices']
fldpc = data['flood_total_percent']
fldpctage = data['flood_percent']
ebbpc = data['ebb_total_percent']
ebbpctage = data['ebb_percent']
if tec is not None:
tec_name = tec.model
cut_in_speed = tec.cutInSpeed
rated_speed = tec.ratedSpeed
# Percentages of time breakdown
x1 = cut_in_speed/2.0
x2 = (cut_in_speed+rated_speed)/2.0
x3 = (rated_speed+maxspd)/2.0
figHnd = plt.figure(figsize=(12,5))
nvals = len(fldindices)
wghts = np.ones(fldindices.shape,dtype=float)/nvals
# Flood
(fldax,ebbax) = figHnd.subplots(1,2)
fldax.hist(magn[fldindices],np.arange(0,maxspd,dspd),weights=wghts,orientation='horizontal',color=[0.5,0.5,0.5],edgecolor='black',rwidth=1.0)
fldax.set_position([0.07,0.11,0.43,0.75])
if tec is not None:
fldax.plot([0,fraclim],np.asarray([1,1])*cut_in_speed,color=(0,0.3,0),linestyle='--')
fldax.plot([0,fraclim],np.asarray([1,1])*rated_speed,color=(0,0.3,0),linestyle='--')
fldax.text(txtx,cut_in_speed-0.1,'CUT IN SPEED',color=(0,0.3,0),fontsize=6,\
horizontalalignment='center',verticalalignment='center')
fldax.text(txtx,cut_in_speed+0.1,'['+tec_name+']',color=(0,0.3,0),fontsize=6,\
horizontalalignment='center',verticalalignment='center')
fldax.text(txtx,rated_speed-0.1,'RATED SPEED',color=(0,0.3,0),fontsize=6,\
horizontalalignment='center',verticalalignment='center')
fldax.text(txtx,rated_speed+0.1,'['+tec_name+']',color=(0,0.3,0),fontsize=6,\
horizontalalignment='center',verticalalignment='center')
fldax.set_xlim([0,fraclim])
fldax.set_ylim([0,maxspd])
fldax.grid('on')
if tec is not None:
fldax.text(txtx,x1,"{0:+.1f}".format(fldpctage[0])+'%',fontsize=8,\
horizontalalignment='center',verticalalignment='center')
fldax.text(txtx,x2,"{0:+.1f}".format(fldpctage[1])+'%',fontsize=8,\
horizontalalignment='center',verticalalignment='center')
fldax.text(txtx,x3,"{0:+.1f}".format(fldpctage[2])+'%',fontsize=8,\
horizontalalignment='center',verticalalignment='center')
fldax.set_title('Flood Tide [ '+"{0:.1f}".format(fldpc)+'% of time ]',fontsize=9)
fldax.set_ylabel('Flow Speed [ $\mathregular{m s^{-1}}$ ]',fontsize=9)
fldax.set_xlabel('Fraction of Time',fontsize=9)
tstr = dataset+' - Fraction of Time in Speed Bins [ Bin Width = '+"{0:.1f}".format(binwidth)+' $\mathregular{m s^{-1}}$ ]'
fldax.text(fraclim*1.1,maxspd*1.075,tstr,\
fontsize=10,fontweight='bold',\
horizontalalignment='center',verticalalignment='bottom')
# Ebb
nvals = len(ebbindices)
wghts = np.ones(ebbindices.shape,dtype=float)/nvals
ebbax.hist(magn[ebbindices],np.arange(0,maxspd,dspd),weights=wghts,orientation='horizontal',color=[0.5,0.5,0.5],edgecolor='black',rwidth=1.0)
ebbax.set_position([0.54,0.11,0.43,0.75])
if tec is not None:
ebbax.plot([0,fraclim],np.asarray([1,1])*cut_in_speed,color=(0,0.3,0),linestyle='--')
ebbax.plot([0,fraclim],np.asarray([1,1])*rated_speed,color=(0,0.3,0),linestyle='--')
ebbax.text(txtx,cut_in_speed-0.1,'CUT IN SPEED',color=(0,0.3,0),fontsize=6,\
horizontalalignment='center',verticalalignment='center')
ebbax.text(txtx,cut_in_speed+0.1,'['+tec_name+']',color=(0,0.3,0),fontsize=6,\
horizontalalignment='center',verticalalignment='center')
ebbax.text(txtx,rated_speed-0.1,'RATED SPEED',color=(0,0.3,0),fontsize=6,\
horizontalalignment='center',verticalalignment='center')
ebbax.text(txtx,rated_speed+0.1,'['+tec_name+']',color=(0,0.3,0),fontsize=6,\
horizontalalignment='center',verticalalignment='center')
ebbax.set_xlim([0,fraclim])
ebbax.set_ylim([0,maxspd])
ebbax.grid('on')
if tec is not None:
ebbax.text(txtx,x1,"{0:+.1f}".format(ebbpctage[0])+'%',fontsize=8,\
horizontalalignment='center',verticalalignment='center')
ebbax.text(txtx,x2,"{0:+.1f}".format(ebbpctage[1])+'%',fontsize=8,\
horizontalalignment='center',verticalalignment='center')
ebbax.text(txtx,x3,"{0:+.1f}".format(ebbpctage[2])+'%',fontsize=8,\
horizontalalignment='center',verticalalignment='center')
ebbax.set_title('Ebb Tide [ '+"{0:.1f}".format(ebbpc)+'% of time ]',fontsize=9)
ebbax.set_xlabel('Fraction of Time',fontsize=9)
return figHnd
[docs]def directional_spread_plot(data):
figHnd = plt.figure(figsize=(12,5))
(fldax,ebbax) = figHnd.subplots(1,2)
# Flood
nvals = len(data['del_flood_dir'])
wghts = np.ones(data['del_flood_dir'].shape,dtype=float)/nvals
fldax.hist(data['del_flood_dir'],np.arange(-90,90,2.5),weights=wghts,color=[0.5,0.5,0.5],edgecolor='black')
fldax.set_position([0.07,0.11,0.43,0.75])
fldax.plot([0,0],[0,0.35],'k--')
fldax.set_xlim([-90,90])
fldax.grid('on')
fldax.set_title('Flood Tide [ '+"{0:.1f}".format(data['flood_total_percent'])+'% of time ]',fontsize=9)
fldax.set_xlabel('Flow Direction Anomaly [ $\mathregular{^{\circ}}$ ]',fontsize=9)
fldax.set_ylabel('Fraction in Direction',fontsize=9)
dsetstr = data['dataset']
rng = fldax.get_ylim()
fldax.text(100.0,rng[1]*1.075,\
dsetstr+' - Directional Spread Around Peak Flow [ Bin Width = 2.5$\mathregular{^{\circ}}$ ]',\
fontsize=10,fontweight='bold',\
horizontalalignment='center',verticalalignment='bottom')
flood_pc_lt_5deg = 100.0*len(np.where((data['del_flood_dir'] >= -5.0) & (data['del_flood_dir'] <= 5.0))[0])/len(data['del_flood_dir'])
print('% Flood Direction within 5 degrees of peak = '+"{0:.1f}".format(flood_pc_lt_5deg))
# Ebb
nvals = len(data['del_ebb_dir'])
wghts = np.ones(data['del_ebb_dir'].shape,dtype=float)/nvals
ebbax.hist(data['del_ebb_dir'],np.arange(-90,90,2.5),weights=wghts,color=[0.5,0.5,0.5],edgecolor='black')
ebbax.set_position([0.54,0.11,0.43,0.75])
ebbax.plot([0,0],[0,0.35],'k--')
ebbax.set_xlim([-90,90])
ebbax.grid('on')
ebbax.set_title('Ebb Tide [ '+"{0:.1f}".format(data['ebb_total_percent'])+'% of time ]',fontsize=9)
ebbax.set_xlabel('Flow Direction Anomaly [ $\mathregular{^{\circ}}$ ]',fontsize=9)
ebb_pc_lt_5deg = 100.0*len(np.where((data['del_ebb_dir'] >= -5.0) & (data['del_ebb_dir'] <= 5.0))[0])/len(data['del_flood_dir'])
print('% Ebb Direction within 5 degrees of peak = '+"{0:.1f}".format(ebb_pc_lt_5deg))
return figHnd
#--------------------------------------------------------------------------
# MAIN PROCESS
#--------------------------------------------------------------------------
#--------------------------------------------------------------------------
# INTERFACE
#--------------------------------------------------------------------------