# -*- coding: utf-8 -*-
"""
Preprocessor for extracting flow characterisaton information from a timeseries of tidal flow data.
"""
#Flow Characterisation tools
#
#Chris Old
#IES, School of Engineering, University of Edinburgh
#Aug 2022
# ----------------------------------------------------------------------------
# IMPORTS
# ----------------------------------------------------------------------------
# Standard Python Dependencies
import numpy as np
# Non-Standard Python Dependencies
from sklearn.cluster import KMeans
from scipy.signal import convolve
# Local Module Dependencies
from .timeFuncs import dateNumFromDateStr, dateStrFromDateNum
# Other Dependencies
# ----------------------------------------------------------------------------
# GLOBAL VARIABLES
# ----------------------------------------------------------------------------
# ----------------------------------------------------------------------------
# CLASS DEFINITIONS
# ----------------------------------------------------------------------------
# ----------------------------------------------------------------------------
# FUNCTION DEFINITIONS
# ----------------------------------------------------------------------------
[docs]
def moving_average(x, w):
"""Apply moving average to a series of data with a defined sample window width using convolution.
"""
return convolve(x, np.ones(w), method='direct') / w
[docs]
def movmean(v,n):
"""Apply moving mean to a series of data with a defined sample window width using running mean.
"""
ret = np.nancumsum(v, dtype=float)
ret[n:] = ret[n:] - ret[:-n]
mavg = ret[n:] / n
return mavg
[docs]
def chord_area(R,h):
"""Calculate the area of a chord for a circle of radius R and a chord height h from the origin.
"""
F1 = 0.5*np.pi*R*R
F2 = R*R*np.arcsin(h/R)
F3 = h * np.sqrt(R*R - h*h)
Ac = F1-F2-F3
return Ac
[docs]
def pwra_vel(z,v,zmin,zmax,nz=3):
"""Calculate the power-weighted rotor avearge velocity.
"""
nsegs = nz-1
dz = (zmax-zmin)/(nsegs)
zr = np.arange(nz)*dz
vz = np.interp(zr+zmin,z,v)
dv = (vz[0:-1]+vz[1:])/2.0
r = (zmax-zmin)/2.0
vr = 0.0
A = 0.0
for seg in np.arange(nsegs):
h1 = zr[seg]-r
h2 = zr[seg+1]-r
Ac1 = chord_area(r,h1)
Ac2 = chord_area(r,h2)
dA = np.abs(Ac2-Ac1)
vr = vr + dA * np.power(dv[seg],3)
A = A + dA
pwra = np.power(np.abs(vr)/A,1.0/3.0)
return pwra
[docs]
def rotor_average_vel(z,v,zmin,zmax,nz=3):
"""Calculate the rotor averaged velocity.
"""
dz = (zmax-zmin)/(nz-1)
zr = zmin + np.arange(nz)*dz
vr = np.interp(zr,z,v)
rotav = np.nanmean(vr)
return rotav
[docs]
def tec_power(vel_rot,rdia,Cp=0.33,rho_sw=1025.0):
"""Calculate the TEC power based on an a timeseries of representative hu-height velocity.
"""
A = np.pi*np.square(rdia/2.0)
Pwr = 0.5*rho_sw*A*np.power(vel_rot,3)*Cp
return Pwr
[docs]
def streamwise_flow(ppdata):
"""Calculate the flood and ebb streamwise velocities from pre-processed data.
"""
u = ppdata['vel_east']
v = ppdata['vel_north']
ebbi = ppdata['ebb_indices']
ebb_angle = ppdata['peak_dir']['dir'][ppdata['peak_id']['ebb']]
ebb_uv = [np.cos(np.radians(ebb_angle)),np.sin(np.radians(ebb_angle))]
ebbu = np.tile(ebb_uv,(len(ebbi),1)).transpose()
fldi = ppdata['flood_indices']
fld_angle = ppdata['peak_dir']['dir'][ppdata['peak_id']['flood']]
fld_uv = [np.cos(np.radians(fld_angle)),np.sin(np.radians(fld_angle))]
fldu = np.tile(fld_uv,(len(fldi),1)).transpose()
vel_strm = np.empty(u.shape)
vel_strm.fill(np.nan)
print(' Constructing streamwise velocity data...')
vel_lay = np.vstack((u[ebbi],v[ebbi]))
vel_strm[ebbi] = np.sum(vel_lay*ebbu,0)
del vel_lay
vel_lay = np.vstack((u[fldi],v[fldi]))
vel_strm[fldi] = np.sum(vel_lay*fldu,0)
return vel_strm
[docs]
def flow_acceleration(data):
"""Calculate the flow acceleration from a timeseries of velocity data exttracted at a specific depth.
"""
t = data['times']
u = data['vel_east']
v = data['vel_north']
dt = np.median(np.diff(t))*24.0*60.0*60.0
dudt = np.gradient(u)/dt
dvdt = np.gradient(v)/dt
magn = np.sqrt(np.power(dudt,2) + np.power(dvdt,2))
dirn = np.rad2deg(np.arctan2(dvdt,dudt))
flw_acc = {}
flw_acc['times'] = t
flw_acc['dudt'] = dudt
flw_acc['dvdt'] = dvdt
flw_acc['magn'] = magn
flw_acc['dirn'] = dirn
return flw_acc
[docs]
def get_flood_ebb_indices(thta,pkDir,fld_id):
"""Get the ebb/flood indices from a timeseries of flow direction at a specific depth.
"""
mean_dir = pkDir['mean_dir']
fld_peak_dir = pkDir['dir'][fld_id]
# Flood and Ebb angle ranges
if (fld_peak_dir - mean_dir) < 0:
fldangles = [mean_dir-180.0,mean_dir]
if fldangles[0] < -180.0:
fldangles = [[-180.0,mean_dir],[fldangles[0]+360.0, 180.0]]
ebbangles = [mean_dir,mean_dir+180.0]
if ebbangles[1] > 180.0:
ebbangles = [[mean_dir,180.0],[-180.0,ebbangles[1]-360.0]]
else:
fldangles = [mean_dir,mean_dir+180.0]
if fldangles[1] > 180.0:
fldangles = [[mean_dir,180.0],[-180.0,fldangles[1]-360.0]]
ebbangles = [mean_dir-180.0,mean_dir]
if ebbangles[0] < -180.0:
ebbangles = [[-180.0,mean_dir],[ebbangles[0]+360.0, 180.0]]
# Get indices
if len(np.asarray(fldangles).shape) == 1:
fldind = np.where((thta > fldangles[0]) & (thta < fldangles[1]))[0]
else:
fldind = np.where(((thta > fldangles[0][0]) & (thta < fldangles[0][1])) | ((thta > fldangles[1][0]) & (thta < fldangles[1][1])))[0]
flood_indices = fldind
if len(np.asarray(ebbangles).shape) == 1:
ebbind = np.where((thta > ebbangles[0]) & (thta < ebbangles[1]))[0]
else:
ebbind = np.where(((thta > ebbangles[0][0]) & (thta < ebbangles[0][1])) | ((thta > ebbangles[1][0]) & (thta < ebbangles[1][1])))[0]
ebb_indices = ebbind
return [flood_indices, ebb_indices]
[docs]
def id_tidal_flow_peaks(thta,elv,pkDir):
"""Identify peak directions as ebb and flood.
"""
peak_id = {}
gdindx = np.where(~np.isnan(elv))[0]
fidx = np.where(elv[gdindx] > 0.25*np.nanmax(elv))[0]
if fidx.size == 0:
peak_id['flood'] = 0
peak_id['ebb'] = 1
else:
fldang = np.nanmedian(thta[fidx])
delang = np.abs(pkDir['dir']-fldang)
peak_id['flood'] = np.where(delang == np.min(delang))[0][0]
if peak_id['flood'] == 0:
peak_id['ebb'] = 1
else:
peak_id['ebb'] = 0
return peak_id
[docs]
def tidal_flow_peak_directions(thta,magn):
"""Determine peak flow directions based on a timeseries of flow magnitude and direction.
"""
gdindx = np.where(~np.isnan(magn))[0]
kmeans = KMeans(n_clusters=2,random_state=42,n_init="auto")
X = np.vstack((thta[gdindx],magn[gdindx])).transpose()
kmeans.fit(X)
km = kmeans.labels_
pk1 = np.where(km == 0)[0]
spdrng = np.where(magn[gdindx[pk1]] > 0.5*np.nanmax(magn))[0]
if spdrng.size == 0:
spdrng = np.where(magn[gdindx[pk1]] > 0.25*np.nanmax(magn))[0]
if spdrng.size == 0:
spdrng = np.where(magn[gdindx[pk1]] >= 0.0)[0]
pk1dir = np.median(thta[gdindx[pk1[spdrng]]])
if pk1dir > 0:
pk1bearing = 90.0 - pk1dir
else:
pk1bearing = 90.0 - (360.0 + pk1dir)
if pk1bearing < 0:
pk1bearing = 360.0 + pk1bearing
pk2 = np.where(km == 1)[0]
spdrng = np.where(magn[gdindx[pk2]] > 0.5*np.nanmax(magn))[0]
if spdrng.size == 0:
spdrng = np.where(magn[gdindx[pk2]] > 0.25*np.nanmax(magn))[0]
if spdrng.size == 0:
spdrng = np.where(magn[gdindx[pk2]] >= 0.0)[0]
pk2dir = np.median(thta[gdindx[pk2[spdrng]]])
if pk2dir > 0:
pk2bearing = 90.0 - pk2dir
else:
pk2bearing = 90.0 - (360.0 + pk2dir)
if pk2bearing < 0:
pk2bearing = 360.0 + pk2bearing
deldir = np.abs(pk2dir-pk1dir)
if deldir > 180.0:
deldir = 360.0-deldir
meandir = (pk1dir+pk2dir)/2.0
pkDir = {}
pkDir['dir'] = [pk1dir,pk2dir]
pkDir['bearing'] = [pk1bearing,pk2bearing]
pkDir['del_dir'] = deldir
pkDir['mean_dir'] = meandir
return pkDir
[docs]
def getFlowPercentage(magn,cutInSpeed,ratedSpeed,indices):
"""
Calculate percentage of time in flow speed ranges (below cut-in, operating, above rated)
"""
pctage = [100.0*len(np.where(magn[indices] < cutInSpeed)[0])/len(indices), \
100.0*len(np.where((magn[indices] >= cutInSpeed) & (magn[indices] <= ratedSpeed))[0])/len(indices), \
100.0*len(np.where(magn[indices] > ratedSpeed)[0])/len(indices)]
return pctage
[docs]
def preprocess_ts(tseries,tec=None,ndays=None,tstart=None):
"""
Preprocess a timeseries of tidal data for a single location and specific depth in the water column.
"""
# Convert velocity to magnitude and direction
t = tseries['times'].copy()
epoch = tseries['epoch']
#epoch_off = dateNumFromStr(epoch.replace('/','-'))
epoch_off = dateNumFromDateStr(epoch.replace('/','-'))
if ndays is None:
# third parameter does not exist, so default it to something
#ndays = np.floor(t[-1]-t[0])+1
ndays = np.floor(t[-2]-t[0])+1
start_date = dateStrFromDateNum(np.floor(t[0]+epoch_off))
end_date = dateStrFromDateNum(t[-1]+epoch_off)
npts = len(t)
indx1 = 0
indx2 = npts
else:
dt = np.nanmedian(np.diff(t))
ndaysoffset = np.int64(np.ceil(ndays/dt))
if tstart is None:
start_date = dateStrFromDateNum(np.floor(t[0]+epoch_off))
indx1 = 1
end_date = dateStrFromDateNum(np.floor(t[ndaysoffset]+epoch_off))
else:
start_date = dateStrFromDateNum(tstart)
indx_t = np.where(t >= tstart-epoch_off)[0]
indx1 = indx_t[0]
end_date = dateStrFromDateNum(np.floor(t[indx1+ndaysoffset]+epoch_off))
npts = ndaysoffset
indx2 = indx1 + ndaysoffset
deploystr = '( Data covers '+str(ndays)+' days from '+start_date+' )'
if 'depth' in tseries.keys():
dth = tseries['depth'][indx1:indx2].copy()
else:
dth = tseries['pressure'][indx1:indx2].copy()
if np.sum(np.diff(dth)) == 0:
dth[:] = np.nan
msl = np.nanmedian(dth)
elv = dth - msl
baddth = np.where(np.abs(elv) > 1.5*np.nanmax(elv))[0]
dth[baddth] = np.nan
msl = np.nanmedian(dth)
elv = dth - msl
u = tseries['vel_east'][indx1:indx2]
v = tseries['vel_north'][indx1:indx2]
w = tseries['vel_up'][indx1:indx2]
u[baddth] = np.nan
v[baddth] = np.nan
w[baddth] = np.nan
magn = np.sqrt(np.power(u[0:npts],2) + np.power(v[0:npts],2))
thta = np.rad2deg(np.arctan2(v[0:npts],u[0:npts]))
res_magn = np.sqrt(np.power(np.nanmean(u[0:npts]),2) + np.power(np.nanmean(v[0:npts]),2))
res_thta = np.rad2deg(np.arctan2(np.nanmean(v[0:npts]),np.nanmean(u[0:npts])))
if res_thta < -180.0:
res_thta = res_thta + 360.0
elif res_thta > 180.0:
res_thta = res_thta - 360.0
# Find tidal flow direction peaks
pkdir = tidal_flow_peak_directions(thta,magn)
# Identify flood peak
peak_id = id_tidal_flow_peaks(thta,elv,pkdir)
# Get Flood and Ebb Indices
[fldindices, ebbindices] = get_flood_ebb_indices(thta,pkdir,peak_id['flood'])
# Calculate percentage of time in flow speed ranges (below cut-in, operating, above rated)
if tec is not None:
fldpctage = getFlowPercentage(magn,tec.cutInSpeed,tec.ratedSpeed,fldindices)
ebbpctage = getFlowPercentage(magn,tec.cutInSpeed,tec.ratedSpeed,ebbindices)
else:
fldpctage = [np.nan, np.nan, np.nan]
ebbpctage = [np.nan, np.nan, np.nan]
ngood = len(fldindices) + len(ebbindices)
fldpc = 100.0*len(fldindices)/ngood
ebbpc = 100.0*len(ebbindices)/ngood
# Calculate percentage of time flow in direction bins around peak
flddir = pkdir['dir'][peak_id['flood']]
ebbdir = pkdir['dir'][peak_id['ebb']]
dfldang = thta[fldindices]-flddir
dfldang[dfldang < -180.0] = dfldang[dfldang < -180.0] + 360.0
dfldang[dfldang > 180.0] = dfldang[dfldang > 180.0] - 360.0
debbang = thta[ebbindices]-ebbdir
debbang[debbang < -180.0] = debbang[debbang < -180.0] + 360.0
debbang[debbang > 180.0] = debbang[debbang > 180.0] - 360.0
# Collect Information
data = {}
data['dataSet'] = tseries['dataSet']
data['dataSrc'] = tseries['dataSrc']
data['deploystr'] = deploystr
data['dataLoc'] = tseries['dataLoc']
data['crs'] = tseries['crs']
data['hab'] = tseries['hab']
data['zLoc'] = tseries['zLoc']
data['times'] = t[indx1:indx2]
data['epoch'] = epoch
data['number_of_days'] = ndays
data['start_date'] = start_date
data['end_date'] = end_date
data['vel_east'] = u
data['vel_north'] = v
data['vel_up'] = w
data['thta'] = thta
data['magn'] = magn
data['res_thta'] = res_thta
data['res_magn'] = res_magn
data['peak_dir'] = pkdir
data['peak_id'] = peak_id
data['flood_indices'] = fldindices
data['ebb_indices'] = ebbindices
data['del_flood_dir'] = dfldang
data['del_ebb_dir'] = debbang
data['flood_total_percent'] = fldpc
data['flood_percent'] = fldpctage
data['ebb_total_percent'] = ebbpc
data['ebb_percent'] = ebbpctage
data['depth'] = dth
data['msl'] = msl
data['surface_elev'] = elv
if 'DataFile' in tseries.keys():
data['datafile'] = tseries['DataFile'].copy()
data['data_index_1'] = tseries['DataIndex1']
data['data_index_2'] = tseries['DataIndex2']
# Calculate flow acceleration
flw_acc = flow_acceleration(data)
# Add flow acceleration
data['flw_acc'] = flw_acc
return data
# ----------------------------------------------------------------------------
# MAIN PROCESS
# ----------------------------------------------------------------------------
#--------------------------------------------------------------------------
# INTERFACE
#--------------------------------------------------------------------------