Source code for ocmw.dataproc.ocmw_extract
# -*- coding: utf-8 -*-
"""
Functions to extract timeseries of data at specific locations in the water column from model data location extracts.
Chris Old
IES, School of Engineering, University of Edinburgh
Aug 2022
"""
# ----------------------------------------------------------------------------
# IMPORTS
# ----------------------------------------------------------------------------
# Standard Python Dependencies
import os
# Non-Standard Python Dependencies
import numpy as np
import scipy.io as sio
# Local Module Dependencies
from ocmw.core import flowChar as fc
from ocmw.core.fileTools import loadmat
from ocmw.core.graphics import plot_speedbinned_velocities
# Other Dependencies
# ----------------------------------------------------------------------------
# GLOBAL VARIABLES
# ----------------------------------------------------------------------------
ocmwExtn = '.mat'
# ----------------------------------------------------------------------------
# CLASS DEFINITIONS
# ----------------------------------------------------------------------------
# ----------------------------------------------------------------------------
# FUNCTION DEFINITIONS
# ----------------------------------------------------------------------------
# General form of extraction functions assuming OCMW format input data
[docs]
def extract_ts_at_hub(datapath,filename,dataset,datasrc,zhub,hab=True):
"""
Calculate TS of velocity data at depth zhub from OTM extracted data
Parameters
----------
datapath : str
Path to data file to be processed.
filename : str
Filename of data file to be processed.
dataset : str
Short description of data set being processed.
datasrc : str
Short descrition of how the data set was generated.
zhub : float
Vertical height where time series is to be extracted.
hab : bool, optional
Toggle to indicate whether zhub is height above bed or depth below
surface. The default is True.
Returns
-------
t_series : dict
Dictionary structure containing the time series of velocities at zhub.
"""
# Initialise output data structure
t_series = {}
t_series['dataSet'] = dataset
t_series['dataSrc'] = datasrc
vel = loadmat(datapath+'/'+filename)
dataloc = vel['dataLoc']
crs = vel['crs']
epoch = vel['epoch']
ts = vel['times']
dpth = vel['depth']
u = vel['vel_east']
v = vel['vel_north']
w = vel['vel_up']
z = vel['z']
if u.shape[0] != ts.shape[0]:
u = u.transpose()
v = v.transpose()
w = w.transpose()
z = z.transpose()
nrecs = len(u)
uhab = np.empty((nrecs,))
uhab.fill(np.nan)
vhab = np.empty((nrecs,))
vhab.fill(np.nan)
whab = np.empty((nrecs,))
whab.fill(np.nan)
ugrd = np.empty((nrecs,))
ugrd.fill(np.nan)
vgrd = np.empty((nrecs,))
vgrd.fill(np.nan)
wgrd = np.empty((nrecs,))
wgrd.fill(np.nan)
for irec in np.arange(nrecs):
ht = z[irec,:]-z[irec,0]
if not hab:
zh = ht[-1]-zhub
else:
zh = zhub
dudz = np.gradient(u[irec,:],ht)
dvdz = np.gradient(v[irec,:],ht)
dwdz = np.gradient(w[irec,:],ht)
uhab[irec] = np.interp(zh,ht,u[irec,:])
vhab[irec] = np.interp(zh,ht,v[irec,:])
whab[irec] = np.interp(zh,ht,w[irec,:])
ugrd[irec] = np.interp(zh,ht,dudz)
vgrd[irec] = np.interp(zh,ht,dvdz)
wgrd[irec] = np.interp(zh,ht,dwdz)
t_series['dataLoc'] = dataloc
t_series['crs'] = crs
t_series['hab'] = hab
t_series['zLoc'] = zhub
if hab:
t_series['mooring'] = 'bottom'
else:
t_series['mooring'] = 'floating'
t_series['times'] = ts
t_series['epoch'] = epoch
t_series['depth'] = dpth
t_series['vel_east'] = uhab
t_series['vel_north'] = vhab
t_series['vel_up'] = whab
t_series['dudz'] = ugrd
t_series['dvdz'] = vgrd
t_series['dwdz'] = wgrd
return t_series
[docs]
def extract_depthavg_ts(datapath,filename,dataset,datasrc):
"""
Extract TS depth averaged velocity from OTM extracted data
Parameters
----------
datapath : str
Path to data file to be processed.
filename : str
Filename of data file to be processed.
dataset : str
Short description of data set being processed.
datasrc : str
Short descrition of how the data set was generated.
Returns
-------
t_series : dict
Dictionary structure containing the time series of depth-averagd velocities.
"""
# Initialise output data structure
t_series = {}
t_series['dataSet'] = dataset
t_series['dataSrc'] = datasrc
vel = loadmat(datapath+'/'+filename)
dataloc = vel['dataLoc']
crs = vel['crs']
epoch = vel['epoch']
ts = vel['times']
dpth = vel['depth']
udav = vel['depAvg_U']
vdav = vel['depAvg_V']
t_series['dataLoc'] = dataloc
t_series['crs'] = crs
t_series['hab'] = True
t_series['zLoc'] = np.abs(np.nanmean(dpth))
t_series['mooring'] = 'bottom'
t_series['times'] = ts
t_series['epoch'] = epoch
t_series['depth'] = dpth
t_series['vel_east'] = udav
t_series['vel_north'] = vdav
t_series['vel_up'] = np.zeros(udav.shape)
return t_series
[docs]
def extract_rotav_ts_at_hub(datapath,filename,dataset,datasrc,ppdata,rdia,zhub,hab=True,pwra=True,Cp=0.33):
"""
Calculate TS of rotor averaged velocity at depth zhub from model extracted data
Parameters
----------
datapath : str
Path to data file to be processed.
filename : str
Filename of data file to be processed.
dataset : str
Short description of data set being processed.
datasrc : str
Short descrition of how the data set was generated.
ppdata : dict
Pre-processed time series data generated using "flowChar.preprocess_ts".
rdia : float
rotor diameter.
zhub : float
Vertical height where time series is to be extracted.
hab : bool, optional
Toggle to indicate whether zhub is height above bed or depth below
surface. The default is True.
pwra : bool, optional
Toggle to switch between power weighted rotor average and depth average. The default is True.
Cp : float, optional
Turbine power coefficient. The default is 0.33.
Returns
-------
t_series : TYPE
Dictionary structure containing the time series of rotor averaged velocities at zhub.
"""
# Set rotor z range
if hab:
zmin = zhub - rdia/2.0
zmax = zhub + rdia/2.0
else:
zmin = zhub + rdia/2.0
zmax = zhub - rdia/2.0
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()
# Initialise output data structure
t_series = {}
t_series['dataSet'] = dataset
t_series['dataSrc'] = datasrc
vel = loadmat(datapath+'/'+filename)
dataloc = vel['dataLoc']
crs = vel['crs']
epoch = vel['epoch']
ts = vel['times'].copy()
dpth = vel['depth'].copy()
u = vel['vel_east'].copy()
v = vel['vel_north'].copy()
z = vel['z'].copy()
if u.shape[0] != ts.shape[0]:
u = u.transpose()
v = v.transpose()
z = z.transpose()
nrecs = u.shape[0]
nlays = u.shape[1]
vel_strm = np.empty(u.shape)
vel_strm.fill(np.nan)
print(' Constructing streamwise velocity data...')
for ilay in np.arange(nlays):
vel_lay = np.vstack((u[ebbi,ilay],v[ebbi,ilay]))
vel_strm[ebbi,ilay] = np.sum(vel_lay*ebbu,0)
del vel_lay
vel_lay = np.vstack((u[fldi,ilay],v[fldi,ilay]))
vel_strm[fldi,ilay] = np.sum(vel_lay*fldu,0)
print(' Calculating rotor average velocity..')
vel_rot = np.empty((nrecs,))
for irec in np.arange(nrecs):
ht = z[irec,:]-z[irec,0]
if not hab:
zminp = ht[-1]-zmin
zmaxp = ht[-1]-zmax
else:
zminp = zmin
zmaxp = zmax
if pwra:
vel_rot[irec] = fc.pwra_vel(ht,vel_strm[irec,:],zminp,zmaxp,11)
else:
vel_rot[irec] = fc.rotor_average_vel(ht,vel_strm[irec,:],zminp,zmaxp,11)
tidal_pwr = fc.tec_power(vel_rot,rdia,Cp=Cp)
print(' Collating results into dictionary structure...')
t_series['dataLoc'] = dataloc
t_series['crs'] = crs
t_series['hab'] = hab
t_series['zLoc'] = zhub
if hab:
t_series['mooring'] = 'bottom'
else:
t_series['mooring'] = 'floating'
t_series['rdia'] = rdia
t_series['ebb_dir'] = ebb_angle
t_series['ebb_indices'] = ebbi
t_series['flood_dir'] = fld_angle
t_series['flood_indices'] = fldi
t_series['times'] = ts
t_series['epoch'] = epoch
t_series['depth'] = dpth
t_series['magn'] = vel_rot
t_series['tec_pwr'] = tidal_pwr
return t_series
[docs]
def extract_model_ts_at_surface(datapath,filename,dataset,datasrc):
"""
Extract TS of surface velocity from model extracted data
datapath : str
Path to data file to be processed.
filename : str
Filename of data file to be processed.
dataset : str
Short description of data set being processed.
datasrc : str
Short descrition of how the data set was generated.
Returns
-------
t_series : dict
Dictionary structure containing the time series of surface velocities.
"""
# Set extrac to 0.0m depth below surfac
zhub = 0.0
hab = False
# Load time series
vel = loadmat(datapath+'/'+filename)
# Create output data structure
t_series = {}
t_series['dataSet'] = dataset
t_series['dataSrc'] = datasrc
dataloc = vel['dataLoc']
crs = vel['crs']
epoch = vel['epoch']
ts = vel['times']
dpth = vel['depth']
if vel['vel_east'].shape[0] != ts.shape[0]:
usrf = vel['vel_east'][-1,:]
vsrf = vel['vel_north'][-1,:]
wsrf = vel['vel_up'][-1,:]
else:
usrf = vel['vel_east'][:,-1]
vsrf = vel['vel_north'][:,-1]
wsrf = vel['vel_up'][:,-1]
t_series['dataLoc'] = dataloc
t_series['crs'] = crs
t_series['hab'] = hab
t_series['zLoc'] = zhub
if hab:
t_series['mooring'] = 'bottom'
else:
t_series['mooring'] = 'floating'
t_series['times'] = ts
t_series['epoch'] = epoch
t_series['depth'] = dpth
t_series['vel_east'] = usrf
t_series['vel_north'] = vsrf
t_series['vel_up'] = wsrf
return t_series
[docs]
def checkOCMWExtn(fileName):
fpath,fname = os.path.split(fileName)
fn, fe = os.path.splitext(fname)
if fe not in [ocmwExtn]:
fname = fn + ocmwExtn
return fpath,fname
[docs]
def load_ocmw_file(dataPath,fileName):
"""
Reader for loading the OCMW formated extract files
Parameters
----------
datapath : str
Path to data file to be processed.
filename : str
Filename of data file to be processed.
Returns
-------
data : dict
Dictionary structure containing the data from the OCMW file.
"""
data = loadmat(dataPath+'/'+fileName)
return data
[docs]
def save_ocmw_file(dataPath,ocmwFileName,dataDict):
"""
Write data dictionary to OCMW formated output file.
Parameters
----------
dataPath : str
Path to where data are to be saved.
ocmwFileName : str
Filename ofoutput OCMW data file.
dataDict: dict
Dictionary containing data to be written to OCMW output file
Returns
-------
none
"""
ocmwFile = dataPath+'/'+ocmwFileName
sio.savemat(ocmwFile,dataDict)
return
[docs]
def speedBinIndices(dataPath,ocmwPPFile,binWidth):
"""
Calculate speed bin indices based on OCMW timeseries extracted at a
fixed height above the seabed.
"""
pp = load_ocmw_file(dataPath,ocmwPPFile)
m = pp['magn']
nbins = (np.max(m)//binWidth)+1
bins = np.arange(0.0,(nbins+1)*binWidth,binWidth)
indices = np.searchsorted(bins,m)
binCntrs = (bins[1:]+bins[:-1])/2.0
return indices, binCntrs
[docs]
def speedBinning(rootPath,dataFile,ocmwPPFile,binWidth):
"""
Speed bin OCMW extracted velocity data using a pre-processed OCMW data
file to define the speed bin indices.
"""
rdPath = rootPath+'/EXTRACT'
ppPath = rootPath+'/ANALYSIS'
indices, cntrs = speedBinIndices(ppPath,ocmwPPFile,binWidth)
rd = load_ocmw_file(rdPath,dataFile)
magn3D = np.sqrt(np.square(rd['vel_east'])+np.square(rd['vel_north']))
z = rd['z']
fig = plot_speedbinned_velocities(magn3D, z, indices, cntrs)
return fig
# ----------------------------------------------------------------------------
# MAIN PROCESS
# ----------------------------------------------------------------------------
#--------------------------------------------------------------------------
# INTERFACE
#--------------------------------------------------------------------------