Binning Tools

Velocity Analysis

Analysis in DOLfYN is primarily handled through the VelBinner class. Below is a list of functions that can be called from VelBinner.

`VelBinner`	This is the base binning (averaging) tool.
`reshape`	Reshape the array arr to shape (...,n,n_bin+n_pad).
`detrend`	Reshape the array arr to shape (...,n,n_bin+n_pad) and remove the best-fit trend line from each bin.
`demean`	Reshape the array arr to shape (...,n,n_bin+n_pad) and remove the mean from each bin.
`mean`	Reshape the array arr to shape (...,n,n_bin+n_pad) and take the mean of each bin along the specified axis.
`var`	Reshape the array arr to shape (...,n,n_bin+n_pad) and take the variance of each bin along the specified axis.
`std`	Reshape the array arr to shape (...,n,n_bin+n_pad) and take the standard deviation of each bin along the specified axis.
`do_avg`	Bin the dataset and calculate the ensemble averages of each variable.
`do_var`	Bin the dataset and calculate the ensemble variances of each variable.
`calc_coh`	Calculate coherence between veldat1 and veldat2.
`calc_phase_angle`	Calculate the phase difference between two signals as a function of frequency (complimentary to coherence).
`calc_acov`	Calculate the auto-covariance of the raw-signal veldat
`calc_xcov`	Calculate the cross-covariance between arrays veldat1 and veldat2
`calc_tke`	Calculate the turbulent kinetic energy (TKE) (variances of u,v,w).
`calc_psd`	Calculate the power spectral density of velocity.
`calc_freq`	Calculate the ordinary or radial frequency vector for the PSDs
`calc_psd_base`	Calculate power spectral density of dat
`calc_csd_base`	Calculate the cross power spectral density of dat.

Turbulence Analysis

Functions for analyzing ADV data via the ADVBinner class, beyond those described in VelBinner.

`ADVBinner`	A class that builds upon VelBinner for calculating turbulence statistics and velocity spectra from ADV data
`calc_turbulence`	Functional version of ADVBinner that computes a suite of turbulence statistics for the input dataset, and returns a binned data object.
`calc_csd`	Calculate the cross-spectral density of velocity components.
`calc_stress`	Calculate the stresses (covariances of u,v,w)
`calc_doppler_noise`	Calculate bias due to Doppler noise using the noise floor of the velocity spectra.
`check_turbulence_cascade_slope`	This function calculates the slope of the PSD, the power spectra of velocity, within the given frequency range.
`calc_epsilon_LT83`	Calculate the dissipation rate from the PSD
`calc_epsilon_SF`	Calculate dissipation rate using the "structure function" (SF) method
`calc_epsilon_TE01`	Calculate the dissipation rate according to TE01.
`calc_L_int`	Calculate integral length scales.

Functions for analyzing ADCP data via the ADPBinner class, beyond those described in VelBinner.

`ADPBinner`	A class for calculating turbulence statistics from ADCP data
`calc_dudz`	The shear in the first velocity component.
`calc_dvdz`	The shear in the second velocity component.
`calc_dwdz`	The shear in the third velocity component.
`calc_shear2`	The horizontal shear squared.
`calc_doppler_noise`	Calculate bias due to Doppler noise using the noise floor of the velocity spectra.
`calc_stress_4beam`	Calculate the stresses from the covariance of along-beam velocity measurements
`calc_stress_5beam`	Calculate the stresses from the covariance of along-beam velocity measurements
`calc_total_tke`	Calculate magnitude of turbulent kinetic energy from 5-beam ADCP.
`check_turbulence_cascade_slope`	This function calculates the slope of the PSD, the power spectra of velocity, within the given frequency range.
`calc_dissipation_LT83`	Calculate the TKE dissipation rate from the velocity spectra.
`calc_dissipation_SF`	Calculate TKE dissipation rate from ADCP along-beam velocity using the "structure function" (SF) method.
`calc_ustar_fit`	Approximate friction velocity from shear stress using a logarithmic profile.

class dolfyn.binned.TimeBinner(n_bin, fs, n_fft=None, n_fft_coh=None, noise=[0, 0, 0])[source]

Bases: object

Initialize an averaging object

Parameters:

n_bin (int) – Number of data points to include in a ‘bin’ (ensemble), not the number of bins
fs (int) – Instrument sampling frequency in Hz
n_fft (int) – Number of data points to use for fft (n_fft`<=`n_bin). Default: n_fft`=`n_bin
n_fft_coh (int) – Number of data points to use for coherence and cross-spectra ffts Default: n_fft_coh`=`n_fft
noise (float, list or numpy.ndarray) – Instrument’s doppler noise in same units as velocity

reshape(arr, n_pad=0, n_bin=None)[source]

Reshape the array arr to shape (…,n,n_bin+n_pad).

Parameters:

arr (numpy.ndarray) – Input data
n_pad (int) – Is used to add n_pad/2 points from the end of the previous ensemble to the top of the current, and n_pad/2 points from the top of the next ensemble to the bottom of the current. Zeros are padded in the upper-left and lower-right corners of the matrix (beginning/end of timeseries). In this case, the array shape will be (…,`n`,`n_pad`+`n_bin`)
n_bin (int (default: self.n_bin)) – Override this binner’s n_bin.

Returns:

out (numpy.ndarray) – Data in reshaped format, where the last axis is of length int(n_bin).

Notes

n_bin can be non-integer, in which case the output array size will be n_pad`+`n_bin, and the decimal will cause skipping of some data points in arr. In particular, every mod(n_bin,1) bins will have a skipped point. For example: - for n_bin=2048.2 every 1/5 bins will have a skipped point. - for n_bin=4096.9 every 9/10 bins will have a skipped point.

detrend(arr, axis=-1, n_pad=0, n_bin=None)[source]

Reshape the array arr to shape (…,n,n_bin+n_pad) and remove the best-fit trend line from each bin.

Parameters:

arr (numpy.ndarray) – Input data
axis (int (default: -1)) – Axis along which to take mean
n_pad (int (default: 0)) – Is used to add n_pad/2 points from the end of the previous ensemble to the top of the current, and n_pad/2 points from the top of the next ensemble to the bottom of the current. Zeros are padded in the upper-left and lower-right corners of the matrix (beginning/end of timeseries). In this case, the array shape will be (…,`n`,`n_pad`+`n_bin`)
n_bin (int (default: self.n_bin)) – Override this binner’s n_bin.

Returns:

out (numpy.ndarray) – Detrended data, where the last axis is of length int(n_bin).

demean(arr, axis=-1, n_pad=0, n_bin=None)[source]

Reshape the array arr to shape (…,n,n_bin+n_pad) and remove the mean from each bin.

Parameters:

arr (numpy.ndarray) – Input data
axis (int (default: -1)) – Axis along which to take mean
n_pad (int (default: 0)) – Is used to add n_pad/2 points from the end of the previous ensemble to the top of the current, and n_pad/2 points from the top of the next ensemble to the bottom of the current. Zeros are padded in the upper-left and lower-right corners of the matrix (beginning/end of timeseries). In this case, the array shape will be (…,`n`,`n_pad`+`n_bin`)
n_bin (int (default: self.n_bin)) – Override this binner’s n_bin.

Returns:

out (numpy.ndarray) – Demeaned data, where the last axis is of length int(n_bin).

mean(arr, axis=-1, n_bin=None)[source]

Reshape the array arr to shape (…,n,n_bin+n_pad) and take the mean of each bin along the specified axis.

Parameters:

arr (numpy.ndarray) – Input data
axis (int (default: -1)) – Axis along which to take mean
n_bin (int (default: self.n_bin)) – Override this binner’s n_bin.

Returns:

out (numpy.ndarray)

var(arr, axis=-1, n_bin=None)[source]

Reshape the array arr to shape (…,n,n_bin+n_pad) and take the variance of each bin along the specified axis.

Parameters:

arr (numpy.ndarray) – Input data
axis (int (default: -1)) – Axis along which to take variance
n_bin (int (default: self.n_bin)) – Override this binner’s n_bin.

Returns:

out (numpy.ndarray)

std(arr, axis=-1, n_bin=None)[source]

Reshape the array arr to shape (…,n,n_bin+n_pad) and take the standard deviation of each bin along the specified axis.

Parameters:

arr (numpy.ndarray) – Input data
axis (int (default: -1)) – Axis along which to take std dev
n_bin (int (default: self.n_bin)) – Override this binner’s n_bin.

Returns:

out (numpy.ndarray)

calc_psd_base(dat, fs=None, window='hann', noise=0, n_bin=None, n_fft=None, n_pad=None, step=None)[source]

Calculate power spectral density of dat

Parameters:

dat (xarray.DataArray) – The raw dataArray of which to calculate the psd.
fs (float (optional)) – The sample rate (Hz).
window (str) – String indicating the window function to use (default: ‘hanning’).
noise (float) – The white-noise level of the measurement (in the same units as dat).
n_bin (int) – n_bin of veldat2, number of elements per bin if ‘None’ is taken from VelBinner
n_fft (int) – n_fft of veldat2, number of elements per bin if ‘None’ is taken from VelBinner
n_pad (int (optional)) – The number of values to pad with zero (default: 0)
step (int (optional)) – Controls amount of overlap in fft (default: the step size is chosen to maximize data use, minimize nens, and have a minimum of 50% overlap.).

Returns:

out (numpy.ndarray) – The power spectral density of dat

Notes

PSD’s are calculated based on sample rate units

calc_csd_base(dat1, dat2, fs=None, window='hann', n_fft=None, n_bin=None)[source]

Calculate the cross power spectral density of dat.

Parameters:

dat1 (numpy.ndarray) – The first (shorter, if applicable) raw dataArray of which to calculate the cpsd.
dat2 (numpy.ndarray) – The second (the shorter, if applicable) raw dataArray of which to calculate the cpsd.
fs (float (optional)) – The sample rate (rad/s or Hz).
window (str) – String indicating the window function to use (default: ‘hanning’).
n_fft (int) – n_fft of veldat2, number of elements per bin if ‘None’ is taken from VelBinner
n_bin (int) – n_bin of veldat2, number of elements per bin if ‘None’ is taken from VelBinner

Returns:

out (numpy.ndarray) – The cross power spectral density of dat1 and dat2

Notes

PSD’s are calculated based on sample rate units

The two velocity inputs do not have to be perfectly synchronized, but they should have the same start and end timestamps

calc_freq(fs=None, units='rad/s', n_fft=None, coh=False)[source]

Calculate the ordinary or radial frequency vector for the PSDs

Parameters:

fs (float (optional)) – The sample rate (Hz).
units (string) – Frequency units in either Hz or rad/s (f or omega)
coh (bool) – Calculate the frequency vector for coherence/cross-spectra (default: False) i.e. use self.n_fft_coh instead of self.n_fft.
n_fft (int) – n_fft of veldat2, number of elements per bin if ‘None’ is taken from VelBinner

Returns:

out (numpy.ndarray) – Spectrum frequency array in units of ‘Hz’ or ‘rad/s’

class dolfyn.velocity.Velocity(ds, *args, **kwargs)[source]

Bases: object

All ADCP and ADV xarray datasets wrap this base class. The turbulence-related attributes defined within this class assume that the 'tke_vec' and 'stress_vec' data entries are included in the dataset. These are typically calculated using a VelBinner tool, but the method for calculating these variables can depend on the details of the measurement (instrument, it’s configuration, orientation, etc.).

See also

VelBinner

rotate2(out_frame='earth', inplace=True)[source]

Rotate the dataset to a new coordinate system.

Parameters:

out_frame (string {'beam', 'inst', 'earth', 'principal'}) – The coordinate system to rotate the data into.
inplace (bool (default: True)) – When True the existing data object is modified. When False a copy is returned.

Returns:

ds (xarray.Dataset or None) – Returns the rotated dataset when ``inplace=False``, otherwise returns None.

Notes

This function rotates all variables in ds.attrs['rotate_vars'].
To rotate to the ‘principal’ frame, a value of ds.attrs['principal_heading'] must exist. The function calc_principal_heading is recommended for this purpose, e.g.:
```
ds.attrs['principal_heading'] = dolfyn.calc_principal_heading(ds['vel'].mean(range))
```
where here we are using the depth-averaged velocity to calculate the principal direction.

set_declination(declin, inplace=True)[source]

Set the magnetic declination

Parameters:

declination (float) – The value of the magnetic declination in degrees (positive values specify that Magnetic North is clockwise from True North)
inplace (bool (default: True)) – When True the existing data object is modified. When False a copy is returned.

Returns:

ds (xarray.Dataset or None) – Returns the rotated dataset when ``inplace=False``, otherwise returns None.

Notes

This method modifies the data object in the following ways:

If the dataset is in the earth reference frame at the time of

setting declination, it will be rotated into the “True-East, True-North, Up” (hereafter, ETU) coordinate system

dat['orientmat'] is modified to be an ETU to

instrument (XYZ) rotation matrix (rather than the magnetic-ENU to XYZ rotation matrix). Therefore, all rotations to/from the ‘earth’ frame will now be to/from this ETU coordinate system.

The value of the specified declination will be stored in

dat.attrs['declination']

dat['heading'] is adjusted for declination

(i.e., it is relative to True North).

If dat.attrs['principal_heading'] is set, it is

adjusted to account for the orientation of the new ‘True’ earth coordinate system (i.e., calling set_declination on a data object in the principal coordinate system, then calling dat.rotate2(‘earth’) will yield a data object in the new ‘True’ earth coordinate system)

set_inst2head_rotmat(rotmat, inplace=True)[source]

Set the instrument to head rotation matrix for the Nortek ADV if it hasn’t already been set through a ‘.userdata.json’ file.

Parameters:

rotmat (float) – 3x3 rotation matrix
inplace (bool (default: True)) – When True the existing data object is rotated. When False a copy is returned that is rotated.

Returns:

ds (xarray.Dataset or None) – Returns the rotated dataset when ``inplace=False``, otherwise returns None.

Notes

If the data object is in earth or principal coords, it is first rotated to ‘inst’ before assigning inst2head_rotmat, it is then rotated back to the coordinate system in which it was input. This way the inst2head_rotmat gets applied correctly (in inst coordinate system).

save(filename, **kwargs)[source]

Save the data object (underlying xarray dataset) as netCDF (.nc).

Parameters:

filename (str) – Filename and/or path with the ‘.nc’ extension
**kwargs (these are passed directly to xarray.Dataset.to_netcdf().) –

Notes

See DOLfYN’s save function for additional details.

property variables: A sorted list of the variable names in the dataset.

property attrs: The attributes in the dataset.

property coords: The coordinates in the dataset.

property u

The first velocity component.

This is simply a shortcut to self[‘vel’][0]. Therefore, depending on the coordinate system of the data object (self.attrs[‘coord_sys’]), it is:

beam: beam1
inst: x
earth: east
principal: streamwise

property v

The second velocity component.

This is simply a shortcut to self[‘vel’][1]. Therefore, depending on the coordinate system of the data object (self.attrs[‘coord_sys’]), it is:

beam: beam2
inst: y
earth: north
principal: cross-stream

property w: The third velocity component. This is simply a shortcut to self[‘vel’][2]. Therefore, depending on the coordinate system of the data object (self.attrs[‘coord_sys’]), it is: - beam: beam3 - inst: z - earth: up - principal: up

property U: Horizontal velocity as a complex quantity

property U_mag: Horizontal velocity magnitude

property U_dir: Angle of horizontal velocity vector. Direction is ‘to’, as opposed to ‘from’. This function calculates angle as “degrees CCW from X/East/streamwise” and then converts it to “degrees CW from X/North/streamwise”.

property E_coh

Coherent turbulence energy

Niel Kelley’s ‘coherent turbulence energy’, which is the RMS of the Reynold’s stresses.

See: NREL Technical Report TP-500-52353

property I_tke

Turbulent kinetic energy intensity.

Ratio of sqrt(tke) to horizontal velocity magnitude.

property I

Turbulence intensity.

Ratio of standard deviation of horizontal velocity to horizontal velocity magnitude.

property tke: Turbulent kinetic energy (sum of the three components)

property upvp_: u’v’bar Reynolds stress

property upwp_: u’w’bar Reynolds stress

property vpwp_: v’w’bar Reynolds stress

property upup_: u’u’bar component of the tke

property vpvp_: v’v’bar component of the tke

property wpwp_: w’w’bar component of the tke

class dolfyn.velocity.VelBinner(n_bin, fs, n_fft=None, n_fft_coh=None, noise=[0, 0, 0])[source]

Bases: TimeBinner

This is the base binning (averaging) tool. All DOLfYN binning tools derive from this base class.

Examples

The VelBinner class is used to compute averages and turbulence statistics from ‘raw’ (not averaged) ADV or ADP measurements, for example:

# First read or load some data.
rawdat = dolfyn.read_example('BenchFile01.ad2cp')

# Now initialize the averaging tool:
binner = dolfyn.VelBinner(n_bin=600, fs=rawdat.fs)

# This computes the basic averages
avg = binner.do_avg(rawdat)

Initialize an averaging object

Parameters:

n_bin (int) – Number of data points to include in a ‘bin’ (ensemble), not the number of bins
fs (int) – Instrument sampling frequency in Hz
n_fft (int) – Number of data points to use for fft (n_fft`<=`n_bin). Default: n_fft`=`n_bin
n_fft_coh (int) – Number of data points to use for coherence and cross-spectra ffts Default: n_fft_coh`=`n_fft
noise (float, list or numpy.ndarray) – Instrument’s doppler noise in same units as velocity

tke = <xarray.DataArray 'tke' (tke: 3)> array(['upup_', 'vpvp_', 'wpwp_'], dtype='<U5') Dimensions without coordinates: tke Attributes: units: 1 long_name: Turbulent Kinetic Energy Vector Components coverage_content_type: coordinate

tau = <xarray.DataArray 'tau' (tau: 3)> array(['upvp_', 'upwp_', 'vpwp_'], dtype='<U5') Dimensions without coordinates: tau Attributes: units: 1 long_name: Reynolds Stress Vector Components coverage_content_type: coordinate

S = <xarray.DataArray 'S' (S: 3)> array(['Sxx', 'Syy', 'Szz'], dtype='<U3') Dimensions without coordinates: S Attributes: units: 1 long_name: Power Spectral Density Vector Components coverage_content_type: coordinate

C = <xarray.DataArray 'C' (C: 3)> array(['Cxy', 'Cxz', 'Cyz'], dtype='<U3') Dimensions without coordinates: C Attributes: units: 1 long_name: Cross-Spectral Density Vector Components coverage_content_type: coordinate

do_avg(raw_ds, out_ds=None, names=None)[source]

Bin the dataset and calculate the ensemble averages of each variable.

Parameters:

raw_ds (xarray.Dataset) – The raw data structure to be binned
out_ds (xarray.Dataset) – The bin’d (output) data object to which averaged data is added.
names (list of strings) – The names of variables to be averaged. If names is None, all data in raw_ds will be binned.

Returns:

out_ds (xarray.Dataset) – The new (or updated when out_ds is not None) dataset with the averages of all the variables in raw_ds.

:raises AttributeError : when out_ds is supplied as input (not None): :raises and the values in out_ds.attrs are inconsistent with: :raises raw_ds.attrs or the properties of this VelBinner (n_bin,: :raises n_fft, fs, etc.):

Notes

raw_ds.attrs are copied to out_ds.attrs. Inconsistencies between the two (when out_ds is specified as input) raise an AttributeError.

do_var(raw_ds, out_ds=None, names=None, suffix='_var')[source]

Bin the dataset and calculate the ensemble variances of each variable. Complementary to do_avg().

Parameters:

raw_ds (xarray.Dataset) – The raw data structure to be binned.
out_ds (xarray.Dataset) – The binned (output) dataset to which variance data is added, nominally dataset output from do_avg()
names (list of strings) – The names of variables of which to calculate variance. If names is None, all data in raw_ds will be binned.

Returns:

out_ds (xarray.Dataset) – The new (or updated when out_ds is not None) dataset with the variance of all the variables in raw_ds.

:raises AttributeError : when out_ds is supplied as input (not None): :raises and the values in out_ds.attrs are inconsistent with: :raises raw_ds.attrs or the properties of this VelBinner (n_bin,: :raises n_fft, fs, etc.):

Notes

raw_ds.attrs are copied to out_ds.attrs. Inconsistencies between the two (when out_ds is specified as input) raise an AttributeError.

calc_coh(veldat1, veldat2, window='hann', debias=True, noise=(0, 0), n_fft_coh=None, n_bin=None)[source]

Calculate coherence between veldat1 and veldat2.

Parameters:

veldat1 (xarray.DataArray) – The first (the longer, if applicable) raw dataArray of which to calculate coherence
veldat2 (xarray.DataArray) – The second (the shorter, if applicable) raw dataArray of which to calculate coherence
window (str) – String indicating the window function to use (default: ‘hanning’)
noise (float) – The white-noise level of the measurement (in the same units as veldat).
n_fft_coh (int) – n_fft of veldat2, number of elements per bin if ‘None’ is taken from VelBinner
n_bin (int) – n_bin of veldat2, number of elements per bin if ‘None’ is taken from VelBinner

Returns:

da (xarray.DataArray) – The coherence between signal veldat1 and veldat2.

Notes

The two velocity inputs do not have to be perfectly synchronized, but they should have the same start and end timestamps.

calc_phase_angle(veldat1, veldat2, window='hann', n_fft_coh=None, n_bin=None)[source]

Calculate the phase difference between two signals as a function of frequency (complimentary to coherence).

Parameters:

veldat1 (xarray.DataArray) – The first (the longer, if applicable) raw dataArray of which to calculate phase angle
veldat2 (xarray.DataArray) – The second (the shorter, if applicable) raw dataArray of which to calculate phase angle
window (str) – String indicating the window function to use (default: ‘hanning’).
n_fft (int) – Number of elements per bin if ‘None’ is taken from VelBinner
n_bin (int) – Number of elements per bin from veldat2 if ‘None’ is taken from VelBinner

Returns:

da (xarray.DataArray) – The phase difference between signal veldat1 and veldat2.

Notes

The two velocity inputs do not have to be perfectly synchronized, but they should have the same start and end timestamps.

calc_acov(veldat, n_bin=None)[source]

Calculate the auto-covariance of the raw-signal veldat

Parameters:

veldat (xarray.DataArray) – The raw dataArray of which to calculate auto-covariance
n_bin (float) – Number of data elements to use

Returns:

da (xarray.DataArray) – The auto-covariance of veldat

Notes

As opposed to calc_xcov, which returns the full cross-covariance between two arrays, this function only returns a quarter of the full auto-covariance. It computes the auto-covariance over half of the range, then averages the two sides (to return a ‘quartered’ covariance).

This has the advantage that the 0 index is actually zero-lag.

calc_xcov(veldat1, veldat2, npt=1, n_bin=None, normed=False)[source]

Calculate the cross-covariance between arrays veldat1 and veldat2

Parameters:

veldat1 (xarray.DataArray) – The first raw dataArray of which to calculate cross-covariance
veldat2 (xarray.DataArray) – The second raw dataArray of which to calculate cross-covariance
npt (int) – Number of timesteps (lag) to calculate covariance
n_fft (int) – n_fft of veldat2, number of elements per bin if ‘None’ is taken from VelBinner
n_bin (int) – n_bin of veldat2, number of elements per bin if ‘None’ is taken from VelBinner

Returns:

da (xarray.DataArray) – The cross-covariance between signal veldat1 and veldat2.

Notes

The two velocity inputs must be the same length

calc_tke(veldat, noise=None, detrend=True)[source]

Calculate the turbulent kinetic energy (TKE) (variances of u,v,w).

Parameters:

veldat (xarray.DataArray) – Velocity data array from ADV or single beam from ADCP. The last dimension is assumed to be time.
noise (float or array-like) – A vector of the noise levels of the velocity data with the same first dimension as the velocity vector.
detrend (bool (default: False)) – Detrend the velocity data (True), or simply de-mean it (False), prior to computing tke. Note: the psd routines use detrend, so if you want to have the same amount of variance here as there use detrend=True.

Returns:

tke_vec (xarray.DataArray) – dataArray containing u’u’_, v’v’_ and w’w’_

calc_psd(veldat, freq_units='rad/s', fs=None, window='hann', noise=None, n_bin=None, n_fft=None, n_pad=None, step=None)[source]

Calculate the power spectral density of velocity.

Parameters:

veldat (xr.DataArray) – The raw velocity data (of dims ‘dir’ and ‘time’).
freq_units (string) – Frequency units of the returned spectra in either Hz or rad/s (f or \(\omega\))
fs (float (optional)) – The sample rate (default: from the binner).
window (string or array) – Specify the window function. Options: 1, None, ‘hann’, ‘hamm’
noise (float or array-like) – A vector of the noise levels of the velocity data with the same first dimension as the velocity vector.
n_bin (int (optional)) – The bin-size (default: from the binner).
n_fft (int (optional)) – The fft size (default: from the binner).
n_pad (int (optional)) – The number of values to pad with zero (default: 0)
step (int (optional)) – Controls amount of overlap in fft (default: the step size is chosen to maximize data use, minimize nens, and have a minimum of 50% overlap.).

Returns:

psd (xarray.DataArray (3, M, N_FFT)) – The spectra in the ‘u’, ‘v’, and ‘w’ directions.

class dolfyn.adv.turbulence.ADVBinner(n_bin, fs, n_fft=None, n_fft_coh=None, noise=[0, 0, 0])[source]

Bases: VelBinner

A class that builds upon VelBinner for calculating turbulence statistics and velocity spectra from ADV data

Parameters:

n_bin (int) – The length of each bin, in number of points, for this averaging operator.
fs (int) – Instrument sampling frequency in Hz
n_fft (int (optional, default: n_fft = n_bin)) – The length of the FFT for computing spectra (must be <= n_bin)
n_fft_coh (int (optional, default: n_fft_coh`=`n_fft)) – Number of data points to use for coherence and cross-spectra ffts
noise (float, list or numpy.ndarray) – Instrument’s doppler noise in same units as velocity

Initialize an averaging object

Parameters:

n_bin (int) – Number of data points to include in a ‘bin’ (ensemble), not the number of bins
fs (int) – Instrument sampling frequency in Hz
n_fft (int) – Number of data points to use for fft (n_fft`<=`n_bin). Default: n_fft`=`n_bin
n_fft_coh (int) – Number of data points to use for coherence and cross-spectra ffts Default: n_fft_coh`=`n_fft
noise (float, list or numpy.ndarray) – Instrument’s doppler noise in same units as velocity

__call__(ds, freq_units='rad/s', window='hann')[source]: Call self as a function.

calc_stress(veldat, detrend=True)[source]

Calculate the stresses (covariances of u,v,w)

Parameters:

veldat (xr.DataArray) – Velocity data array from ADV data. The last dimension is assumed to be time.
detrend (bool (default: True)) – detrend the velocity data (True), or simply de-mean it (False), prior to computing stress. Note: the psd routines use detrend, so if you want to have the same amount of variance here as there use detrend=True.

Returns:

out (xarray.DataArray)

calc_csd(veldat, freq_units='rad/s', fs=None, window='hann', n_bin=None, n_fft_coh=None)[source]

Calculate the cross-spectral density of velocity components.

Parameters:

veldat (xarray.DataArray) – The raw 3D velocity data.
freq_units (string) – Frequency units of the returned spectra in either Hz or rad/s (f or \(\omega\))
fs (float (optional)) – The sample rate (default: from the binner).
window (string or array) – Specify the window function. Options: 1, None, ‘hann’, ‘hamm’
n_bin (int (optional)) – The bin-size (default: from the binner).
n_fft_coh (int (optional)) – The fft size (default: n_fft_coh from the binner).

Returns:

csd (xarray.DataArray (3, M, N_FFT)) – The first-dimension of the cross-spectrum is the three different cross-spectra: ‘uv’, ‘uw’, ‘vw’.

calc_doppler_noise(psd, pct_fN=0.8)[source]

Calculate bias due to Doppler noise using the noise floor of the velocity spectra.

Parameters:

psd (xarray.DataArray (dir, time, f)) – The ADV velocity spectra
pct_fN (float) – Percent of Nyquist frequency to calculate characeristic frequency

Returns:

doppler_noise (xarray.DataArray) – Doppler noise level in units of m/s

Notes

Approximates bias from

where :math: sigma_{noise} is the bias due to Doppler noise, N is the constant variance or spectral density, and f_{c} is the characteristic frequency.

The characteristic frequency is then found as

where f_{s}/2 is the Nyquist frequency.

Richard, Jean-Baptiste, et al. “Method for identification of Doppler noise levels in turbulent flow measurements dedicated to tidal energy.” International Journal of Marine Energy 3 (2013): 52-64.

Thiébaut, Maxime, et al. “Investigating the flow dynamics and turbulence at a tidal-stream energy site in a highly energetic estuary.” Renewable Energy 195 (2022): 252-262.

check_turbulence_cascade_slope(psd, freq_range=[6.28, 12.57])[source]

This function calculates the slope of the PSD, the power spectra of velocity, within the given frequency range. The purpose of this function is to check that the region of the PSD containing the isotropic turbulence cascade decreases at a rate of \(f^{-5/3}\).

Parameters:

psd (xarray.DataArray ([time,] freq)) – The power spectral density (1D or 2D)
freq_range (iterable(2) (default: [6.28, 12.57])) – The range over which the isotropic turbulence cascade occurs, in units of the psd frequency vector (Hz or rad/s)

Returns:

(m, b) (tuple (slope, y-intercept)) – A tuple containing the coefficients of the log-adjusted linear regression between PSD and frequency

Notes

Calculates slope based on the standard formula for dissipation:

\[S(k) = \alpha \epsilon^{2/3} k^{-5/3} + N\]

The slope of the isotropic turbulence cascade, which should be equal to \(k^{-5/3}\) or \(f^{-5/3}\), where k and f are the wavenumber and frequency vectors, is estimated using linear regression with a log transformation:

\[log10(y) = m*log10(x) + b\]

Which is equivalent to

\[y = 10^{b} x^{m}\]

Where \(y\) is S(k) or S(f), \(x\) is k or f, \(m\) is the slope (ideally -5/3), and \(10^{b}\) is the intercept of y at x^m=1.

calc_epsilon_LT83(psd, U_mag, freq_range=[6.28, 12.57])[source]

Calculate the dissipation rate from the PSD

Parameters:

psd (xarray.DataArray (...,time,freq)) – The power spectral density
U_mag (xarray.DataArray (...,time)) – The bin-averaged horizontal velocity [m s-1] (from dataset shortcut)
freq_range (iterable(2) (default: [6.28, 12.57])) – The range over which to integrate/average the spectrum, in units of the psd frequency vector (Hz or rad/s)

Returns:

epsilon (xarray.DataArray (…,n_time)) – dataArray of the dissipation rate

Notes

This uses the standard formula for dissipation:

\[S(k) = \alpha \epsilon^{2/3} k^{-5/3} + N\]

where \(\alpha = 0.5\) (1.5 for all three velocity components), k is wavenumber, S(k) is the turbulent kinetic energy spectrum, and `N’ is the doppler noise level associated with the TKE spectrum.

With \(k \rightarrow \omega / U\), then – to preserve variance – \(S(k) = U S(\omega)\), and so this becomes:

\[S(\omega) = \alpha \epsilon^{2/3} \omega^{-5/3} U^{2/3} + N\]

With \(k \rightarrow (2\pi f) / U\), then

\[S(\omega) = \alpha \epsilon^{2/3} f^{-5/3} (U/(2\pi))^{2/3} + N\]

LT83 : Lumley and Terray, “Kinematics of turbulence convected by a random wave field”. JPO, 1983, vol13, pp2000-2007.

calc_epsilon_SF(vel_raw, U_mag, fs=None, freq_range=[2.0, 4.0])[source]

Calculate dissipation rate using the “structure function” (SF) method

Parameters:

vel_raw (xarray.DataArray) – The raw velocity data (1D dimension time) upon which to perform the SF technique.
U_mag (xarray.DataArray) – The bin-averaged horizontal velocity (from dataset shortcut)
fs (float) – The sample rate of vel_raw [Hz]
freq_range (iterable(2) (default: [2., 4.])) – The frequency range over which to compute the SF [Hz] (i.e. the frequency range within which the isotropic turbulence cascade falls)

Returns:

epsilon (xarray.DataArray) – dataArray of the dissipation rate

calc_epsilon_TE01(dat_raw, dat_avg, freq_range=[6.28, 12.57])[source]

Calculate the dissipation rate according to TE01.

Parameters:

dat_raw (xarray.Dataset) – The raw (off the instrument) adv dataset
dat_avg (xarray.Dataset) – The bin-averaged adv dataset (calc’d from ‘calc_turbulence’ or ‘do_avg’). The spectra (psd) and basic turbulence statistics (‘tke_vec’ and ‘stress_vec’) must already be computed.
freq_range (iterable(2) (default: [6.28, 12.57])) – The range over which to integrate/average the spectrum, in units of the psd frequency vector (Hz or rad/s)

Notes

TE01 : Trowbridge, J and Elgar, S, “Turbulence measurements in the Surf Zone”. JPO, 2001, vol31, pp2403-2417.

calc_L_int(a_cov, U_mag, fs=None)[source]

Calculate integral length scales.

Parameters:

a_cov (xarray.DataArray) – The auto-covariance array (i.e. computed using calc_acov).
U_mag (xarray.DataArray) – The bin-averaged horizontal velocity (from dataset shortcut)
fs (float) – The raw sample rate

Returns:

L_int (|np.ndarray| (…, n_time)) – The integral length scale (T_int*U_mag).

Notes

The integral time scale (T_int) is the lag-time at which the auto-covariance falls to 1/e.

If T_int is not reached, L_int will default to ‘0’.

dolfyn.adv.turbulence.calc_turbulence(ds_raw, n_bin, fs, n_fft=None, freq_units='rad/s', window='hann')[source]

Functional version of ADVBinner that computes a suite of turbulence statistics for the input dataset, and returns a binned data object.

Parameters:

ds_raw (xarray.Dataset) – The raw adv datset to bin, average and compute turbulence statistics of.
freq_units (string (default: rad/s)) – Frequency units of the returned spectra in either Hz or rad/s (f or \(\omega\))
window (1, None, 'hann', 'hamm') – The window to use for calculating power spectral densities

Returns:

ds (xarray.Dataset) – Returns an ‘binned’ (i.e. ‘averaged’) data object. All fields (variables) of the input data object are averaged in n_bin chunks. This object also computes the following items over those chunks:

tke_vec : The energy in each component, each components is alternatively accessible as: upup_, vpvp_, wpwp_)
stress_vec : The Reynolds stresses, each component is alternatively accessible as: upwp_, vpwp_, upvp_)
U_std : The standard deviation of the horizontal velocity U_mag.
psd : DataArray containing the spectra of the velocity in radial frequency units. The data-array contains: - vel : the velocity spectra array (m^2/s/rad)) - omega : the radial frequncy (rad/s)

class dolfyn.adp.turbulence.ADPBinner(n_bin, fs, n_fft=None, n_fft_coh=None, noise=None, orientation='up', diff_style='centered_extended')[source]

Bases: VelBinner

A class for calculating turbulence statistics from ADCP data

Parameters:

n_bin (int) – Number of data points to include in a ‘bin’ (ensemble), not the number of bins
fs (int) – Instrument sampling frequency in Hz
n_fft (int) – Number of data points to use for fft (n_fft`<=`n_bin). Default: n_fft`=`n_bin
n_fft_coh (int) – Number of data points to use for coherence and cross-spectra ffts Default: n_fft_coh`=`n_fft
noise (float, list or numpy.ndarray) – Instrument’s doppler noise in same units as velocity
orientation (str, default='up') – Instrument’s orientation, either ‘up’ or ‘down’
diff_style (str, default='centered_extended') –
Style of numerical differentiation using Newton’s Method. Either ‘first’ (first difference), ‘centered’ (centered difference), or ‘centered_extended’ (centered difference with first and last points

extended using a first difference).

calc_dudz(vel, orientation=None)[source]

The shear in the first velocity component.

Parameters:

vel (xarray.DataArray) – ADCP raw velocity
orientation (str, default=ADPBinner.orientation) – Direction ADCP is facing (‘up’ or ‘down’)

Notes

The derivative direction is along the profiler’s ‘z’ coordinate (‘dz’ is actually diff(self[‘range’])), not necessarily the ‘true vertical’ direction.

calc_dvdz(vel)[source]

The shear in the second velocity component.

Parameters:: vel (xarray.DataArray) – ADCP raw velocity

Notes

The derivative direction is along the profiler’s ‘z’ coordinate (‘dz’ is actually diff(self[‘range’])), not necessarily the ‘true vertical’ direction.

calc_dwdz(vel)[source]

The shear in the third velocity component.

Parameters:: vel (xarray.DataArray) – ADCP raw velocity

Notes

The derivative direction is along the profiler’s ‘z’ coordinate (‘dz’ is actually diff(self[‘range’])), not necessarily the ‘true vertical’ direction.

calc_shear2(vel)[source]

The horizontal shear squared.

Parameters:: vel (xarray.DataArray) – ADCP raw velocity

Notes

This is actually (dudz)^2 + (dvdz)^2. So, if those variables are not actually vertical derivatives of the horizontal velocity, then this is not the ‘horizontal shear squared’.