Reading ADV Data with MHKiT

This example presents a simplified workflow for analyzing Acoustic Doppler Velocimetry (ADV) data using MHKiT. MHKiT incorporates the DOLfYN codebase as a module to handle ADV and Acoustic Doppler Current Profiler (ADCP) data.

A standard ADV data analysis workflow can be segmented into the following steps:

1. Raw Data Review: Evaluate the original data by verifying timestamps and assessing the quality of velocity data, specifically looking for any data spikes.

2. Data Cleaning: Identify and eliminate any spurious data points. If needed, bad data points can be replaced with interpolated values.

3. Data Rotation: Transform the data into the principal flow coordinates, which are the streamwise, cross-stream, and vertical directions.

4. Data Averaging: Aggregate the data into bins or ensembles, each of which spans a predefined time length, typically between 5 and 10 minutes.

5. Statistical Analysis: Compute turbulence statistics such as turbulence intensity, Turbulent Kinetic Energy (TKE), and Reynolds stresses for the observed flow field.

Start your analysis by importing the necessary tools:

from mhkit import dolfyn
from mhkit.dolfyn.adv import api

c:\Users\mcve343\Anaconda3\lib\site-packages\xarray\backends\cfgrib_.py:29: UserWarning: Failed to load cfgrib - most likely there is a problem accessing the ecCodes library. Try `import cfgrib` to get the full error message
  warnings.warn(

Read Raw Instrument Data

DOLfYN currently only carries support for the Nortek Vector ADV. The example loaded here is a short clip of data from a test deployment to show DOLfYN’s capabilities.

Start by reading in the raw datafile downloaded from the instrument. The dolfyn.read function reads the raw file and dumps the information into an xarray Dataset, which contains three groups of variables:

1. Velocity, amplitude, and correlation of the Doppler velocimetry

2. Measurements of the instrument’s bearing and environment

3. Orientation matrices DOLfYN uses for rotating through coordinate frames.

ds = dolfyn.read("data/dolfyn/vector_data01.VEC")

Reading file data/dolfyn/vector_data01.VEC ...

There are two ways to see what’s in a Dataset. The first is to simply type the dataset’s name to see the standard xarray output. To access a particular variable in a dataset, use dict-style (ds['vel']) or attribute-style syntax (ds.vel). See the xarray docs for more details on how to use the xarray format.

# print the dataset
ds

<xarray.Dataset>
Dimensions:              (x1: 3, x2: 3, time: 122912, dir: 3, beam: 3,
                          earth: 3, inst: 3)
Coordinates:
  * x1                   (x1) int32 1 2 3
  * x2                   (x2) int32 1 2 3
  * time                 (time) datetime64[ns] 2012-06-12T12:00:02.968749046 ...
  * dir                  (dir) <U1 'X' 'Y' 'Z'
  * beam                 (beam) int32 1 2 3
  * earth                (earth) <U1 'E' 'N' 'U'
  * inst                 (inst) <U1 'X' 'Y' 'Z'
Data variables: (12/15)
    beam2inst_orientmat  (x1, x2) float64 2.709 -1.34 -1.364 ... -0.3438 -0.3499
    batt                 (time) float32 13.2 13.2 13.2 13.2 ... nan nan nan nan
    c_sound              (time) float32 1.493e+03 1.493e+03 ... nan nan
    heading              (time) float32 5.6 10.5 10.51 10.52 ... nan nan nan nan
    pitch                (time) float32 -31.5 -31.7 -31.69 ... nan nan nan
    roll                 (time) float32 0.4 4.2 4.253 4.306 ... nan nan nan nan
    ...                   ...
    orientation_down     (time) bool True True True True ... True True True True
    vel                  (dir, time) float32 -1.002 -1.008 -0.944 ... nan nan
    amp                  (beam, time) uint8 104 110 111 113 108 ... 0 0 0 0 0
    corr                 (beam, time) uint8 97 91 97 98 90 95 95 ... 0 0 0 0 0 0
    pressure             (time) float64 5.448 5.436 5.484 5.448 ... 0.0 0.0 0.0
    orientmat            (earth, inst, time) float32 0.0832 0.155 ... -0.7065
Attributes: (12/39)
    inst_make:                   Nortek
    inst_model:                  Vector
    inst_type:                   ADV
    rotate_vars:                 ['vel']
    n_beams:                     3
    profile_mode:                continuous
    ...                          ...
    recorder_size_bytes:         4074766336
    vel_range:                   normal
    firmware_version:            3.34
    fs:                          32.0
    coord_sys:                   inst
    has_imu:                     0

xarray.Dataset

• Dimensions:
  • x1: 3
  • x2: 3
  • time: 122912
  • dir: 3
  • beam: 3
  • earth: 3
  • inst: 3
• Coordinates: (7)
  • x1
    (x1)
    int32
    1 2 3

    coverage_content_type :

    coordinate

    array([1, 2, 3])

  • x2
    (x2)
    int32
    1 2 3

    coverage_content_type :

    coordinate

    array([1, 2, 3])

  • time
    (time)
    datetime64[ns]
    2012-06-12T12:00:02.968749046 .....

    coverage_content_type :

    coordinate

    units :

    seconds since 1970-01-01 00:00:00

    long_name :

    Time

    standard_name :

    time

    array(['2012-06-12T12:00:02.968749046', '2012-06-12T12:00:02.999999046',
           '2012-06-12T12:00:03.031249046', ..., '2012-06-12T13:04:03.875000715',
           '2012-06-12T13:04:03.906250715', '2012-06-12T13:04:03.937500715'],
          dtype='datetime64[ns]')

  • dir
    (dir)
    <U1
    'X' 'Y' 'Z'

    units :

    1

    long_name :

    Reference Frame

    coverage_content_type :

    coordinate

    ref_frame :

    inst

    array(['X', 'Y', 'Z'], dtype='<U1')

  • beam
    (beam)
    int32
    1 2 3

    units :

    1

    long_name :

    Beam Reference Frame

    coverage_content_type :

    coordinate

    array([1, 2, 3])

  • earth
    (earth)
    <U1
    'E' 'N' 'U'

    units :

    1

    long_name :

    Earth Reference Frame

    coverage_content_type :

    coordinate

    array(['E', 'N', 'U'], dtype='<U1')

  • inst
    (inst)
    <U1
    'X' 'Y' 'Z'

    units :

    1

    long_name :

    Instrument Reference Frame

    coverage_content_type :

    coordinate

    array(['X', 'Y', 'Z'], dtype='<U1')

• Data variables: (15)
  • beam2inst_orientmat
    (x1, x2)
    float64
    2.709 -1.34 ... -0.3438 -0.3499

    units :

    1

    long_name :

    Rotation Matrix

    array([[ 2.70874023, -1.34008789, -1.36376953],
           [ 0.0144043 , -2.38354492,  2.39550781],
           [-0.3371582 , -0.34375   , -0.34985352]])

  • batt
    (time)
    float32
    13.2 13.2 13.2 13.2 ... nan nan nan

    units :

    V

    long_name :

    Battery Voltage

    array([13.2, 13.2, 13.2, ...,  nan,  nan,  nan], dtype=float32)

  • c_sound
    (time)
    float32
    1.493e+03 1.493e+03 ... nan nan

    units :

    m s-1

    long_name :

    Speed of Sound

    standard_name :

    speed_of_sound_in_sea_water

    array([1492.6, 1492.6, 1492.6, ...,    nan,    nan,    nan], dtype=float32)

  • heading
    (time)
    float32
    5.6 10.5 10.51 ... nan nan nan

    units :

    degree

    long_name :

    Heading

    standard_name :

    platform_orientation

    array([ 5.6   , 10.5   , 10.5125, ...,     nan,     nan,     nan],
          dtype=float32)

  • pitch
    (time)
    float32
    -31.5 -31.7 -31.69 ... nan nan nan

    units :

    degree

    long_name :

    Pitch

    standard_name :

    platform_pitch

    array([-31.5     , -31.7     , -31.690626, ...,        nan,        nan,
                  nan], dtype=float32)

  • roll
    (time)
    float32
    0.4 4.2 4.253 4.306 ... nan nan nan

    units :

    degree

    long_name :

    Roll

    standard_name :

    platform_roll

    array([0.4      , 4.2000003, 4.253125 , ...,       nan,       nan,
                 nan], dtype=float32)

  • temp
    (time)
    float32
    12.67 12.68 12.68 ... nan nan nan

    units :

    degree_C

    long_name :

    Temperature

    standard_name :

    sea_water_temperature

    array([12.67    , 12.679999, 12.679999, ...,       nan,       nan,
                 nan], dtype=float32)

  • error
    (time)
    uint8
    0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0

    units :

    1

    long_name :

    Error Code

    array([0, 0, 0, ..., 0, 0, 0], dtype=uint8)

  • status
    (time)
    uint8
    117 117 0 0 0 0 0 ... 0 0 0 0 0 0 0

    units :

    1

    long_name :

    Status Code

    array([117, 117,   0, ...,   0,   0,   0], dtype=uint8)

  • orientation_down
    (time)
    bool
    True True True ... True True True

    units :

    1

    long_name :

    Orientation of ADV Communication Cable

    array([ True,  True,  True, ...,  True,  True,  True])

  • vel
    (dir, time)
    float32
    -1.002 -1.008 -0.944 ... nan nan

    units :

    m s-1

    long_name :

    Water Velocity

    coverage_content_type :

    physicalMeasurement

    array([[-1.0020001 , -1.008     , -0.94400007, ...,         nan,
                    nan,         nan],
           [ 0.097     ,  0.068     ,  0.066     , ...,         nan,
                    nan,         nan],
           [ 0.115     ,  0.12400001,  0.12200001, ...,         nan,
                    nan,         nan]], dtype=float32)

  • amp
    (beam, time)
    uint8
    104 110 111 113 108 ... 0 0 0 0 0

    units :

    1

    long_name :

    Acoustic Signal Amplitude

    standard_name :

    signal_intensity_from_multibeam_acoustic_doppler_velocity_sensor_in_sea_water

    coverage_content_type :

    physicalMeasurement

    array([[104, 110, 111, ...,   0,   0,   0],
           [109, 111, 118, ...,   0,   0,   0],
           [111, 107, 103, ...,   0,   0,   0]], dtype=uint8)

  • corr
    (beam, time)
    uint8
    97 91 97 98 90 95 ... 0 0 0 0 0 0

    units :

    %

    long_name :

    Acoustic Signal Correlation

    coverage_content_type :

    physicalMeasurement

    array([[97, 91, 97, ...,  0,  0,  0],
           [97, 97, 98, ...,  0,  0,  0],
           [96, 97, 97, ...,  0,  0,  0]], dtype=uint8)

  • pressure
    (time)
    float64
    5.448 5.436 5.484 ... 0.0 0.0 0.0

    units :

    dbar

    long_name :

    Pressure

    standard_name :

    sea_water_pressure

    array([5.448, 5.436, 5.484, ..., 0.   , 0.   , 0.   ])

  • orientmat
    (earth, inst, time)
    float32
    0.0832 0.155 ... -0.7065 -0.7065

    units :

    1

    long_name :

    Orientation Matrix

    array([[[ 0.08320302,  0.155048  ,  0.15524618, ..., -0.8397447 ,
             -0.8397447 , -0.8397447 ],
            [ 0.8485708 ,  0.83656424,  0.8366149 , ...,  0.21873687,
              0.21873687,  0.21873687],
            [-0.52249855, -0.5254717 , -0.52533245, ..., -0.49697396,
             -0.49697396, -0.49697396]],
    
           [[ 0.9948472 ,  0.97360104,  0.9733992 , ...,  0.48448867,
              0.48448867,  0.48448867],
            [-0.10121067, -0.21958625, -0.22025378, ...,  0.7150815 ,
              0.7150815 ,  0.7150815 ],
            [-0.00595242, -0.06231184, -0.06310495, ..., -0.5039139 ,
             -0.5039139 , -0.5039139 ]],
    
           [[-0.05793348, -0.1675142 , -0.16850099, ...,  0.24515234,
              0.24515234,  0.24515234],
            [-0.51931095, -0.50193846, -0.5015614 , ..., -0.6639372 ,
             -0.6639372 , -0.6639372 ],
            [-0.85261935, -0.84852624, -0.84855384, ..., -0.7064614 ,
             -0.7064614 , -0.7064614 ]]], dtype=float32)

• Attributes: (39)

  inst_make :

  Nortek

  inst_model :

  Vector

  inst_type :

  ADV

  rotate_vars :

  ['vel']

  n_beams :

  3

  profile_mode :

  continuous

  burst_mode :

  True

  power_level :

  1

  sync_out_pos :

  end

  sample_on_sync :

  False

  start_on_sync :

  False

  compass_update_rate :

  1

  n_bins :

  1

  deployment_name :

  APLUW_

  wrap_mode :

  False

  analog_in :

  0

  software_version :

  1.32.00

  salinity :

  30.0

  comments :

  APL-UW vector on Tidal Turbulence Mooring in Admiralty, times PDT

  user_specified_sound_speed :

  False

  analog_output :

  False

  output_format :

  Vector

  serial_output :

  False

  power_output_analog :

  False

  n_pings_per_burst :

  1

  pressure_sensor :

  yes

  compass :

  yes

  tilt_sensor :

  yes

  carrier_freq_kHz :

  6000

  serial_number :

  VEC 9062

  ProLogFWver :

  4.08

  PIC_version :

  0

  hardware_rev :

  4

  recorder_size_bytes :

  4074766336

  vel_range :

  normal

  firmware_version :

  3.34

  fs :

  32.0

  coord_sys :

  inst

  has_imu :

  0

A second way provided to look at data is through the DOLfYN view. This view has several convenience methods, shortcuts, and functions built-in. It includes an alternate – and somewhat more informative/compact – description of the data object when in interactive mode. This can be accessed using

ds_dolfyn = ds.velds
ds_dolfyn

<ADV data object>:  Nortek Vector
  . 1.07 hours (started: Jun 12, 2012 12:00)
  . inst-frame
  . (122912 pings @ 32.0Hz)
  Variables:
  - time ('time',)
  - vel ('dir', 'time')
  - orientmat ('earth', 'inst', 'time')
  - heading ('time',)
  - pitch ('time',)
  - roll ('time',)
  - temp ('time',)
  - pressure ('time',)
  - amp ('beam', 'time')
  - corr ('beam', 'time')
  ... and others (see `<obj>.variables`)

Quality Control

ADV velocity data tends to have spikes due to Doppler noise, and the common way to “despike” the data is by using the phase-space algorithm by Goring and Nikora (2002). DOLfYN integrates this function using a 2-step approach: create a logical mask where True corresponds to a spike detection, and then utilize an interpolation function to replace the spikes.

# Clean the file using the Goring+Nikora method:
mask = api.clean.GN2002(ds.vel, npt=5000)
# Replace bad datapoints via cubic spline interpolation
ds["vel"] = api.clean.clean_fill(ds["vel"], mask, npt=12, method="cubic", maxgap=None)

print("Percent of data containing spikes: {0:.2f}%".format(100 * mask.mean()))

# If interpolation isn't desired:
ds_nan = ds.copy(deep=True)
ds_nan.coords["mask"] = (("dir", "time"), ~mask)
ds_nan["vel"] = ds_nan["vel"].where(ds_nan["mask"])

Percent of data containing spikes: 0.73%

Coordinate Rotations

Now that the data has been cleaned, the next step is to rotate the velocity data into true East, North, Up (ENU) coordinates.

ADVs use an internal compass or magnetometer to determine magnetic ENU directions. The set_declination function takes the user supplied magnetic declination (which can be looked up online for specific coordinates) and adjusts the orientation matrix saved within the dataset.

Instruments save vector data in the coordinate system specified in the deployment configuration file. To make the data useful, it must be rotated through coordinate systems (“beam”<->”inst”<->”earth”<->”principal”), done through the rotate2 function. If the “earth” (ENU) coordinate system is specified, DOLfYN will automatically rotate the dataset through the necessary coordinate systems to get there. The inplace set as true will alter the input dataset “in place”, a.k.a. it not create a new dataset.

# First set the magnetic declination
dolfyn.set_declination(
    ds, declin=10, inplace=True
)  # declination points 10 degrees East

# Rotate that data from the instrument to earth frame (ENU):
dolfyn.rotate2(ds, "earth", inplace=True)

Once in the true ENU frame of reference, we can calculate the principal flow direction for the velocity data and rotate it into the principal frame of reference (streamwise, cross-stream, vertical). Principal flow directions are aligned with and orthogonal to the flow streamlines at the measurement location.

First, the principal flow direction must be calculated through calc_principal_heading. As a standard for DOLfYN functions, those that begin with “calc_*” require the velocity data for input. This function is different from others in DOLfYN in that it requires place the output in an attribute called “principal_heading”, as shown below.

Again we use rotate2 to change coordinate systems.

ds.attrs["principal_heading"] = dolfyn.calc_principal_heading(ds["vel"])
dolfyn.rotate2(ds, "principal", inplace=True)

Averaging Data

The next step in ADV analysis is to average the velocity data into time bins (ensembles) and calculate turbulence statistics. These averaged values are then used to calculate turbulence statistics. There are two distinct methods for performing this operation, both of which utilize the same variable inputs and produce identical datasets.

1. Object-Oriented Approach (standard): Define an ‘averaging object’, create a dataset binned in time, and calculate basic turbulence statistics. This is accomplished by initiating an object from the ADVBinner class and then feeding that object with our dataset.

2. Functional Approach (simple): The same operations can be performed using the functional counterpart of ADVBinner, turbulence_statistics.

Function inputs shown here are the dataset itself: - n_bin: the number of elements in each bin; - fs: the ADV’s sampling frequency in Hz; - n_fft: optional, the number of elements per FFT for spectral analysis; - freq_units: optional, either in Hz or rad/s, of the calculated spectral frequency vector.

All of the variables in the returned dataset have been bin-averaged, where each average is computed using the number of elements specified in n_bins. Additional variables in this dataset include the turbulent kinetic energy (TKE) vector (“ds_binned.tke_vec”), the Reynold’s stresses (“ds_binned.stress”), and the power spectral densities (“ds_binned.psd”), calculated for each bin.

# Option 1 (standard)
binner = api.ADVBinner(n_bin=ds.fs * 600, fs=ds.fs, n_fft=1024)
ds_binned = binner.bin_average(ds)

# Option 2 (simple)
# ds_binned = api.calc_turbulence(ds, n_bin=ds.fs*600, fs=ds.fs, n_fft=1024, freq_units="Hz")

The benefit to using ADVBinner is that one has access to all of the velocity and turbulence analysis functions that DOLfYN contains. If basic analysis will suffice, the turbulence_statistics function is the most convienent. Either option can still utilize DOLfYN’s shortcuts.

See the DOLfYN API for the full list of functions and shortcuts. A few examples are shown below.

Some things to know: - All functions operate bin-by-bin. - Some functions will fail if there are NaN’s in the data stream (Notably the PSD functions) - “Shorcuts”, as referred to in DOLfYN, are functions accessible by the xarray accessor velds, as shown below. The list of “shorcuts” available through velds are listed here. Some shorcut variables require the raw dataset, some an averaged dataset.

For instance, - bin_variance calculates the binned-variance of each variable in the raw dataset, the complementary to bin_average. Variables returned by this function contain a “_var” suffix to their name. - cross_spectral_density calculates the cross spectral power density between each direction of the supplied DataArray. Note that inputs specified in creating the ADVBinner object can be overridden or additionally specified for a particular function call. - velds.I is the shortcut for turbulence intensity. This particular shortcut requires a dataset created by bin_average, because it requires bin-averaged data to calculate.

# Calculate the variance of each variable in the dataset and add to the averaged dataset
ds_binned = binner.bin_variance(ds, out_ds=ds_binned)

# Calculate the power spectral density
ds_binned["auto_spectra"] = binner.power_spectral_density(ds["vel"], freq_units="Hz")
# Calculate dissipation rate from isotropic turbulence cascade
ds_binned["dissipation"] = binner.dissipation_rate_LT83(
    ds_binned["auto_spectra"], ds_binned.velds.U_mag, freq_range=[0.5, 1]
)

# Calculate the cross power spectral density
ds_binned["cross_spectra"] = binner.cross_spectral_density(
    ds["vel"], freq_units="Hz", n_fft_coh=512
)

# Calculated the turbulence intensity (requires a binned dataset)
ds_binned["TI"] = ds_binned.velds.I

Plotting can be performed using matplotlib. As an example, the mean spectrum in the streamwise direction is plotted here. This spectrum shows the mean energy density in the flow at a particular flow frequency.

import matplotlib.pyplot as plt

%matplotlib inline

plt.figure()
plt.loglog(ds_binned["freq"], ds_binned["auto_spectra"].sel(S="Sxx").mean(dim="time"))
plt.xlabel("Frequency [Hz]")
plt.ylabel("Energy Density $\mathrm{[m^2/s^s/Hz]}$")
plt.title("Streamwise Direction")

Text(0.5, 1.0, 'Streamwise Direction')

Saving and Loading DOLfYN datasets

Datasets can be saved and reloaded using the save and load functions. Xarray is saved natively in netCDF format, hence the “.nc” extension.

Note: DOLfYN datasets cannot be saved using xarray’s native ds.to_netcdf; however, DOLfYN datasets can be opened using xarray.open_dataset.

# Uncomment these lines to save and load to your current working directory
# dolfyn.save(ds, 'your_data.nc')
# ds_saved = dolfyn.load('your_data.nc')