Skip to content

Commit

Permalink
Various Updates (#127)
Browse files Browse the repository at this point in the history
* Add function

* Use better data for test

* Ensure that noise subtraction is working properly

* Correct error

* Update test files with noise-subtracted variables

* Should be using along-beam data for adcp turbulence test

* Restore one line

* Fix doppler noise variable coordinate names

* Need to handle 'time' or 'time_b5'

* Update changelog

* Fix window check

* Update documentation

* Workarounds for dual profiling instrument

* dual profile final fixes

* Refactor 'create_dataset' function

* Reinstate some skipped ping logic

* Add ability to read ID 31, clean up altimeter attrs

* Handle variable bottom track beams_cy

* Update test attributes

* Update changelog

* Cleanup

* Update notebook

* Revert test file

* Revert environment change

* Attempt to fix things

* Don't add extra dims

* Add test data

* Small fixes

* Remove unused code

* Update dependency

* More clarity on tke in notebooks

* Git lfs pointer warning

* Fix test data

* Add requirements to conda env
  • Loading branch information
jmcvey3 authored May 8, 2024
1 parent 146049a commit fa25579
Show file tree
Hide file tree
Showing 28 changed files with 4,545 additions and 2,460 deletions.
6 changes: 6 additions & 0 deletions changelog.md
Original file line number Diff line number Diff line change
Expand Up @@ -13,10 +13,16 @@ and this project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.
- Fix netCDF4 compression encoding
- Retain prior netCDF4 variable encoding
- Fix bug in reading raw Nortek Signature altimeter data
- Fix bug where noise input wasn't being subtracted from auto-spectra
- Fix bug that would error out when entering custom FFT window

- API/Useability
- Updates to support python 3.10 and 3.11
- Added ability to read Nortek AWAC waves data
- Added ability to subtract Doppler noise in TKE dissipation rate functions
- Added function to calculate turbulence intensity and remove noise
- Add ability to read Nortek dual profiling instruments
- Add ability to read ID 31 (initial altimeter scan for averaged altimeter measurements)

- Nortek Vectrino (.vno)
- Add support for Nortek Vectrino (.vno) files.
Expand Down
5,549 changes: 3,967 additions & 1,582 deletions docs/ADCP_Example.ipynb

Large diffs are not rendered by default.

611 changes: 41 additions & 570 deletions docs/ADV_Example.ipynb

Large diffs are not rendered by default.

6 changes: 3 additions & 3 deletions docs/apidoc/dolfyn.binners.rst
Original file line number Diff line number Diff line change
Expand Up @@ -21,12 +21,12 @@ Below is a list of functions that can be called from `VelBinner`.
~dolfyn.binned.TimeBinner.std
~dolfyn.velocity.VelBinner.do_avg
~dolfyn.velocity.VelBinner.do_var
~dolfyn.velocity.VelBinner.calc_ti
~dolfyn.velocity.VelBinner.calc_psd
~dolfyn.velocity.VelBinner.calc_coh
~dolfyn.velocity.VelBinner.calc_phase_angle
~dolfyn.velocity.VelBinner.calc_acov
~dolfyn.velocity.VelBinner.calc_xcov
~dolfyn.velocity.VelBinner.calc_tke
~dolfyn.velocity.VelBinner.calc_psd
~dolfyn.binned.TimeBinner.calc_freq
~dolfyn.binned.TimeBinner.calc_psd_base
~dolfyn.binned.TimeBinner.calc_csd_base
Expand All @@ -43,6 +43,7 @@ Functions for analyzing ADV data via the `ADVBinner` class, beyond those describ
~dolfyn.adv.turbulence.ADVBinner
~dolfyn.adv.turbulence.calc_turbulence
~dolfyn.adv.turbulence.ADVBinner.calc_csd
~dolfyn.velocity.VelBinner.calc_tke
~dolfyn.adv.turbulence.ADVBinner.calc_stress
~dolfyn.adv.turbulence.ADVBinner.calc_doppler_noise
~dolfyn.adv.turbulence.ADVBinner.check_turbulence_cascade_slope
Expand All @@ -64,7 +65,6 @@ Functions for analyzing ADCP data via the `ADPBinner` class, beyond those descri
~dolfyn.adp.turbulence.ADPBinner.calc_doppler_noise
~dolfyn.adp.turbulence.ADPBinner.calc_stress_4beam
~dolfyn.adp.turbulence.ADPBinner.calc_stress_5beam
~dolfyn.adp.turbulence.ADPBinner.calc_total_tke
~dolfyn.adp.turbulence.ADPBinner.check_turbulence_cascade_slope
~dolfyn.adp.turbulence.ADPBinner.calc_dissipation_LT83
~dolfyn.adp.turbulence.ADPBinner.calc_dissipation_SF
Expand Down
74 changes: 26 additions & 48 deletions dolfyn/adp/turbulence.py
Original file line number Diff line number Diff line change
Expand Up @@ -238,9 +238,10 @@ def calc_doppler_noise(self, psd, pct_fN=0.8):
N2 = psd.sel(freq=f_range) * psd.freq.sel(freq=f_range)
noise_level = np.sqrt(N2.mean(dim='freq'))

time_coord = psd.dims[0] # no reason this shouldn't be time or time_b5
return xr.DataArray(
noise_level.values.astype('float32'),
dims=['time'],
coords={time_coord: psd.coords[time_coord]},
attrs={'units': 'm s-1',
'long_name': 'Doppler Noise Level',
'description': 'Doppler noise level calculated '
Expand Down Expand Up @@ -439,7 +440,7 @@ def calc_stress_5beam(self, ds, noise=None, orientation=None, beam_angle=None, t
in pitch and roll. u'v'_ cannot be directly calculated by a 5-beam ADCP,
so it is approximated by the covariance of `u` and `v`. The uncertainty
introduced by using this approximation is small if deviations from pitch
and roll are small (< 10 degrees).
and roll are small (<= 5 degrees).
Dewey, R., and S. Stringer. "Reynolds stresses and turbulent kinetic
energy estimates from various ADCP beam configurations: Theory." J. of
Expand All @@ -455,7 +456,7 @@ def calc_stress_5beam(self, ds, noise=None, orientation=None, beam_angle=None, t

# Run through warnings
b_angle, noise = self._stress_func_warnings(
ds, beam_angle, noise, tilt_thresh=10)
ds, beam_angle, noise, tilt_thresh=5)

# Fetch beam order
beam_order, phi2, phi3 = self._check_orientation(
Expand Down Expand Up @@ -522,47 +523,6 @@ def calc_stress_5beam(self, ds, noise=None, orientation=None, beam_angle=None, t

return tke_vec, stress_vec

def calc_total_tke(self, ds, noise=None, orientation=None, beam_angle=None):
"""Calculate magnitude of turbulent kinetic energy from 5-beam ADCP.
Parameters
----------
ds : xarray.Dataset
Raw dataset in beam coordinates
ds_avg : xarray.Dataset
Binned dataset in final coordinate reference frame
noise : int or xarray.DataArray, default=0 (time)
Doppler noise level in units of m/s
orientation : str, default=ds.attrs['orientation']
Direction ADCP is facing ('up' or 'down')
beam_angle : int, default=ds.attrs['beam_angle']
ADCP beam angle in units of degrees
Returns
-------
tke : xarray.DataArray
Turbulent kinetic energy magnitude
Notes
-----
This function is a wrapper around 'calc_stress_5beam' that then
combines the TKE components.
Warning: the integral length scale of turbulence captured by the
ADCP measurements (i.e. the size of turbulent structures) increases
with increasing range from the instrument.
"""

tke_vec = self.calc_stress_5beam(
ds, noise, orientation, beam_angle, tke_only=True)

tke = tke_vec.sum('tke') / 2
tke.attrs['units'] = 'm2 s-2'
tke.attrs['long_name'] = 'TKE Magnitude',
tke.attrs['standard_name'] = 'specific_turbulent_kinetic_energy_of_sea_water'

return tke.astype('float32')

def check_turbulence_cascade_slope(self, psd, freq_range=[0.2, 0.4]):
"""This function calculates the slope of the PSD, the power spectra
of velocity, within the given frequency range. The purpose of this
Expand Down Expand Up @@ -623,7 +583,7 @@ def check_turbulence_cascade_slope(self, psd, freq_range=[0.2, 0.4]):

return m, b

def calc_dissipation_LT83(self, psd, U_mag, freq_range=[0.2, 0.4]):
def calc_dissipation_LT83(self, psd, U_mag, freq_range=[0.2, 0.4], noise=None):
"""Calculate the TKE dissipation rate from the velocity spectra.
Parameters
Expand All @@ -633,8 +593,12 @@ def calc_dissipation_LT83(self, psd, U_mag, freq_range=[0.2, 0.4]):
U_mag : xarray.DataArray (time)
The bin-averaged horizontal velocity (a.k.a. speed) from a single depth bin (range)
f_range : iterable(2)
The range over which to integrate/average the spectrum, in units
The range over which to integrate/average the spectrum, in units
of the psd frequency vector (Hz or rad/s)
noise : float or array-like
Instrument noise level in same units as velocity. Typically
found from `adp.turbulence.calc_doppler_noise`.
Default: None.
Returns
-------
Expand Down Expand Up @@ -669,6 +633,20 @@ def calc_dissipation_LT83(self, psd, U_mag, freq_range=[0.2, 0.4]):
raise Exception('PSD should be 2-dimensional (time, frequency)')
if len(U_mag.shape) != 1:
raise Exception('U_mag should be 1-dimensional (time)')
if not hasattr(freq_range, "__iter__") or len(freq_range) != 2:
raise ValueError("`freq_range` must be an iterable of length 2.")
if noise is not None:
if np.shape(noise)[0] != np.shape(psd)[0]:
raise Exception(
'Noise should have same first dimension as PSD')
else:
noise = np.array(0)

# Noise subtraction from binner.TimeBinner.calc_psd_base
psd = psd.copy()
if noise is not None:
psd -= noise**2 / (self.fs / 2)
psd = psd.where(psd > 0, np.min(np.abs(psd)) / 100)

freq = psd.freq
idx = np.where((freq_range[0] < freq) & (freq < freq_range[1]))
Expand Down Expand Up @@ -853,7 +831,7 @@ def calc_ustar_fit(self, ds_avg, upwp_, z_inds=slice(1, 5), H=None):
"""

if not H:
H = ds_avg.depth.values
H = ds_avg["depth"].values
z = ds_avg['range'].values
upwp_ = upwp_.values

Expand All @@ -863,6 +841,6 @@ def calc_ustar_fit(self, ds_avg, upwp_, z_inds=slice(1, 5), H=None):

return xr.DataArray(
u_star.astype('float32'),
coords={'time': ds_avg.time},
coords={'time': ds_avg["time"]},
attrs={'units': 'm s-1',
'long_name': 'Friction Velocity'})
35 changes: 28 additions & 7 deletions dolfyn/adv/turbulence.py
Original file line number Diff line number Diff line change
Expand Up @@ -21,8 +21,10 @@ class ADVBinner(VelBinner):
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
noise : float or array-like
Instrument noise level in same units as velocity. Typically
found from `adv.turbulence.calc_doppler_noise`.
Default: None.
"""

def __call__(self, ds, freq_units='rad/s', window='hann'):
Expand All @@ -49,7 +51,7 @@ def __call__(self, ds, freq_units='rad/s', window='hann'):

def calc_stress(self, veldat, detrend=True):
"""
Calculate the stresses (covariances of u,v,w)
Calculate the stresses (covariances of u,v,w in m^2/s^2)
Parameters
----------
Expand Down Expand Up @@ -165,7 +167,7 @@ def calc_csd(self, veldat,
'time': time,
'coh_freq': coh_freq},
dims=['C', 'time', 'coh_freq'],
attrs={'units': units,
attrs={'units': units,
'n_fft_coh': n_fft,
'long_name': 'Cross Spectral Density'})
csd['coh_freq'].attrs['units'] = freq_units
Expand Down Expand Up @@ -230,7 +232,7 @@ def calc_doppler_noise(self, psd, pct_fN=0.8):

return xr.DataArray(
noise_level.values.astype('float32'),
dims=['dir', 'time'],
coords={'S': psd['S'], 'time': psd['time']},
attrs={'units': 'm/s',
'long_name': 'Doppler Noise Level',
'description': 'Doppler noise level calculated '
Expand Down Expand Up @@ -296,7 +298,7 @@ def check_turbulence_cascade_slope(self, psd, freq_range=[6.28, 12.57]):

return m, b

def calc_epsilon_LT83(self, psd, U_mag, freq_range=[6.28, 12.57]):
def calc_epsilon_LT83(self, psd, U_mag, freq_range=[6.28, 12.57], noise=None):
"""Calculate the dissipation rate from the PSD
Parameters
Expand All @@ -308,6 +310,10 @@ def calc_epsilon_LT83(self, psd, U_mag, freq_range=[6.28, 12.57]):
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)
noise : float or array-like
Instrument noise level in same units as velocity. Typically
found from `adv.turbulence.calc_doppler_noise`.
Default: None.
Returns
-------
Expand Down Expand Up @@ -339,8 +345,23 @@ def calc_epsilon_LT83(self, psd, U_mag, freq_range=[6.28, 12.57]):
"""

# Ensure time has been averaged
if len(psd.time)!=len(U_mag.time):
if len(psd.time) != len(U_mag.time):
raise Exception("`U_mag` should be from ensembled-averaged dataset")
if not hasattr(freq_range, "__iter__") or len(freq_range) != 2:
raise ValueError("`freq_range` must be an iterable of length 2.")

if noise is not None:
if np.shape(noise)[0] != 3:
raise Exception(
'Noise should have same first dimension as velocity')
else:
noise = np.array([0, 0, 0])[:, None, None]

# Noise subtraction from binner.TimeBinner.calc_psd_base
psd = psd.copy()
if noise is not None:
psd -= (noise**2 / (self.fs / 2))
psd = psd.where(psd > 0, np.min(np.abs(psd)) / 100)

freq = psd.freq
idx = np.where((freq_range[0] < freq) & (freq < freq_range[1]))
Expand Down
2 changes: 1 addition & 1 deletion dolfyn/binned.py
Original file line number Diff line number Diff line change
Expand Up @@ -401,7 +401,7 @@ def calc_psd_base(self, dat, fs=None, window='hann', noise=0,
for slc in slice1d_along_axis(dat.shape, -1):
out[slc] = psd(dat[slc], n_fft, fs,
window=window, step=step)
if noise != 0:
if np.any(noise):
out -= noise**2 / (fs/2)
# Make sure all values of the PSD are >0 (but still small):
out[out < 0] = np.min(np.abs(out)) / 100
Expand Down
3 changes: 3 additions & 0 deletions dolfyn/example_data/dual_profile.ad2cp
Git LFS file not shown
8 changes: 7 additions & 1 deletion dolfyn/io/api.py
Original file line number Diff line number Diff line change
Expand Up @@ -192,7 +192,7 @@ def save(ds, filename,

if compression:
# New netcdf4-c cannot compress variable length strings
if isinstance(ds[ky].data[0], str):
if ds[ky].size <= 1 or isinstance(ds[ky].data[0], str):
continue
enc[ky].update(dict(zlib=True, complevel=1))

Expand All @@ -219,6 +219,12 @@ def load(filename):
"""

filename = _check_file_ext(filename, 'nc')

file_type = _get_filetype(filename)
if file_type == '<GIT-LFS pointer>':
raise IOError("File '{}' looks like a git-lfs pointer. You may need to "
"install and initialize git-lfs. See https://git-lfs.github.com"
" for details.".format(filename))

ds = xr.load_dataset(filename, engine='netcdf4')

Expand Down
Loading

0 comments on commit fa25579

Please sign in to comment.