ADCP Example

The following example shows a typical workflow for analyzing ADCP data using DOLfYN’s tools.

A typical ADCP data workflow is broken down into 1. Review the raw data - Check timestamps - Calculate/check that the depth bin locations are correct - Look at velocity, beam amplitude and/or beam correlation data quality 2. Remove data located above the water surface or below the seafloor 3. Check for spurious datapoints and remove if necessary 4. If not already done within the instrument, average the data into bins of a set time length (normally 5 to 10 min) 5. Conduct further analysis as required

Start by importing the necessary DOLfYN tools through MHKiT:

[1]:

# Import core DOLfYN functions
import dolfyn as dlfn
# Import ADCP-specific API tools
from dolfyn.adp import api

First Steps and QC’ing Data

1.) Set deployment height

Because this is a Nortek instrument, the deployment software doesn’t take into account the deployment height, aka where in the water column the ADCP is. The center of the first depth bin is located at a distance = deployment height + blanking distance + cell size, so the range coordinate needs to be corrected so that ‘0’ corresponds to the seafloor. This can be done in DOLfYN using the set_range_offset function. This same function can be used to account for the depth of down-facing instrument below the water surface.

Note, TRDI’s deployment software asks the user to enter the deployment height/depth, which, if needed, can be adjusted after-the-fact using set_range_offset as well.

[4]:

# The ADCP transducers were measured to be 0.6 m from the feet of the lander
dat = api.clean.set_range_offset(dat, 0.6)

So, the center of bin 1 is located at 1.2 m:

[5]:

dat.range

[5]:

<xarray.DataArray 'range' (range: 28)>
array([ 1.2,  1.7,  2.2,  2.7,  3.2,  3.7,  4.2,  4.7,  5.2,  5.7,  6.2,  6.7,
        7.2,  7.7,  8.2,  8.7,  9.2,  9.7, 10.2, 10.7, 11.2, 11.7, 12.2, 12.7,
       13.2, 13.7, 14.2, 14.7])
Coordinates:
  * range    (range) float64 1.2 1.7 2.2 2.7 3.2 ... 12.7 13.2 13.7 14.2 14.7
Attributes:
    units:    m

2.) Remove data beyond surface level

To reduce the amount of data the code must run through, we can remove all data at and above the water surface. Because the instrument was looking up, we can use the pressure sensor data and the function surface_from_P. This does require that the pressure sensor was ‘zeroed’ prior to deployment. If the instrument is looking down or lacks pressure data, use the function find_surface to detect the seabed or water surface.

ADCPs don’t record water salinity, so it will need to be added to the function. The returned dataset contains the an additional variable “depth”, which is the distance of the water surface away from the seafloor.

After the water surface has been found, data in depth bins at and above the physical water surface can be removed using nan_beyond_surface.

[6]:

dat = api.clean.surface_from_P(dat, salinity=31)
dat = api.clean.nan_beyond_surface(dat)

3.) Correlation filter

Once beyond-surface bins have been removed, ADCP data is typically filtered by correlation clear out spurious velocity datapoints (caused by bubbles, kelp, fish, etc moving through one or multiple beams).

We can take a quick look at the data to see about where this value should be using xarray’s built-in plotting. In the following line of code, we use xarray’s slicing capabilities to show data from beam 1 between a range of 0 to 10 m from the ADCP.

[7]:

%matplotlib inline
dat.corr.sel(beam=1, range=slice(0,10)).plot()

[7]:

<matplotlib.collections.QuadMesh at 0x203ec9024f0>

It’s a good idea to check the other beams as well. Much of this data is high quality, and to not lose data will low correlation caused by natural variation, we’ll use the correlation_filter to set velocity values corresponding to correlations below 50% to NaN.

Note that this threshold is dependent on the deployment environment and instrument, and it isn’t uncommon to use a value as low as 30%, or ~50 counts.

[8]:

dat = api.clean.correlation_filter(dat, thresh=50)

Review the Data

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

ADCPs 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. (Note: the “heading” variable will not change).

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”), 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.

Because this ADCP data was already in the “earth” coordinate system, rotate2 will return the input dataset. set_declination will run correctly no matter the coordinate system.

[9]:

dat = dlfn.set_declination(dat, 15.8) # 15.8 deg East

dat = dlfn.rotate2(dat, 'earth')

Data is already in the earth coordinate system

Because this deployment was set up in ‘burst mode’ (as opposed to an ‘averaging mode’) the next standard step in this analysis is to average the velocity data into time bins. If an instrument was set up to record velocity data in an averaging mode (a specific profiling and averaging interval, e.g. take 5 minutes of data every 30 minutes), this step was completed within the ADCP during deployment and can be skipped.

To average the data into time bins (aka ensembles), start by initiating the binning tool VelBinner. “n_bin” is the number of data points in each ensemble, in this case 300 seconds worth of data, and “fs” is the sampling frequency, which is 1 Hz for this deployment. Once initiated, average the data using the binning tool’s do_avg function.

[10]:

avg_tool = api.VelBinner(n_bin=dat.fs*300, fs=dat.fs)
dat_avg = avg_tool.do_avg(dat)

Two more variables not automatically provided that may be of interest are the horizontal velocity magnitude and its direction, U_mag and U_dir. There are included as “shortcut” functions, and are accessed through the keyword Veldata, as shown in the code block below. The full list of “shorcuts” are listed here, but note all of them might not be available to ADCP data.

Because of the way ADCP coordinate systems were established, when in ENU, U_dir is in “degrees CCW from East” (aka “Cartesian angle”), as opposed to “degrees CW from North”, which is typical for navigation. The convert_degrees function will transform from one to the other and vice versa, and is used below to save U_dir in place.

[11]:

dat_avg['U_mag'] = dat_avg.Veldata.U_mag

dat_avg['U_dir'] = dat_avg.Veldata.U_dir
dat_avg['U_dir'].values = dlfn.tools.misc.convert_degrees(dat_avg['U_dir'].values)

Plotting can be accomplished through the user’s preferred package. Matplotlib is shown here for simplicity, and the horizontal velocity magnitude and direction are plotted below with a blue line delineating the water surface level.

[12]:

%matplotlib inline
from matplotlib import pyplot as plt
import matplotlib.dates as dt

# Convert time to a datetime format
t = dlfn.time.epoch2date(dat_avg.time)

ax = plt.figure(figsize=(12,8)).add_axes([.14, .14, .8, .74])
# Plot the horizontal velocity magnitude
plt.pcolormesh(t, dat_avg.range, dat_avg['U_mag'], cmap='Blues', shading='nearest')
# Plot the water surface
ax.plot(t, dat_avg.depth)

# Set up time on x-axis
ax.set_xlabel('Time')
ax.xaxis.set_major_formatter(dt.DateFormatter('%H:%M'))

ax.set_ylabel('Altitude [m]')
ax.set_ylim([0, 12]);
plt.colorbar(label='Horizontal Vel [m/s]');

[13]:

# Plot horizontal velocity direction
ax = plt.figure(figsize=(12,8)).add_axes([.14, .14, .8, .74])
plt.pcolormesh(t, dat_avg.range, dat_avg['U_dir'], cmap='twilight', shading='nearest')
# Plot the water surface
ax.plot(t, dat_avg.depth)

# set up time on x-axis
ax.set_xlabel('Time')
ax.xaxis.set_major_formatter(dt.DateFormatter('%H:%M'))

ax.set_ylabel('Altitude [m]')
ax.set_ylim([0, 12]);
plt.colorbar(label='Horizontal Vel Dir [deg CW from true N]');

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.

[14]:

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