# -*- coding: utf-8 -*-
"""
Functions to generate standardiszed graphics from data generated by various OCMW processes.
"""
#Standardised graphics generating functions
#
#Chris Old
#IES, School of Engineering, University of Edinburgh
#Aug2022
# ----------------------------------------------------------------------------
# IMPORTS
# ----------------------------------------------------------------------------
# Standard Python Dependencies
import copy
# Non-Standard Python Dependencies
import numpy as np
from matplotlib import tri, colors
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import matplotlib.cm as cm
from mpl_toolkits.axes_grid1 import make_axes_locatable
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
# Local Module Dependencies
from .timeFuncs import datetimeFromNum
# Other Dependencies
# ----------------------------------------------------------------------------
# GLOBAL VARIABLES
# ----------------------------------------------------------------------------
parula_vals = [[0.2422 , 0.1504 , 0.6603 , 1. ],
[0.25039048, 0.16499524, 0.70761429, 1. ],
[0.25777143, 0.18178095, 0.7511381 , 1. ],
[0.26472857, 0.19775714, 0.79521429, 1. ],
[0.27064762, 0.21467619, 0.83637143, 1. ],
[0.27511429, 0.2342381 , 0.87098571, 1. ],
[0.2783 , 0.25587143, 0.89907143, 1. ],
[0.28033333, 0.27823333, 0.9221 , 1. ],
[0.2813381 , 0.30059524, 0.94137619, 1. ],
[0.28101429, 0.32275714, 0.95788571, 1. ],
[0.27946667, 0.34467143, 0.97167619, 1. ],
[0.27597143, 0.36668095, 0.98290476, 1. ],
[0.26991429, 0.3892 , 0.9906 , 1. ],
[0.26024286, 0.41232857, 0.99515714, 1. ],
[0.24403333, 0.43583333, 0.99883333, 1. ],
[0.22064286, 0.46025714, 0.99728571, 1. ],
[0.19633333, 0.48471905, 0.98915238, 1. ],
[0.18340476, 0.50737143, 0.97979524, 1. ],
[0.17864286, 0.52885714, 0.96815714, 1. ],
[0.1764381 , 0.54990476, 0.95201905, 1. ],
[0.16874286, 0.5702619 , 0.93587143, 1. ],
[0.154 , 0.5902 , 0.9218 , 1. ],
[0.14602857, 0.60911905, 0.90785714, 1. ],
[0.13802381, 0.62762857, 0.89729048, 1. ],
[0.12481429, 0.64592857, 0.88834286, 1. ],
[0.11125238, 0.6635 , 0.87631429, 1. ],
[0.09520952, 0.67982857, 0.85978095, 1. ],
[0.06887143, 0.69477143, 0.83935714, 1. ],
[0.02966667, 0.70816667, 0.81633333, 1. ],
[0.00357143, 0.72026667, 0.7917 , 1. ],
[0.00665714, 0.73121429, 0.76601429, 1. ],
[0.04332857, 0.74109524, 0.73940952, 1. ],
[0.09639524, 0.75 , 0.7120381 , 1. ],
[0.14077143, 0.7584 , 0.68415714, 1. ],
[0.1717 , 0.7669619 , 0.65544286, 1. ],
[0.19376667, 0.77576667, 0.6251 , 1. ],
[0.21608571, 0.7843 , 0.5923 , 1. ],
[0.24695714, 0.79179524, 0.55674286, 1. ],
[0.29061429, 0.79729048, 0.51882857, 1. ],
[0.34064286, 0.8008 , 0.47885714, 1. ],
[0.3909 , 0.80287143, 0.43544762, 1. ],
[0.44562857, 0.80241905, 0.39091905, 1. ],
[0.5044 , 0.7993 , 0.348 , 1. ],
[0.5615619 , 0.79423333, 0.30448095, 1. ],
[0.61739524, 0.78761905, 0.2612381 , 1. ],
[0.67198571, 0.77927143, 0.2227 , 1. ],
[0.7242 , 0.76984286, 0.19102857, 1. ],
[0.77383333, 0.75980476, 0.16460952, 1. ],
[0.82031429, 0.74981429, 0.15352857, 1. ],
[0.86343333, 0.7406 , 0.15963333, 1. ],
[0.90354286, 0.73302857, 0.17741429, 1. ],
[0.93925714, 0.72878571, 0.20995714, 1. ],
[0.97275714, 0.72977143, 0.23944286, 1. ],
[0.99564762, 0.74337143, 0.23714762, 1. ],
[0.99698571, 0.76585714, 0.21994286, 1. ],
[0.99520476, 0.78925238, 0.2027619 , 1. ],
[0.9892 , 0.81356667, 0.18853333, 1. ],
[0.97862857, 0.83862857, 0.17655714, 1. ],
[0.96764762, 0.8639 , 0.16429048, 1. ],
[0.96100952, 0.88901905, 0.15367619, 1. ],
[0.95967143, 0.91345714, 0.14225714, 1. ],
[0.96279524, 0.9373381 , 0.12650952, 1. ],
[0.96911429, 0.96062857, 0.1063619 , 1. ],
[0.9769 , 0.9839 , 0.0805 , 1. ]]
# ----------------------------------------------------------------------------
# CLASS DEFINITIONS
# ----------------------------------------------------------------------------
#--------------------------------------------------------------------------
# FUNCTION DEFINITIONS
#--------------------------------------------------------------------------
[docs]
def gen_parula_cmap():
"""
Generate matlab parula Colour map for matplotlib
"""
parula = colors.ListedColormap(parula_vals)
return parula
#--------------------------------------
# Model mesh plots
#--------------------------------------
[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 plotMeshField(mesh,vals,figsize=(12.0,4.8),vmin=None,vmax=None,
titleStr="",unitStr="", colormap=None,
show_cbar=False, cbar_width="1.5%",
sub_reg=None, ismesh=False):
"""
Plot a 2D flow parameter defined at either the mesh nodes or element centroids as a colourised mesh field.
"""
x = mesh['meshx']
y = mesh['meshy']
elems = mesh['elems']
centroid = False
if len(vals) > len(x):
centroid = True
if vmin is None:
vmin = np.min(vals)
if vmax is None:
vmax = np.max(vals)
if colormap is None:
if ismesh:
cmap = "Purples"
else:
cmap = "PuOr_r"
else:
cmap = colormap
msh = tri.Triangulation(x,y,elems)
fig = plt.figure(figsize=figsize)
ax = fig.gca()
if sub_reg is None:
divider = make_axes_locatable(ax)
ax_cb = divider.new_horizontal(size=cbar_width, pad=0.1)
fig = ax.get_figure()
plt.title(titleStr)
fig.add_axes(ax_cb)
if centroid:
im = ax.tripcolor(msh,facecolors=vals,cmap=cmap,shading='flat',vmin=vmin,vmax=vmax)
else:
im = ax.tripcolor(msh,vals,cmap=cmap,shading='flat',vmin=vmin,vmax=vmax)
ax.tick_params(axis='both', which='major', labelsize=12)
plt.setp(ax.get_yticklabels(), rotation=90, ha="center", rotation_mode="anchor")
ax.set_aspect(1)
cbar = plt.colorbar(im, cax=ax_cb)
cbar.ax.tick_params(labelsize=12)
cbar.ax.set_ylabel('[ '+unitStr+' ]',fontsize=12)
plt.draw()
else:
plt.title(titleStr)
if centroid:
im = ax.tripcolor(msh,facecolors=vals,cmap=cmap,shading='flat',vmin=vmin,vmax=vmax)
else:
im = ax.tripcolor(msh,vals,cmap=cmap,shading='flat',vmin=vmin,vmax=vmax)
if ismesh:
ax.triplot(msh,linewidth=0.5,color='gray')
ax.set_xlim((sub_reg[0],sub_reg[1]))
ax.set_ylim((sub_reg[2],sub_reg[3]))
ax.set_aspect(1)
plt.setp(ax.get_yticklabels(), rotation=90, ha="center", rotation_mode="anchor")
if show_cbar:
divider = make_axes_locatable(ax)
ax_cb = divider.new_horizontal(size=cbar_width, pad=0.1)
fig.add_axes(ax_cb)
cbar = plt.colorbar(im, cax=ax_cb)
cbar.ax.tick_params(labelsize=12)
cbar.ax.set_ylabel('[ '+unitStr+' ]',fontsize=12)
plt.draw()
if show_cbar or (sub_reg is None):
return fig, ax, ax_cb
else:
return fig, ax
#--------------------------------------
# Flow validation plots
#--------------------------------------
[docs]
def plot_validation_stats(obs,mod,metrics,paramStr,unitStr,min_val,max_val,info):
"""
Plot validation statistics as a correlation plot plus a list of metrics.
"""
# 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
#--------------------------------------
# Flow characterisation plots
#--------------------------------------
[docs]
def flow_class_diagram(data,tec=None,velmax=None):
"""
Plot a flow ccharacterisation diagram based on the preprocessed timeseries data.
"""
# 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 Characterisation - '+dataset+' [ '+datasrc+' ] at HAB = '+"{0:+.1f}".format(zloc)+'m\n' + deploystr + '\n')
else:
plt.title('Flow Characterisation - '+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,gethist=False):
"""
Plot a flow speed bin histograms for flood and ebb tides.
"""
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)
fld_hist = 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)
if binwidth < 0.1:
tstr = dataset+' - Fraction of Time in Speed Bins [ Bin Width = '+"{0:.2f}".format(binwidth)+' $\mathregular{m s^{-1}}$ ]'
else:
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
ebb_hist = 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)
if gethist:
return figHnd, fld_hist, ebb_hist
else:
return figHnd
[docs]
def directional_spread_plot(data):
"""
Plot a flow directional spread histograms for flood and ebb tides.
"""
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
#---------------------------------------------
# WW3 model validation plots
#---------------------------------------------
[docs]
def plot_time_series(t_mod, v_mod, t_obs, v_obs,
obsVar, obsPlatform, modVar, modfname,
year, month,
xlims=None, ylims=None):
"""
Overlaid plot of model and observed time series of specified parameter for WW3 model validation.
"""
i_b = 'INSITU_GLO_WAVE'
obsDates = [datetimeFromNum(t) for t in t_obs]
modDates = [datetimeFromNum(t) for t in t_mod]
figHnd = plt.figure(figsize=(10,5)) # create first figure
label = obsPlatform+' | year: '+str(year)+' | month: '+str(month)
plt.plot(obsDates, v_obs, 'bo-', linewidth=1, markersize=3, label=label)
label= 'RSCD file: '+modfname
plt.plot(modDates, v_mod, 'ro--', linewidth=1, markersize=3, label=label)
plt.legend()
plt.grid(True)
plt.title('Buoy source: '+i_b,fontsize=12)
plt.xlabel('Time (dates)',fontsize=10)
plt.ylabel(obsVar+' / '+modVar+' [m]',fontsize=10)
plt.tight_layout()
#figHnd.show()
return figHnd
[docs]
def plot_correlation(obs, model, obsVar, obsPlatform, modVar, modfname,
year, month, metrics):
"""
Plot a WW3 validation statistics as a correlation plot plus metrics.
"""
figHnd = plt.figure(figsize=(6,6)) # create first figure
sc = plt.scatter(model, obs, s=3)
#plt.colorbar(sc)
if max(obs) > max(model):
max_plot_val = max(obs)
else:
max_plot_val = max(model)
plt.plot([0, max_plot_val+1], [0, max_plot_val+1], color='k', label='1:1')
plt.legend(loc='upper right')
plt.title(obsVar+' vs '+modVar, fontsize=12)
plt.xlabel(modfname, fontsize=10)
plt.ylabel('Buoy: '+obsPlatform, fontsize=10)
ax = plt.gca()
plt.text(0.05,0.9,
'Corr. : '+str(round(metrics['R'], 3))+\
' RMSE : '+str(round(metrics['RMSE'], 3))+' [m]'+\
' NRMSE : '+str(round(metrics['NRMSE'], 3))+' [%]'+\
'\nSI : '+str(round(metrics['SI'], 3))+\
' NMB : '+str(round(metrics['NMB'], 3))+' [%]',\
bbox=dict(fc=(1, 1, 1)), fontsize=10,
transform=ax.transAxes)
plt.grid(True)
plt.tight_layout()
#figHnd.show()
return figHnd
[docs]
def plot_close(figHnd):
plt.close(figHnd)
return
[docs]
def Regional_Weighted(lat: list, lon: list, wght: list, plot_area: list,
clims: list,
plotTitle: str, outpath: str, outname: str,
saveplot: bool):
"""
Regional map of parameters at specific locations colourised by a
wieghting factor.
"""
# plot_area = [-15, 15, 35, 70]
lon_div = [plot_area[0],
plot_area[0]+5,
plot_area[0]+10,
0,
plot_area[1]-10,
plot_area[1]-5,
plot_area[1]]
lat_div = [plot_area[2],
plot_area[2]+5,
plot_area[2]+10,
plot_area[2]+15,
plot_area[3]-15,
plot_area[3]-10,
plot_area[3]-5,
plot_area[3]]
#fig = plt.figure(figsize=(11, 10))
fig = plt.figure(figsize=(9, 10))
cm = plt.cm.get_cmap('YlOrRd')
ax = plt.axes(projection=ccrs.PlateCarree())
ax.add_feature(cfeature.LAND,color='grey',alpha=0.7)
ax.add_feature(cfeature.OCEAN,color='grey',alpha=0.15)
ax.coastlines(resolution='10m') # Options are 110, 50m, or 10m
sc = plt.scatter(lon, lat, c=wght, transform=ccrs.PlateCarree(),
s=80,alpha=0.8, cmap=cm, zorder = 3)
if clims[0] is not None:
plt.clim(clims[0],clims[1])
#plt.colorbar(sc, fraction=0.05, aspect=30, pad=0.1, shrink=0.95)
plt.colorbar(sc, fraction=0.05, aspect=30, pad=0.1, shrink=0.80)
ax.set_extent(plot_area)
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=1, color='gray', alpha=0.8)
gl.xlabel_style = {'size': 12,'color': 'gray'}#, 'weight': 'bold'}
gl.ylabel_style = {'size': 12,'color': 'gray'}#, 'weight': 'bold'}
gl.xlocator = mticker.FixedLocator(lon_div)
gl.ylocator = mticker.FixedLocator(lat_div)
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
plt.title(plotTitle)
plt.tight_layout(pad=5)
if saveplot:
fig.savefig(outpath+'/'+outname+'.png')
return
[docs]
def Regional(lat: list, lon: list,
plot_area: list, lat_div: list, lon_div: list,
outpath: str, outname: str,
saveplot: bool):
"""
Regional map of specific locations.
"""
# plot_area = [-15, 15, 35, 70]
'''
lon_div = [plot_area[0],
plot_area[0]+5,
plot_area[0]+10,
0,
plot_area[1]-10,
plot_area[1]-5,
plot_area[1]]
lat_div = [plot_area[2],
plot_area[2]+5,
plot_area[2]+10,
plot_area[2]+15,
plot_area[3]-15,
plot_area[3]-10,
plot_area[3]-5,
plot_area[3]]
'''
fig = plt.figure(figsize=(9.5, 10))
cm = plt.cm.get_cmap('YlOrRd')
ax = plt.axes(projection=ccrs.PlateCarree())
ax.add_feature(cfeature.LAND,color='grey',alpha=0.7)
ax.add_feature(cfeature.OCEAN,color='grey',alpha=0.15)
ax.coastlines(resolution='10m') # Options are 110, 50m, or 10m
sc = plt.scatter(lon, lat, c='c', transform=ccrs.PlateCarree(),
s=20,alpha=0.8, cmap=cm,zorder = 3)
ax.set_extent(plot_area)
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=1, color='gray', alpha=0.8)
gl.xlabel_style = {'size': 12,'color': 'gray'}#, 'weight': 'bold'}
gl.ylabel_style = {'size': 12,'color': 'gray'}#, 'weight': 'bold'}
gl.xlocator = mticker.FixedLocator(lon_div)
gl.ylocator = mticker.FixedLocator(lat_div)
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
plt.tight_layout(pad=5)
if saveplot:
fig.savefig(outpath+'/'+outname+'.png')
return
[docs]
def plot_1d_wave_spec(freq, spec1D, titleStr:str=None,
outpath:str='./', outname:str='spec_1d', saveplot:bool=False):
"""
Plot of a 1-D wave frequency spectrum
"""
fhnd = plt.figure(figsize=(10,5))
plt.plot(freq,spec1D,'k')
plt.grid()
plt.xlabel('Frequency [Hz]')
plt.ylabel('S(f) [$m^2 s$]')
plt.title(titleStr)
if saveplot:
fhnd.savefig(outpath+'/'+outname+'.png')
return
[docs]
def plot_2d_wave_spec(freq, ntheta, spec2D, fmax:float=0.25, cmax:float=2.0, titleStr:str=None):
"""
Plot of a 2-D wave frequency spectrum
"""
ifrq = np.where(freq <= fmax)[0]
f = freq[ifrq]
nfreq = len(f)
sp = spec2D[ifrq,:]
sp = 100.0 * sp / np.nansum(sp)
sp = np.append(sp,sp[:,0].reshape(nfreq,1),1)
theta = np.linspace(0.0,360.0,ntheta+1)
r, t = np.meshgrid(f, np.radians(theta))
clevels = np.arange(0.1*cmax/2.0,cmax,0.05*cmax/2.0)
fig, ax = plt.subplots(subplot_kw=dict(projection='polar'), figsize=(10,10))
ax.contourf(t, r, sp.transpose(), clevels, cmap=cm.twilight)
clevels = np.arange(0.1*cmax/2.0,cmax,0.05*cmax/2.0)
ax.contour(t, r, sp.transpose(), clevels, colors=['k'], linewidths=0.2)
if titleStr is not None:
plt.title(titleStr)
return fig
[docs]
def plot_speedbinned_velocities(magn3D,z,indices,cntrs):
"""
Plot speed binned velocity data
"""
fig = plt.figure(figsize=(14,7))
for indx in np.arange(cntrs.shape[0]):
if np.mod(indx,3)==0:
binIndx = np.where(indices==indx+1)[0]
plt.plot(np.nanmean(magn3D[binIndx,:],0).transpose(),np.nanmean(z[binIndx,:],0).transpose(),'k')
plt.plot(np.nanmin(magn3D[binIndx,:],0).transpose(),np.nanmean(z[binIndx,:],0).transpose(),'b')
plt.plot(np.nanmax(magn3D[binIndx,:],0).transpose(),np.nanmean(z[binIndx,:],0).transpose(),'b')
ax = plt.gca()
ax.set_position([0.07,0.1,0.9,0.8])
plt.xlabel('Depth below MSL [ m ]')
plt.ylabel('Speed [ m/s ]')
return fig
#--------------------------------------------------------------------------
# MAIN PROCESS
#--------------------------------------------------------------------------
#--------------------------------------------------------------------------
# INTERFACE
#--------------------------------------------------------------------------