Comparing Transect Data at the Tanana River Test Site

The TRTS is located near Nenana, Alaska, and was first investigated as a test site in 2008 after a CEC deployment in Eagle Alaska ended early due to debris-related issues. Since then, it has become an established test site with modular deployment infrastructure that can be adapted for different CEC deployments. The TRTS allows CECs to be tested in a harsh environment with cold temperatures, sediment, and heavy debris, representing the conditions of many remote Alaska river communities that could greatly benefit from ME. The following example will familiarize the user with the MHKiT DOLfYN and Delft3D modules by performing a Deflt3d numerical model validation of the Tanana River Test Site (TRTS) against field data.

Start by importing the necessary python packages and MHKiT module.

from matplotlib.pyplot import figure
import scipy.interpolate as interp
import matplotlib.pyplot as plt
import xarray as xr
import pandas as pd
import numpy as np
import matplotlib
import netCDF4
import utm

# MHKiT Imports
from mhkit.dolfyn.rotate import api as ap
from mhkit.dolfyn.adp import api
from mhkit.river.io import d3d
from mhkit import river

1. Preparing the ADCP Field Data

Acoustic Dopple Curent Profiler (ADCP) field data was collected in 2010 on the Tanana river near Nenana. From the picture shown below the flow can be seen enetering the picture from the bottom heading north. The figure further shows the two passes (Transects) of data collected. The data is collected using a barge used with a debris shield infront to protect instruments and test devices. The transect we are looking at are up-steam of the testing barge making them a good repesentiaion of the flow field the Curent Energy Convertes (CEC) are operating in. This transect data was collected using a Teledyne Rio Grand ADCP monted at the Bow of the boat facing downward. The boat faced and moved perpendicular to the river flow across the river as the ADCP was recording.

In the next section we will import the two transects shown here into MHKiT and then interpolate these two passes into a single idealized transect for comparison with the simulated Deflt3D data.

ADCP: Importing Data

The MHKiT DOLfYN api module can import .PDO and other binary formats from ADCP and ADV to an xarray. For this analysis we import two transect passes stored in the ‘data’ folder in .PDO format. For comparison to the numerical model we want to average the two transects into one xarray (‘transect_1_2’).

# Read in the two transect passes
transect_1_raw = api.read(
    "data/river/ADCP_transect/tanana_transects_08_10_10_0_002_10-08-10_142214.PD0"
)
transect_2_raw = api.read(
    "data/river/ADCP_transect/tanana_transects_08_10_10_0_003_10-08-10_143335.PD0"
)
# Create one dataset from the two passes
transect_1_2 = xr.merge([transect_1_raw, transect_2_raw])
# Print the xarray data
transect_1_2

Reading file data/river/ADCP_transect/tanana_transects_08_10_10_0_002_10-08-10_142214.PD0 ...

Reading file data/river/ADCP_transect/tanana_transects_08_10_10_0_003_10-08-10_143335.PD0 ...

<xarray.Dataset> Size: 4MB
Dimensions:              (time: 1227, time_gps: 1227, range: 47, beam: 4,
                          dir: 4, x1: 4, x2: 4, earth: 3, inst: 3)
Coordinates:
  * time                 (time) datetime64[ns] 10kB 2010-08-10T14:28:15.55999...
  * time_gps             (time_gps) datetime64[ns] 10kB 2010-08-10T22:28:17.2...
  * range                (range) float32 188B 0.57 0.82 1.07 ... 11.82 12.07
  * beam                 (beam) int32 16B 1 2 3 4
  * dir                  (dir) <U3 48B 'X' 'Y' 'Z' 'err'
  * x1                   (x1) int64 32B 1 2 3 4
  * x2                   (x2) int64 32B 1 2 3 4
  * earth                (earth) <U1 12B 'E' 'N' 'U'
  * inst                 (inst) <U1 12B 'X' 'Y' 'Z'
Data variables: (12/36)
    number               (time) float64 10kB 3.652e+03 3.653e+03 ... 4.904e+03
    builtin_test_fail    (time) float64 10kB 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    c_sound              (time) float32 5kB 1.466e+03 1.466e+03 ... 1.466e+03
    depth                (time) float32 5kB 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    pitch                (time) float32 5kB -0.1 -0.09 -0.21 ... 0.2 0.07 0.56
    roll                 (time) float32 5kB 3.33 3.32 3.33 ... 2.24 2.65 2.33
    ...                   ...
    rtk_age_gps          (time_gps) float32 5kB 5.0 6.0 2.0 3.0 ... 2.0 3.0 3.0
    speed_over_grnd_gps  (time_gps) float32 5kB 0.03549 0.03395 ... 0.1672
    dir_over_grnd_gps    (time_gps) float32 5kB 115.5 94.08 ... 160.5 182.1
    hdwtime_gps          (time_gps) datetime64[ns] 10kB 2010-08-10T14:28:15.5...
    beam2inst_orientmat  (x1, x2) float32 64B 1.462 -1.462 0.0 ... -1.034 -1.034
    orientmat            (earth, inst, time) float32 44kB -0.9037 ... 1.0
Attributes: (12/38)
    inst_make:                TRDI
    inst_type:                ADCP
    rotate_vars:              ['vel', 'vel_bt']
    has_imu:                  0
    prog_ver:                 10.16
    inst_model:               Rio Grande
    ...                       ...
    false_target_threshold:   50
    transmit_lag_m:           0.08
    bandwidth:                0
    sourceprog:               WINRIVER
    fs:                       1.92
    vel_gps_corrected:        0

xarray.Dataset

• Dimensions:
  • time: 1227
  • time_gps: 1227
  • range: 47
  • beam: 4
  • dir: 4
  • x1: 4
  • x2: 4
  • earth: 3
  • inst: 3
• Coordinates: (9)
  • time
    (time)
    datetime64[ns]
    2010-08-10T14:28:15.559999942 .....

    units :

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

    long_name :

    Time

    standard_name :

    time

    coverage_content_type :

    coordinate

    array(['2010-08-10T14:28:15.559999942', '2010-08-10T14:28:16.319999933',
           '2010-08-10T14:28:17.069999933', ..., '2010-08-10T14:39:59.190000057',
           '2010-08-10T14:39:59.910000085', '2010-08-10T14:40:00.630000114'],
          dtype='datetime64[ns]')

  • time_gps
    (time_gps)
    datetime64[ns]
    2010-08-10T22:28:17.200000047 .....

    units :

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

    long_name :

    Time

    standard_name :

    time

    coverage_content_type :

    coordinate

    array(['2010-08-10T22:28:17.200000047', '2010-08-10T22:28:18.000000000',
           '2010-08-10T22:28:18.599999904', ..., '2010-08-10T22:40:00.799999952',
           '2010-08-10T22:40:01.400000095', '2010-08-10T22:40:02.200000047'],
          dtype='datetime64[ns]')

  • range
    (range)
    float32
    0.57 0.82 1.07 ... 11.82 12.07

    coverage_content_type :

    coordinate

    units :

    m

    long_name :

    Profile Range

    description :

    Distance to the center of each depth bin

    array([ 0.57,  0.82,  1.07,  1.32,  1.57,  1.82,  2.07,  2.32,  2.57,  2.82,
            3.07,  3.32,  3.57,  3.82,  4.07,  4.32,  4.57,  4.82,  5.07,  5.32,
            5.57,  5.82,  6.07,  6.32,  6.57,  6.82,  7.07,  7.32,  7.57,  7.82,
            8.07,  8.32,  8.57,  8.82,  9.07,  9.32,  9.57,  9.82, 10.07, 10.32,
           10.57, 10.82, 11.07, 11.32, 11.57, 11.82, 12.07], dtype=float32)

  • beam
    (beam)
    int32
    1 2 3 4

    units :

    1

    long_name :

    Beam Reference Frame

    coverage_content_type :

    coordinate

    array([1, 2, 3, 4], dtype=int32)

  • dir
    (dir)
    <U3
    'X' 'Y' 'Z' 'err'

    units :

    1

    long_name :

    Reference Frame

    coverage_content_type :

    coordinate

    ref_frame :

    ship

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

  • x1
    (x1)
    int64
    1 2 3 4

    array([1, 2, 3, 4])

  • x2
    (x2)
    int64
    1 2 3 4

    array([1, 2, 3, 4])

  • 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: (36)
  • number
    (time)
    float64
    3.652e+03 3.653e+03 ... 4.904e+03

    units :

    1

    long_name :

    Ensemble Number

    standard_name :

    number_of_observations

    array([3652., 3653., 3654., ..., 4902., 4903., 4904.])

  • builtin_test_fail
    (time)
    float64
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    units :

    1

    long_name :

    Built-In Test Failures

    array([0., 0., 0., ..., 0., 0., 0.])

  • c_sound
    (time)
    float32
    1.466e+03 1.466e+03 ... 1.466e+03

    units :

    m s-1

    long_name :

    Speed of Sound

    standard_name :

    speed_of_sound_in_sea_water

    array([1466., 1466., 1466., ..., 1466., 1466., 1466.], dtype=float32)

  • depth
    (time)
    float32
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    units :

    m

    long_name :

    Depth

    standard_name :

    depth

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

  • pitch
    (time)
    float32
    -0.1 -0.09 -0.21 ... 0.2 0.07 0.56

    units :

    degree

    long_name :

    Pitch

    standard_name :

    platform_pitch

    array([-0.1 , -0.09, -0.21, ...,  0.2 ,  0.07,  0.56], dtype=float32)

  • roll
    (time)
    float32
    3.33 3.32 3.33 ... 2.24 2.65 2.33

    units :

    degree

    long_name :

    Roll

    standard_name :

    platform_roll

    array([3.33, 3.32, 3.33, ..., 2.24, 2.65, 2.33], dtype=float32)

  • heading
    (time)
    float32
    154.6 154.5 155.0 ... 145.8 147.8

    units :

    degree

    long_name :

    Heading

    standard_name :

    platform_orientation

    array([154.65, 154.49, 155.01, ..., 146.32, 145.75, 147.81], dtype=float32)

  • temp
    (time)
    float32
    15.2 15.25 15.23 ... 15.21 15.22

    units :

    degree_C

    long_name :

    Temperature

    standard_name :

    sea_water_temperature

    array([15.2 , 15.25, 15.23, ..., 15.2 , 15.21, 15.22], dtype=float32)

  • salinity
    (time)
    float32
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    units :

    psu

    long_name :

    Salinity

    standard_name :

    sea_water_salinity

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

  • min_preping_wait
    (time)
    float32
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    units :

    s

    long_name :

    Minimum Pre-Ping Wait Time Between Measurements

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

  • heading_std
    (time)
    float32
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    units :

    degree

    long_name :

    Heading Standard Deviation

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

  • pitch_std
    (time)
    float32
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    units :

    degree

    long_name :

    Pitch Standard Deviation

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

  • roll_std
    (time)
    float32
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    units :

    degree

    long_name :

    Roll Standard Deviation

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

  • pressure
    (time)
    float32
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    units :

    dbar

    long_name :

    Pressure

    standard_name :

    sea_water_pressure

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

  • pressure_std
    (time)
    float32
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    units :

    dbar

    long_name :

    Pressure Standard Deviation

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

  • vel
    (dir, range, time)
    float32
    0.057 -0.192 0.121 ... nan nan nan

    units :

    m s-1

    coverage_content_type :

    physicalMeasurement

    long_name :

    Water Velocity

    array([[[ 0.057, -0.192,  0.121, ..., -0.01 , -0.629,  0.134],
            [-0.058, -0.131, -0.119, ...,  0.048, -0.146, -0.078],
            [ 0.059,  0.406,  0.067, ..., -0.123, -0.524,  0.057],
            ...,
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan]],
    
           [[-0.227, -0.047, -0.151, ..., -0.022, -0.488, -0.618],
            [ 0.023, -0.158,  0.204, ..., -0.404, -0.384, -0.396],
            [-0.303, -0.215, -0.019, ..., -0.208, -0.905, -0.434],
            ...,
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan]],
    
           [[-0.01 ,  0.022, -0.119, ...,  0.039,  0.047, -0.032],
            [-0.074, -0.01 , -0.06 , ...,  0.157,  0.089, -0.008],
            [-0.042, -0.013, -0.067, ...,  0.129,  0.067, -0.003],
            ...,
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan]],
    
           [[ 0.26 , -0.043, -0.087, ...,  0.073,  0.187,  0.282],
            [ 0.102,  0.111, -0.204, ...,  0.216,  0.515, -0.105],
            [ 0.324, -0.038,  0.214, ..., -0.08 , -0.142, -0.097],
            ...,
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan]]],
          dtype=float32)

  • amp
    (beam, range, time)
    float32
    224.0 229.0 225.0 ... 63.0 49.0

    units :

    1

    coverage_content_type :

    physicalMeasurement

    long_name :

    Acoustic Signal Amplitude

    standard_name :

    signal_intensity_from_multibeam_acoustic_doppler_velocity_sensor_in_sea_water

    array([[[224., 229., 225., ..., 228., 229., 227.],
            [222., 222., 221., ..., 222., 228., 226.],
            [217., 215., 214., ..., 223., 219., 224.],
            ...,
            [ 53.,  54.,  54., ...,  58.,  61.,  53.],
            [ 52.,  53.,  52., ...,  57.,  58.,  50.],
            [ 54.,  53.,  52., ...,  60.,  56.,  50.]],
    
           [[207., 213., 210., ..., 214., 212., 213.],
            [212., 210., 205., ..., 207., 207., 207.],
            [201., 205., 204., ..., 205., 204., 206.],
            ...,
            [ 46.,  51.,  46., ...,  54.,  63.,  43.],
            [ 43.,  55.,  45., ...,  55.,  53.,  45.],
            [ 45.,  51.,  44., ...,  66.,  49.,  45.]],
    
           [[223., 226., 225., ..., 227., 227., 227.],
            [219., 220., 219., ..., 222., 227., 225.],
            [213., 213., 216., ..., 218., 218., 222.],
            ...,
            [ 49.,  53.,  51., ...,  60.,  60.,  51.],
            [ 52.,  57.,  52., ...,  58.,  62.,  51.],
            [ 50.,  53.,  50., ...,  57.,  78.,  50.]],
    
           [[215., 214., 217., ..., 216., 216., 215.],
            [215., 208., 211., ..., 220., 212., 215.],
            [207., 212., 203., ..., 218., 211., 204.],
            ...,
            [ 51.,  57.,  52., ...,  61.,  92.,  49.],
            [ 52.,  58.,  53., ...,  59.,  66.,  52.],
            [ 51.,  60.,  50., ...,  74.,  63.,  49.]]], dtype=float32)

  • corr
    (beam, range, time)
    float32
    124.0 118.0 130.0 ... 16.0 8.0 21.0

    units :

    1

    coverage_content_type :

    physicalMeasurement

    long_name :

    Acoustic Signal Correlation

    standard_name :

    beam_consistency_indicator_from_multibeam_acoustic_doppler_velocity_profiler_in_sea_water

    array([[[124., 118., 130., ..., 125., 111., 133.],
            [128., 140., 128., ..., 118., 140., 120.],
            [132., 119., 119., ..., 109., 116., 122.],
            ...,
            [ 18.,  17.,  19., ...,  30.,  19.,  36.],
            [ 11.,  21.,  22., ...,  43.,  17.,  17.],
            [ 17.,  11.,  11., ...,  10.,  18.,  15.]],
    
           [[128., 110., 118., ..., 120., 118., 114.],
            [126., 115., 127., ..., 128., 123., 117.],
            [128., 131., 120., ..., 124., 127., 124.],
            ...,
            [ 16.,  14.,  24., ...,  30.,  30.,  20.],
            [ 13.,  11.,  22., ...,   7.,   9.,  21.],
            [ 34.,  18.,  17., ...,  18.,  15.,  16.]],
    
           [[114., 132., 123., ..., 138., 125., 112.],
            [139., 142., 110., ..., 134., 134., 125.],
            [124., 122., 121., ..., 130., 116., 122.],
            ...,
            [ 12.,  16.,  17., ...,  28.,  13.,  14.],
            [ 11.,  17.,  17., ...,  14.,  26.,  13.],
            [ 17.,  17.,  11., ...,  16.,   9.,  15.]],
    
           [[118., 142., 130., ..., 136., 133., 124.],
            [142., 120., 137., ..., 111., 119., 130.],
            [110., 121., 126., ..., 146., 132., 114.],
            ...,
            [ 19.,  13.,  17., ...,  18.,  23.,  11.],
            [  7.,  25.,  16., ...,  26.,  23.,  27.],
            [ 13.,  25.,  13., ...,  16.,   8.,  21.]]], dtype=float32)

  • prcnt_gd
    (beam, range, time)
    float32
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    units :

    %

    coverage_content_type :

    physicalMeasurement

    long_name :

    Percent Good

    standard_name :

    proportion_of_acceptable_signal_returns_from_acoustic_instrument_in_sea_water

    array([[[  0.,   0.,   0., ...,   0.,   0.,   0.],
            [  0.,   0.,   0., ...,   0.,   0.,   0.],
            [  0.,   0.,   0., ...,   0.,   0.,   0.],
            ...,
            [  0.,   0.,   0., ...,   0.,   0.,   0.],
            [  0.,   0.,   0., ...,   0.,   0.,   0.],
            [  0.,   0.,   0., ...,   0.,   0.,   0.]],
    
           [[  0.,   0.,   0., ...,   0.,   0.,   0.],
            [  0.,   0.,   0., ...,   0.,   0.,   0.],
            [  0.,   0.,   0., ...,   0.,   0.,   0.],
            ...,
            [  0.,   0.,   0., ...,   0.,   0.,   0.],
            [  0.,   0.,   0., ...,   0.,   0.,   0.],
            [  0.,   0.,   0., ...,   0.,   0.,   0.]],
    
           [[  0.,   0.,   0., ...,   0.,   0.,   0.],
            [  0.,   0.,   0., ...,   0.,   0.,   0.],
            [  0.,   0.,   0., ...,   0.,   0.,   0.],
            ...,
            [100., 100., 100., ..., 100., 100., 100.],
            [100., 100., 100., ..., 100., 100., 100.],
            [100., 100., 100., ..., 100., 100., 100.]],
    
           [[100., 100., 100., ..., 100., 100., 100.],
            [100., 100., 100., ..., 100., 100., 100.],
            [100., 100., 100., ..., 100., 100., 100.],
            ...,
            [  0.,   0.,   0., ...,   0.,   0.,   0.],
            [  0.,   0.,   0., ...,   0.,   0.,   0.],
            [  0.,   0.,   0., ...,   0.,   0.,   0.]]], dtype=float32)

  • dist_bt
    (beam, time)
    float32
    3.95 3.95 3.96 ... 1.93 1.95 2.0

    units :

    m

    coverage_content_type :

    physicalMeasurement

    long_name :

    Bottom Track Measured Depth

    array([[3.95, 3.95, 3.96, ..., 2.65, 2.64, 2.6 ],
           [2.55, 2.56, 2.57, ..., 1.61, 1.59, 1.59],
           [3.31, 3.32, 3.33, ..., 1.93, 1.95, 1.96],
           [2.87, 2.88, 2.88, ..., 1.93, 1.95, 2.  ]], dtype=float32)

  • vel_bt
    (dir, time)
    float32
    0.063 -0.004 0.039 ... 0.01 0.0

    units :

    m s-1

    coverage_content_type :

    physicalMeasurement

    long_name :

    Water Velocity

    array([[ 6.3000e-02, -4.0000e-03,  3.9000e-02, ...,  2.3000e-02,
             2.7000e-02,  0.0000e+00],
           [-3.7000e-02,  1.9000e-02,  1.5000e-02, ..., -6.9000e-02,
            -1.0100e-01, -1.5800e-01],
           [ 9.0000e-03,  0.0000e+00,  8.0000e-03, ...,  0.0000e+00,
            -1.7000e-02, -2.0000e-03],
           [-1.0000e-02, -6.0000e-03, -3.2768e+01, ...,  7.0000e-03,
             1.0000e-02,  0.0000e+00]], dtype=float32)

  • corr_bt
    (beam, time)
    float32
    231.0 243.0 242.0 ... 249.0 246.0

    units :

    1

    coverage_content_type :

    physicalMeasurement

    long_name :

    Acoustic Signal Correlation

    standard_name :

    beam_consistency_indicator_from_multibeam_acoustic_doppler_velocity_profiler_in_sea_water

    array([[231., 243., 242., ..., 252., 255., 248.],
           [229., 239., 246., ..., 249., 248., 254.],
           [223., 237.,   0., ..., 249., 252., 251.],
           [248., 244., 250., ..., 248., 249., 246.]], dtype=float32)

  • amp_bt
    (beam, time)
    float32
    40.0 40.0 36.0 ... 40.0 64.0 49.0

    units :

    1

    coverage_content_type :

    physicalMeasurement

    long_name :

    Acoustic Signal Amplitude

    standard_name :

    signal_intensity_from_multibeam_acoustic_doppler_velocity_sensor_in_sea_water

    array([[40., 40., 36., ..., 64., 66., 51.],
           [45., 38., 47., ..., 78., 78., 83.],
           [40., 43., 36., ..., 55., 71., 71.],
           [55., 54., 52., ..., 40., 64., 49.]], dtype=float32)

  • prcnt_gd_bt
    (beam, time)
    float32
    0.0 0.0 100.0 ... 100.0 100.0 100.0

    units :

    %

    coverage_content_type :

    physicalMeasurement

    long_name :

    Percent Good

    standard_name :

    proportion_of_acceptable_signal_returns_from_acoustic_instrument_in_sea_water

    array([[  0.,   0., 100., ...,   0.,   0.,   0.],
           [  0.,   0.,   0., ...,   0.,   0.,   0.],
           [  0.,   0.,   0., ...,   0.,   0.,   0.],
           [100., 100.,   0., ..., 100., 100., 100.]], dtype=float32)

  • latitude_gps
    (time_gps)
    float32
    64.56 64.56 64.56 ... 64.56 64.56

    units :

    degrees_north

    long_name :

    Latitude

    standard_name :

    latitude

    array([64.560905, 64.560905, 64.560905, ..., 64.56094 , 64.56094 ,
           64.56094 ], dtype=float32)

  • longitude_gps
    (time_gps)
    float32
    -149.1 -149.1 ... -149.1 -149.1

    units :

    degrees_east

    long_name :

    Longitude

    standard_name :

    longitude

    array([-149.06679, -149.06679, -149.06679, ..., -149.06682, -149.06682,
           -149.0668 ], dtype=float32)

  • fix_gps
    (time_gps)
    float32
    9.0 9.0 9.0 9.0 ... 9.0 9.0 9.0 9.0

    units :

    1

    long_name :

    GPS Fix

    array([9., 9., 9., ..., 9., 9., 9.], dtype=float32)

  • n_sat_gps
    (time_gps)
    float32
    9.0 9.0 9.0 9.0 ... 10.0 10.0 10.0

    units :

    count

    long_name :

    Number of Satellites

    array([ 9.,  9.,  9., ..., 10., 10., 10.], dtype=float32)

  • hdop_gps
    (time_gps)
    float32
    0.9 0.9 0.9 0.9 ... 0.8 0.8 0.8 0.8

    units :

    1

    long_name :

    Horizontal Dilution of Precision

    array([0.9, 0.9, 0.9, ..., 0.8, 0.8, 0.8], dtype=float32)

  • elevation_gps
    (time_gps)
    float32
    108.5 108.2 108.2 ... 108.6 108.6

    units :

    m

    long_name :

    Elevation above MLLW

    array([108.483, 108.196, 108.208, ..., 108.644, 108.629, 108.606],
          dtype=float32)

  • rtk_age_gps
    (time_gps)
    float32
    5.0 6.0 2.0 3.0 ... 6.0 2.0 3.0 3.0

    units :

    s

    long_name :

    Age of Received Real Time Kinetic Signal

    array([5., 6., 2., ..., 2., 3., 3.], dtype=float32)

  • speed_over_grnd_gps
    (time_gps)
    float32
    0.03549 0.03395 ... 0.1096 0.1672

    units :

    m s-1

    long_name :

    Platform Speed over Ground

    standard_name :

    platform_speed_wrt_ground

    array([0.03549382, 0.03395062, 0.03858025, ..., 0.11985596, 0.1095679 ,
           0.16718106], dtype=float32)

  • dir_over_grnd_gps
    (time_gps)
    float32
    115.5 94.08 6.931 ... 160.5 182.1

    units :

    degree

    long_name :

    Platform Direction over Ground

    standard_name :

    platform_course

    array([115.501,  94.077,   6.931, ..., 140.963, 160.468, 182.087],
          dtype=float32)

  • hdwtime_gps
    (time_gps)
    datetime64[ns]
    2010-08-10T14:28:15.559999942 .....

    array(['2010-08-10T14:28:15.559999942', '2010-08-10T14:28:16.319999933',
           '2010-08-10T14:28:17.069999933', ...,
           '2010-08-10T14:39:59.190000057', '2010-08-10T14:39:59.910000085',
           '2010-08-10T14:40:00.630000114'], dtype='datetime64[ns]')

  • beam2inst_orientmat
    (x1, x2)
    float32
    1.462 -1.462 0.0 ... -1.034 -1.034

    units :

    1

    long_name :

    Rotation Matrix

    array([[ 1.4619021 , -1.4619021 ,  0.        ,  0.        ],
           [ 0.        ,  0.        , -1.4619021 ,  1.4619021 ],
           [ 0.26604444,  0.26604444,  0.26604444,  0.26604444],
           [ 1.033721  ,  1.033721  , -1.033721  , -1.033721  ]],
          dtype=float32)

  • orientmat
    (earth, inst, time)
    float32
    -0.9037 -0.9025 -0.9064 ... 1.0 1.0

    units :

    1

    long_name :

    Orientation Matrix

    array([[[-0.9037092 , -0.9025101 , -0.9063815 , ..., -0.8321477 ,
             -0.8265897 , -0.84628606],
            [-0.4281468 , -0.43066862, -0.42246014, ..., -0.554554  ,
             -0.56280506, -0.53272873],
            [-0.        , -0.        , -0.        , ..., -0.        ,
             -0.        , -0.        ]],
    
           [[ 0.4281468 ,  0.43066862,  0.42246014, ...,  0.554554  ,
              0.56280506,  0.53272873],
            [-0.9037092 , -0.9025101 , -0.9063815 , ..., -0.8321477 ,
             -0.8265897 , -0.84628606],
            [ 0.        ,  0.        ,  0.        , ...,  0.        ,
              0.        ,  0.        ]],
    
           [[-0.        , -0.        , -0.        , ..., -0.        ,
             -0.        , -0.        ],
            [ 0.        ,  0.        ,  0.        , ...,  0.        ,
              0.        ,  0.        ],
            [ 1.        ,  1.        ,  1.        , ...,  1.        ,
              1.        ,  1.        ]]], dtype=float32)

• Indexes: (9)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2010-08-10 14:28:15.559999942',
                   '2010-08-10 14:28:16.319999933',
                   '2010-08-10 14:28:17.069999933',
                   '2010-08-10 14:28:17.799999952',
                   '2010-08-10 14:28:18.529999971',
                   '2010-08-10 14:28:19.279999971',
                   '2010-08-10 14:28:20.029999971',
                   '2010-08-10 14:28:20.779999971',
                   '2010-08-10 14:28:21.539999961',
                   '2010-08-10 14:28:22.299999952',
                   ...
                   '2010-08-10 14:39:54.119999885',
                   '2010-08-10 14:39:54.849999904',
                   '2010-08-10 14:39:55.569999933',
                   '2010-08-10 14:39:56.289999961',
                   '2010-08-10 14:39:57.019999980',
                      '2010-08-10 14:39:57.750000',
                   '2010-08-10 14:39:58.480000019',
                   '2010-08-10 14:39:59.190000057',
                   '2010-08-10 14:39:59.910000085',
                   '2010-08-10 14:40:00.630000114'],
                  dtype='datetime64[ns]', name='time', length=1227, freq=None))

  • time_gps
    PandasIndex

    PandasIndex(DatetimeIndex(['2010-08-10 22:28:17.200000047',
                             '2010-08-10 22:28:18',
                   '2010-08-10 22:28:18.599999904',
                   '2010-08-10 22:28:19.400000095',
                   '2010-08-10 22:28:20.200000047',
                   '2010-08-10 22:28:20.799999952',
                   '2010-08-10 22:28:21.599999904',
                   '2010-08-10 22:28:22.400000095',
                   '2010-08-10 22:28:23.200000047',
                   '2010-08-10 22:28:23.799999952',
                   ...
                   '2010-08-10 22:39:55.599999904',
                   '2010-08-10 22:39:56.400000095',
                   '2010-08-10 22:39:57.200000047',
                   '2010-08-10 22:39:57.799999952',
                   '2010-08-10 22:39:58.599999904',
                   '2010-08-10 22:39:59.400000095',
                             '2010-08-10 22:40:00',
                   '2010-08-10 22:40:00.799999952',
                   '2010-08-10 22:40:01.400000095',
                   '2010-08-10 22:40:02.200000047'],
                  dtype='datetime64[ns]', name='time_gps', length=1227, freq=None))

  • range
    PandasIndex

    PandasIndex(Index([0.5699999928474426, 0.8199999928474426, 1.0700000524520874,
           1.3200000524520874, 1.5700000524520874, 1.8200000524520874,
            2.069999933242798,  2.319999933242798,  2.569999933242798,
            2.819999933242798,  3.069999933242798,  3.319999933242798,
            3.569999933242798,  3.819999933242798,  4.070000171661377,
            4.320000171661377,  4.570000171661377,  4.820000171661377,
            5.070000171661377,  5.320000171661377,  5.570000171661377,
            5.820000171661377,  6.070000171661377,  6.320000171661377,
            6.570000171661377,  6.820000171661377,  7.070000171661377,
            7.320000171661377,  7.570000171661377,  7.820000171661377,
            8.069999694824219,  8.319999694824219,  8.569999694824219,
            8.819999694824219,  9.069999694824219,  9.319999694824219,
            9.569999694824219,  9.819999694824219, 10.069999694824219,
           10.319999694824219, 10.569999694824219, 10.819999694824219,
           11.069999694824219, 11.319999694824219, 11.569999694824219,
           11.819999694824219, 12.069999694824219],
          dtype='float32', name='range'))

  • beam
    PandasIndex

    PandasIndex(Index([1, 2, 3, 4], dtype='int32', name='beam'))

  • dir
    PandasIndex

    PandasIndex(Index(['X', 'Y', 'Z', 'err'], dtype='object', name='dir'))

  • x1
    PandasIndex

    PandasIndex(Index([1, 2, 3, 4], dtype='int64', name='x1'))

  • x2
    PandasIndex

    PandasIndex(Index([1, 2, 3, 4], dtype='int64', name='x2'))

  • earth
    PandasIndex

    PandasIndex(Index(['E', 'N', 'U'], dtype='object', name='earth'))

  • inst
    PandasIndex

    PandasIndex(Index(['X', 'Y', 'Z'], dtype='object', name='inst'))

• Attributes: (38)

  inst_make :

  TRDI

  inst_type :

  ADCP

  rotate_vars :

  ['vel', 'vel_bt']

  has_imu :

  0

  prog_ver :

  10.16

  inst_model :

  Rio Grande

  beam_angle :

  20

  freq :

  1200

  beam_pattern :

  convex

  orientation :

  down

  n_beams :

  4

  n_cells :

  47

  pings_per_ensemble :

  1

  cell_size :

  0.25

  blank_dist :

  0.25

  profiling_mode :

  12

  min_corr_threshold :

  64

  n_code_reps :

  4

  min_prcnt_gd :

  0

  max_error_vel :

  1.5

  sec_between_ping_groups :

  0.09

  coord_sys :

  ship

  use_pitchroll :

  yes

  use_3beam :

  yes

  bin_mapping :

  yes

  heading_misalign_deg :

  0.0

  magnetic_var_deg :

  0.0

  sensors_src :

  01011101

  sensors_avail :

  00011101

  bin1_dist_m :

  0.57

  transmit_pulse_m :

  0.3

  water_ref_cells :

  [0, 5]

  false_target_threshold :

  50

  transmit_lag_m :

  0.08

  bandwidth :

  0

  sourceprog :

  WINRIVER

  fs :

  1.92

  vel_gps_corrected :

  0

ADCP: Convert the Coordinate System to UTM

By default, the ADCP coordinates are stored as latitude and longitude coordinates. However, latitude and longitude cannot be directly converted to a distance therefor here we convert from latitude and longitude to \(UTM_x\) and \(UTM_y\). This gives us the ability to calculate the linear line of best fit for the idealized transect.

# Convert Coordiantes to UTM using utm module
utm_x_y = utm.from_latlon(transect_1_2.latitude_gps, transect_1_2.longitude_gps, 6, "W")

# Create a DataFrame from the points
gps = [[x, y] for x, y in zip(utm_x_y[0], utm_x_y[1])]
gps_points = pd.DataFrame(np.array(gps), columns=["utm_x", "utm_y"])

ADCP: Declination

The magnetic declination cause data to be rotated slightly this can be corrected using MHKiT DOLfYN’s rotate api set_declination. If magnetic declination has already been accounted for with the deployment software, it will be saved under the attributes magnetic_var_deg for Teledyne instruments. The variable magnetic_var_deg is unaltered by the function set_declination as it creates a new variable (declination).

# Nenana Alaska is 15.7 deg East
angle = 15.7
ap.set_declination(transect_1_2, angle, inplace=True)

ADCP: Rotation of Data Coordinate System

Based on the deployment setting the coordinate system is set to “beam”,”inst”,”earth”, or “principal”. Using the rotate2 function the coordinate system can by changed to fit the needs of the analysis. The inplace set as true will alter the input dataset “in place”, eg. it does not create a new dataset. This is further explained in the river_example ipython notebook.

# Rotate to 'earth' coordinate system
api.rotate2(transect_1_2, "earth", inplace=True)

ADCP: Ideal Transect

To perform a Delft3D transect comparison with the ADCP data an “ideal transect”, defined as a straight path perpendicular to the flow that the boat was aiming for, was processed from the ADCP field data. Using the combined transect_1_2 dataset a first degree numpy polyfit function (linear regression) was used to find a slope and intercept. Using the equation we can calculate \(UTM_y\) coordinates to create an ideal linear transect. The resultant ideal transect is shown in the figure below.

# Linear regression using first order polyfit
a, b = np.polyfit(gps_points.utm_x, gps_points.utm_y, 1)

# Generate a DataFrame of points from the linear regression
ideal = [[x, y] for x, y in zip(gps_points.utm_x, a * gps_points.utm_x + b)]
ideal_points = pd.DataFrame(np.array(ideal), columns=["utm_x", "utm_y"])

# Repeat UTM corrdinates to match the ADCP points matrix (dir, range, time)
utm_x_points = np.tile(gps_points.utm_x, np.size(transect_1_2.range))
utm_y_points = np.tile(a * gps_points.utm_x + b, np.size(transect_1_2.range))
depth_points = np.repeat(transect_1_2.range, np.size(gps_points.utm_x))

ADCP_ideal_points = {
    "utm_x": utm_x_points,
    "utm_y": utm_y_points,
    "waterdepth": depth_points,
}
ADCP_ideal_points = pd.DataFrame(ADCP_ideal_points)

# Initialize the figure
figure(figsize=(8, 6))
fig, ax = plt.subplots()

# Get data from the original transect in UTM for comparison
transect_1 = utm.from_latlon(
    transect_1_raw.latitude_gps, transect_1_raw.longitude_gps, 6, "W"
)
transect_2 = utm.from_latlon(
    transect_2_raw.latitude_gps, transect_2_raw.longitude_gps, 6, "W"
)

# Plot the original transect data for comparison
plt.plot(transect_1[0], transect_1[1], "b", label="GPS Transect 1")
plt.plot(transect_2[0], transect_2[1], "r--", label="GPS Transect 2")

# Plot the Idealized Transect
plt.plot(
    ADCP_ideal_points.utm_x, ADCP_ideal_points.utm_y, "k-.", label="Ideal Transect"
)
plt.ticklabel_format(style="scientific", useOffset=False)
ax.get_xaxis().set_major_formatter(
    matplotlib.ticker.FuncFormatter(lambda x, p: format(int(x), ","))
)
plt.xticks(rotation=45)

# Plot Settings
plt.legend()
plt.xlabel("$UTM_x (m)$")
plt.ylabel("$UTM_y (m)$")

Text(0, 0.5, '$UTM_y (m)$')

<Figure size 800x600 with 0 Axes>

ADCP: Range Offset

Depending on where the ADCP is deployed the range values might need to be offset by the deployment depth. For the TRTS dataset used here the range was corrected during the deployment setup, therefore, a zero-input used. This step is provided as a reference for other users.

# Adjust the range offset, included here for reference
offset = 0
api.clean.set_range_offset(transect_1_2, offset)

ADCP: Correlation Filter

Correlation is a measure of data quality, and its output is scaled in units such that the expected correlation (given high signal/noise ratio, S/N) is 128. Here data was fileted for a minimum correlation of 40 counts. [https://www.comm-tec.com/Docs/Manuali/RDI/BBPRIME.pdf]

Looking at the plot in the figure below the ADCP is facing downward from the surface of the river and is defined as a range of 0. This means the top of the river is at the bottom of the plot. This ADCP has 4 beams the plot shown is one the 4 beams. From the plot of one of the beams shown below we can see that some of the values below the river bottom are filtered out.

# Apply the correlation filter
min_correlation = 40
transect_1_2 = api.clean.correlation_filter(transect_1_2, thresh=min_correlation)

# Plot the results the (data is displayed upside-down)
transect_1_2.corr.sel(beam=1).plot()

<matplotlib.collections.QuadMesh at 0x1c42c7a4e90>

ADCP: River Bottom Filter

MHKiT DOLfYN has the functionality to find the river bottom based on the along-beam acoustic amplitude data recorded from the ADCP. However, if depth sounder data is available it can be more reliable at representing the river bottom. Here we use the depth sounder (depth_sounder) to create an array that can be used as a filter when multiplied with the velocity data.

# Filtering out depth sounder values above the river surface
depth_sounder = transect_1_2.where(transect_1_2.dist_bt > 0)

# Of the 4 values beams get the shallowest depth value at each location
bottom = np.min(depth_sounder.dist_bt, axis=0)

# River bottom for ideal transect
bottom_avg = interp.griddata(gps_points, bottom, ideal_points, method="linear")

# Create a matrix of depths
bottom_filter = d3d.create_points(
    x=bottom_avg, y=transect_1_2.range.to_numpy(), waterdepth=1
)

# Creating a mask matrix with ones in the area of the river cross section and nan's outside
river_bottom_filter = []
for index, row in bottom_filter.iterrows():
    if row["x"] > row["y"]:
        filter = 1
    else:
        filter = float("nan")
    river_bottom_filter = np.append(river_bottom_filter, filter)

ADCP: Get Original Transect Data

Here we create a DataFrame of the original velocity components from the transect data. This DataFrame will be used to interpolate the velocity data onto the ideal transect.

# Tiling the GPS data for each depth bin
gps_utm_x = np.tile(gps_points.utm_x, np.size(transect_1_2.range))
gps_utm_y = np.tile(gps_points.utm_y, np.size(transect_1_2.range))

# Repeating the depth bins for each GPS point
depth = np.repeat(transect_1_2.range, np.size(gps_points.utm_x))

# Create Dataframe from the calculated points
ADCP_points = pd.DataFrame(
    {"utm_x": gps_utm_x, "utm_y": gps_utm_y, "waterdepth": depth}
)

# Raveling the veocity data to correspond with 'ADCP_points' and filtering out velocity data bellow the river bottom
ADCP_points["east_velocity"] = np.ravel(transect_1_2.vel[0, :, :]) * river_bottom_filter
ADCP_points["north_velocity"] = (
    np.ravel(transect_1_2.vel[1, :, :]) * river_bottom_filter
)
ADCP_points["vertical_velocity"] = (
    np.ravel(transect_1_2.vel[2, :, :]) * river_bottom_filter
)
ADCP_points = ADCP_points.dropna()

# Show points
ADCP_points

	utm_x	utm_y	waterdepth	east_velocity	north_velocity	vertical_velocity
0	400953.784848	7.161133e+06	0.57	-0.094245	0.214233	-0.010
1	400953.784848	7.161133e+06	0.57	0.181185	0.079026	0.022
2	400953.784848	7.161133e+06	0.57	-0.143789	0.129486	-0.119
3	400953.784848	7.161133e+06	0.57	-0.432847	0.519561	-0.095
4	400953.784848	7.161133e+06	0.57	-0.258352	0.393507	-0.106
...	...	...	...	...	...	...
37885	400989.041693	7.161150e+06	8.07	0.074346	1.819609	0.084
39093	400993.719096	7.161153e+06	8.32	0.409135	0.774590	0.055
39094	400993.719096	7.161153e+06	8.32	0.230414	0.237873	0.038
39098	400992.389563	7.161152e+06	8.32	0.381595	0.560519	0.104
39101	400991.707623	7.161151e+06	8.32	0.058718	0.515059	0.152

22018 rows × 6 columns

ADCP: Interpolating Velocities for Ideal Transect

Using the original transect data saved in ADCP_points the velocity values are interpolated onto the ideal transect points using the scipy interpolate griddata function.

# Remove duplicate points
ADCP_points = ADCP_points.drop_duplicates(subset=["utm_x", "utm_y", "waterdepth"])

# Normalize the data to avoid precision issues
min_coords = ADCP_points[["utm_x", "utm_y", "waterdepth"]].min()
max_coords = ADCP_points[["utm_x", "utm_y", "waterdepth"]].max()
ADCP_points[["utm_x", "utm_y", "waterdepth"]] = (ADCP_points[["utm_x", "utm_y", "waterdepth"]] - min_coords) / (max_coords - min_coords)

ADCP_ideal_points[["utm_x", "utm_y", "waterdepth"]] = (ADCP_ideal_points[["utm_x", "utm_y", "waterdepth"]] - min_coords) / (max_coords - min_coords)

C:\Users\sterl\AppData\Local\Temp\ipykernel_4716\1528100513.py:7: SettingWithCopyWarning:
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  ADCP_points[["utm_x", "utm_y", "waterdepth"]] = (ADCP_points[["utm_x", "utm_y", "waterdepth"]] - min_coords) / (max_coords - min_coords)

# Project velocity onto ideal tansect
ADCP_ideal = pd.DataFrame()
ADCP_ideal["east_velocity"] = interp.griddata(
    ADCP_points[["utm_x", "utm_y", "waterdepth"]],
    ADCP_points["east_velocity"],
    ADCP_ideal_points[["utm_x", "utm_y", "waterdepth"]],
    method="linear",
    rescale=True
)
ADCP_ideal["north_velocity"] = interp.griddata(
    ADCP_points[["utm_x", "utm_y", "waterdepth"]],
    ADCP_points["north_velocity"],
    ADCP_ideal_points[["utm_x", "utm_y", "waterdepth"]],
    method="linear",
    rescale=True
)
ADCP_ideal["vertical_velocity"] = interp.griddata(
    ADCP_points[["utm_x", "utm_y", "waterdepth"]],
    ADCP_points["vertical_velocity"],
    ADCP_ideal_points[["utm_x", "utm_y", "waterdepth"]],
    method="linear",
    rescale=True
)

# Calculate the magnitude of the velocity components
ADCP_ideal["magnitude"] = np.sqrt(
    ADCP_ideal.east_velocity**2
    + ADCP_ideal.north_velocity**2
    + ADCP_ideal.vertical_velocity**2
)

# Denormalize the ideal points for plotting
ADCP_points[["utm_x", "utm_y", "waterdepth"]] = (ADCP_points[["utm_x", "utm_y", "waterdepth"]] *  (max_coords - min_coords) + min_coords)
ADCP_ideal_points[["utm_x", "utm_y", "waterdepth"]] = ADCP_ideal_points[["utm_x", "utm_y", "waterdepth"]] * (max_coords - min_coords) + min_coords

C:\Users\sterl\AppData\Local\Temp\ipykernel_4716\1835236727.py:33: SettingWithCopyWarning:
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  ADCP_points[["utm_x", "utm_y", "waterdepth"]] = (ADCP_points[["utm_x", "utm_y", "waterdepth"]] *  (max_coords - min_coords) + min_coords)

Using the calculated magnitude of velocity from the ADCP we can create a contour plot of the velocity magnitude.

# Set the contour color bar bounds
min_plot = 0
max_plot = 3

# The Contour of velocity magnitude from the ADCP transect data
fig, ax = plt.subplots(figsize=(10, 4.4))

contour_plot = plt.tripcolor(
    ADCP_ideal_points.utm_x,
    -ADCP_ideal_points.waterdepth,
    ADCP_ideal.magnitude * river_bottom_filter,
    vmin=min_plot,
    vmax=max_plot,
)

plt.xlabel("$UTM_x (m)$")
plt.ylabel("Water Depth (m)")
cbar = plt.colorbar(contour_plot)
cbar.set_label("velocity [m/s]")
plt.ylim([-8.5, -1])
plt.xlim([400950, 401090])
plt.plot(ideal_points.utm_x, -bottom_avg, "k", label="river bottom")
plt.legend(loc=7)
ax.get_xaxis().set_major_formatter(
    matplotlib.ticker.FuncFormatter(lambda x, p: format(int(x), ","))
)
plt.xticks(rotation=45)

(array([400940., 400960., 400980., 401000., 401020., 401040., 401060.,
        401080., 401100.]),
 [Text(400940.0, 0, '400,940'),
  Text(400960.0, 0, '400,960'),
  Text(400980.0, 0, '400,980'),
  Text(401000.0, 0, '401,000'),
  Text(401020.0, 0, '401,020'),
  Text(401040.0, 0, '401,040'),
  Text(401060.0, 0, '401,060'),
  Text(401080.0, 0, '401,080'),
  Text(401100.0, 0, '401,100')])

ADCP: Ideal Transect Down-sampled

Numerical models such as Delft3D use discretized points to represent the domain. Environmental models tend to have a low number of points to represent large continuous areas. For instance, the bird’s eye view of Delft3D only calculated numerical effects at a roughly 5m by 5m square. Delft3D then discretizes the depth with an even percentage from the water level to the bed by 12 layers as specified for this initial case. As the collected ADCP data had far more resolution than the numerical model a down-sampling was applied using an average to match the locations at which the Delft3D model returned data. This roughly equated to ten velocity measurements spanning the width of the river and 12 measurments spanning the depth.

# Interpolate points by getting min & max first
start_utmx = min(ADCP_ideal_points.utm_x)
start_utmy = min(ADCP_ideal_points.utm_y)

end_utmx = max(ADCP_ideal_points.utm_x)
end_utmy = min(ADCP_ideal_points.utm_y)

# Using N points for x calculate the y values on an ideal transect from the linear regression used earlier
N = 10
utm_x_ideal_downsampeled = np.linspace(start_utmx, end_utmx, N)
utm_y_ideal_downsampeled = (a * utm_x_ideal_downsampeled) + b


# Plot the Idealized Transect for comparison
plt.plot(
    ADCP_ideal_points.utm_x, ADCP_ideal_points.utm_y, ".", ms=1, label="Ideal Transect"
)

# Plot the downsampled transect
plt.plot(
    utm_x_ideal_downsampeled,
    utm_y_ideal_downsampeled,
    "ro",
    label="Down Sampled Ideal Transect",
)


# Plot settings
ax.get_xaxis().set_major_formatter(
    matplotlib.ticker.FuncFormatter(lambda x, p: format(int(x), ","))
)
plt.xticks(rotation=45)
plt.xlabel("$UTM_x$")
plt.ylabel("$UTM_y$")
plt.legend()

<matplotlib.legend.Legend at 0x1c445525450>

With our down-sampled x,y points we can now create a down-sampled z at our x,y locations. As shown in the figure red dots below the river bottom will be excluded from the interpolation.

# Create an idealized depth N layers deep
N_layers = 12
downsampled_depth = np.linspace(
    np.nanmin(transect_1_2.range), np.nanmax(bottom_avg), N_layers
)

# Repeat this over the N points of the DownSampled Ideal Transect above
depth_ideal_points_downsampled = np.repeat(downsampled_depth, N)

# Tile the x, y over the N of layers to add to a DataFrame
utm_x_ideal_points_downsampled = np.tile(utm_x_ideal_downsampeled, N_layers)
utm_y_ideal_points_downsampled = np.tile(utm_y_ideal_downsampeled, N_layers)

# Create a Dataframe of our idealized x,y,depth points
ADCP_ideal_points_downsamples = pd.DataFrame(
    {
        "utm_x": utm_x_ideal_points_downsampled,
        "utm_y": utm_y_ideal_points_downsampled,
        "waterdepth": depth_ideal_points_downsampled,
    }
)

# Plot the Down sampled data points at the x locations
fig, ax = plt.subplots(figsize=(10, 4.4))
plt.plot(
    ADCP_ideal_points_downsamples.utm_x,
    ADCP_ideal_points_downsamples.waterdepth * -1,
    "ro",
)

# Plot the ADCP river bed
plt.plot(ideal_points.utm_x, -bottom_avg, "k", label="river bottom")

# Plot settings
ax.get_xaxis().set_major_formatter(
    matplotlib.ticker.FuncFormatter(lambda x, p: format(int(x), ","))
)
plt.xticks(rotation=45)
plt.title("DownSampled Ideal Transect Depth")
plt.xlabel("$UTM_x [m]$")
plt.ylabel("$ Depth [m]$")

Text(0, 0.5, '$ Depth [m]$')

ADCP: Interpolating the Downsampled Velocites for the Ideal Transect

The above process of creating a river_bottom_filter and interpolating onto the ideal transect is repeated for the down-sampled grid shown in the figure above. The depth sounder data is saved above in the variable 'bottom' and used again here.

# Normalize the data to avoid precision issues
min_coords = ADCP_points[["utm_x", "utm_y", "waterdepth"]].min()
max_coords = ADCP_points[["utm_x", "utm_y", "waterdepth"]].max()
ADCP_points[["utm_x", "utm_y", "waterdepth"]] = (ADCP_points[["utm_x", "utm_y", "waterdepth"]] - min_coords) / (max_coords - min_coords)

ADCP_ideal_points_downsamples[["utm_x", "utm_y", "waterdepth"]] = (ADCP_ideal_points_downsamples[["utm_x", "utm_y", "waterdepth"]] - min_coords) / (max_coords - min_coords)

C:\Users\sterl\AppData\Local\Temp\ipykernel_4716\3546222366.py:4: SettingWithCopyWarning:
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  ADCP_points[["utm_x", "utm_y", "waterdepth"]] = (ADCP_points[["utm_x", "utm_y", "waterdepth"]] - min_coords) / (max_coords - min_coords)

# Project velocity onto ideal tansect
ADCP_ideal_downsamples = pd.DataFrame()
ADCP_ideal_downsamples["east_velocity"] = interp.griddata(
    ADCP_points[["utm_x", "utm_y", "waterdepth"]],
    ADCP_points["east_velocity"],
    ADCP_ideal_points_downsamples[["utm_x", "utm_y", "waterdepth"]],
    method="linear",
)
ADCP_ideal_downsamples["north_velocity"] = interp.griddata(
    ADCP_points[["utm_x", "utm_y", "waterdepth"]],
    ADCP_points["north_velocity"],
    ADCP_ideal_points_downsamples[["utm_x", "utm_y", "waterdepth"]],
    method="linear",
)
ADCP_ideal_downsamples["vertical_velocity"] = interp.griddata(
    ADCP_points[["utm_x", "utm_y", "waterdepth"]],
    ADCP_points["vertical_velocity"],
    ADCP_ideal_points_downsamples[["utm_x", "utm_y", "waterdepth"]],
    method="linear",
)
ADCP_ideal_downsamples["magnitude"] = np.sqrt(
    ADCP_ideal_downsamples.east_velocity**2
    + ADCP_ideal_downsamples.north_velocity**2
    + ADCP_ideal_downsamples.vertical_velocity**2
)
ADCP_ideal_downsamples

# Denormalize the ideal points for plotting
ADCP_points[["utm_x", "utm_y", "waterdepth"]] = (ADCP_points[["utm_x", "utm_y", "waterdepth"]] *  (max_coords - min_coords) + min_coords)
ADCP_ideal_points_downsamples[["utm_x", "utm_y", "waterdepth"]] = ADCP_ideal_points_downsamples[["utm_x", "utm_y", "waterdepth"]] * (max_coords - min_coords) + min_coords

C:\Users\sterl\AppData\Local\Temp\ipykernel_4716\14157162.py:29: SettingWithCopyWarning:
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  ADCP_points[["utm_x", "utm_y", "waterdepth"]] = (ADCP_points[["utm_x", "utm_y", "waterdepth"]] *  (max_coords - min_coords) + min_coords)

ADCP: Plot Down-Sampled Ideal Transect

To plot the contour of the ideal sample points we need to first create a mask to remove data outside of the river bottom because the ideal transect grid (shown as red dots in the figure above) includes points below the river bottom.

# Create a DataFrame of downsampled points
ideal_downsampeled = [
    [x, y] for x, y in zip(utm_x_ideal_downsampeled, utm_y_ideal_downsampeled)
]
ideal_points_downsampled = pd.DataFrame(
    np.array(ideal_downsampeled), columns=["utm_x", "utm_y"]
)

# River bottom for downsampled ideal transect
bottom_avg_downsampled = interp.griddata(
    gps_points, bottom, ideal_points_downsampled, method="linear"
)

# Create a matrix of depths
bottom_filter_downsampled = d3d.create_points(
    x=bottom_avg_downsampled, y=downsampled_depth, waterdepth=1
)

# Creating a mask matrix with ones in the area of the river cross section and nan's outside
river_bottom_filter_downsampled = []
for index, row in bottom_filter_downsampled.iterrows():
    if row["x"] > row["y"]:
        filter = 1
    else:
        filter = float("nan")
    river_bottom_filter_downsampled = np.append(river_bottom_filter_downsampled, filter)

Using our down-sampled velocity magnitude a contour plot is created. For comparison the river bottom from the original ideal transect is shown.

# Plotting
fig, ax = plt.subplots(figsize=(10, 4.4))
contour_plot = plt.tripcolor(
    ADCP_ideal_points_downsamples.utm_x,
    -ADCP_ideal_points_downsamples.waterdepth,
    ADCP_ideal_downsamples.magnitude * river_bottom_filter_downsampled,
    vmin=min_plot,
    vmax=max_plot,
)

# Plot river bottom for comparison
plt.plot(ideal_points.utm_x, -bottom_avg, "k", label="river bottom")

# Plot Settings
plt.xlabel("$UTM_x$ (m)")
plt.ylabel("Water Depth (m)")
cbar = plt.colorbar(contour_plot)
cbar.set_label("Velocity [m/s]")
plt.ylim([-8.5, -1])
plt.xlim([400950, 401090])
plt.legend(loc=7)
ax.get_xaxis().set_major_formatter(
    matplotlib.ticker.FuncFormatter(lambda x, p: format(int(x), ","))
)
plt.xticks(rotation=45)

(array([400940., 400960., 400980., 401000., 401020., 401040., 401060.,
        401080., 401100.]),
 [Text(400940.0, 0, '400,940'),
  Text(400960.0, 0, '400,960'),
  Text(400980.0, 0, '400,980'),
  Text(401000.0, 0, '401,000'),
  Text(401020.0, 0, '401,020'),
  Text(401040.0, 0, '401,040'),
  Text(401060.0, 0, '401,060'),
  Text(401080.0, 0, '401,080'),
  Text(401100.0, 0, '401,100')])

## 2. Delft3D Data

Typically researchers use a modeling tool such as Deflt3D which, relative to field deployments, can perform low-cost simulations where changes to the environment can be implemented and interpreted quickly. This allows researchers to narrow in on a limited design set for field investigation thereby increasing the efficiency of field deployments. To initiate the numerical model, a 2010 bathymetry survey taken by Terrasond under contract to the University of Alaska Fairbanks (UAF) was used to represent the domain. This domain was discretized using the Delft3D built-in mesh generation tool, RGF grid, to roughly 5m x 5m unstructured squares from a bird’s eye view. Delft3D represents depth using shallow water equations with a sigma layer approach which distributes the layers between the specified bottom of the river and the water level. An initial 12 sigma layers were specified to compare to filed data for this study. The domain was simulated in Delft3D at a constant discharge of :math:`1789

m^{3}/s` and a \(2.7 m\) water level, as recorded by the USGS survey station for August \(10^{th}\) 2010. The simulation shown in the example is designed to match the experimental conditions of when the ADCP plotted above was collected. The Delft3D model takes inputs of water level, discharge, and the river bathymetry. The United States Geological Survey (USGS) has a data collection downstream of the TRTS (station:”15515500”). The water level and discharge are taken from the USGS data,

and the bathymetry data was collected in is the summer of 2010. In this example we will show how to get data from the USGS station and then jump to analyzing the output from the Delft3D model.

D3D: USGS Water Level

The United States Geological Survey has a nearby survey station (#15515500) downstream of the TRTS which records discharge and water level. The MHKiT River module has the capability to request data from USGS survey stations as demonstrated in the river example. Using MHKiT’s river module the water level for the USGS station was collected for the same day as the ADCP data was collected. The USGS assigns a parameter number of ‘00065’ for water level or Gage height. The water level was averaged to be roughly 2.7m. This number is directly input into the Delft3D model and isn’t used in this example notebook. Further examples of using the river module are shown in the river_example. Similar to the water level the discharge was found using the same MHKit river module except swapping out the parameter number for ‘00060’ to get the discharge data. The discharge was averaged to be 1789 \(m^3/s\). This value was also directly input to the Delft3D model and isn’t used in this example notebook.

# Use the requests method to obtain 1 day of instantneous gage height data
water_level_USGS_data = river.io.usgs.request_usgs_data(
    station="15515500",
    parameter="00065",
    start_date="2010-08-10",
    end_date="2010-08-10",
    data_type="Instantaneous",
)

# Plot data
water_level_USGS_data.plot()

# Plot Settings
plt.xlabel("Time")
plt.ylabel("Gage Height (feet)")

# Use the requests method to obtain 1 day of instantneous discharge data
discharge_USGS_data = river.io.usgs.request_usgs_data(
    station="15515500",
    parameter="00060",
    start_date="2010-08-10",
    end_date="2010-08-10",
    data_type="Instantaneous",
)

# Print data
discharge_USGS_data.plot()
# Plot Settings
plt.xlabel("Time")
plt.ylabel("Dischage ($f^3/s$)")

Text(0, 0.5, 'Dischage ($f^3/s$)')

D3D: Importing Delft3D Data

Here we are importing the NetCDF data exported by the Delft3D simulation described in the Delft3D introduction.

# Import the simulated data
d3d_data = netCDF4.Dataset("data/river/ADCP_transect/tanana81010_final_map.nc")

# Get the ADCP sample points
ADCP_ideal_points_downsamples_xy = ADCP_ideal_points_downsamples.rename(
    columns={"utm_x": "x", "utm_y": "y"}
)

D3D: Ideal Transect Data

Next, we get the velocity data at the same points as the down-sampled ideal ADCP data. In Delft3D velocity in the y direction is ucy, velocity in the x direction is ucx and velocity in the z direction is ucz. From these velocity components the velocity magnitude is calculated using the root mean squared formula.

# Interpolate the Delft3D simulated data onto the the sample points
variables = ["ucy", "ucx", "ucz"]
D3D = d3d.variable_interpolation(
    d3d_data, variables, points=ADCP_ideal_points_downsamples_xy
)

# Calculate the magnitude of the velocity
D3D["magnitude"] = np.sqrt(D3D.ucy**2 + D3D.ucx**2 + D3D.ucz**2)

points provided

Using our Delft3D velocity magnitude a contour plot is created. The results are plotted with the original ideal ADCP river bottom for reference.

# Plot Delft3D interpolated Data
fig, ax = plt.subplots(figsize=(10, 4.4))
contour_plot = plt.tripcolor(
    D3D.x,
    -D3D.waterdepth,
    D3D.magnitude * river_bottom_filter_downsampled,
    vmin=min_plot,
    vmax=max_plot,
    # shading='gouraud'
    alpha=1,
)

# Plot the river bottom calculated frol ADCP for comparison
plt.plot(ideal_points.utm_x, -bottom_avg, "k", label="river bottom")

# Figure settings
plt.xlabel("$UTM_x (m)$")
plt.ylabel("Water Depth (m)")
cbar = plt.colorbar(contour_plot)
cbar.set_label("velocity [m/s]")
plt.ylim([-8.5, -1])
plt.xlim([400960, 401090])
plt.legend(loc=7)
ax.get_xaxis().set_major_formatter(
    matplotlib.ticker.FuncFormatter(lambda x, p: format(int(x), ","))
)
plt.xticks(rotation=45)

(array([400960., 400980., 401000., 401020., 401040., 401060., 401080.,
        401100.]),
 [Text(400960.0, 0, '400,960'),
  Text(400980.0, 0, '400,980'),
  Text(401000.0, 0, '401,000'),
  Text(401020.0, 0, '401,020'),
  Text(401040.0, 0, '401,040'),
  Text(401060.0, 0, '401,060'),
  Text(401080.0, 0, '401,080'),
  Text(401100.0, 0, '401,100')])

3. Error Calulations

Data from Delft3D and the ADCP were now on an equivalent set of points in our data processing environment. For comparison the \(L_1\), \(L_2\), and \(L_{\infty}\) norms were calculated for the Delft3D transect, and the down-sampled ADCP data. The L1 norm is equal to the absolute value of the difference between modeled and measured velocities, Delft3D and ADCP, respectively, over the measured velocity. This gives a fraction that can be converted into the percent difference between measured and modeled velocities. The mean absolute error (MAE) in is an average value of L1. The error calculation process was repeated to calculate the L2 norm which is commonly referred to as the mean squared error (MSE). \(L_1\), \(L_2\), and \(L_{\infty}\) are dimensionless.

Error: L1

\(L_1= \frac{|D3D-ADCP|}{ADCP}\)

\(L_1\) is a dimensionless value equal to the difference between the modeled and measured value divided by the actual value. In this case the measured ADCP is considered to be the actual value.

# L1
L1_Magnitude = (
    abs(ADCP_ideal_downsamples.magnitude - D3D.magnitude)
    / ADCP_ideal_downsamples.magnitude
)

Similar to the river_bottom_filter this filter replaces extreme values (100%) values with ‘nan’ and is then combined with the river_bottom_filter to create error_filter.

river_bottom_edge_filter_downsampled = []
for i in L1_Magnitude:
    if 1 > i:
        filter = 1
    else:
        filter = float("nan")
    river_bottom_edge_filter_downsampled = np.append(
        river_bottom_edge_filter_downsampled, filter
    )

error_filter = river_bottom_edge_filter_downsampled * river_bottom_filter_downsampled

The mean absolute error (MAE) is the average \(L_1\) error over the transect area. The MAE is displayed with the fil contour plot of \(L_1\). \(N\) is equation is equle to the number of points in the river cross-sectional area.

\(MAE= \frac{\sum{L_1}}{N}\)

# Calculate and priont the Mean Absolute Error
MAE = np.sum(L1_Magnitude * error_filter) / len(L1_Magnitude[L1_Magnitude < 1000])
MAE

np.float64(0.33569684683914963)

# Set the min and max error values
max_plot_error = 1
min_plot_error = 0

# Plotting the L1 error
fig, ax = plt.subplots(figsize=(10, 4.4))
contour_plot_L1 = plt.tripcolor(
    D3D.x,
    -D3D.waterdepth,
    L1_Magnitude * error_filter,
    vmin=min_plot_error,
    vmax=max_plot_error,
)

# Plot the river bottom for comparison
plt.plot(ideal_points.utm_x, -bottom_avg, "k", label="river bottom")

# Plot settings
plt.xlim([400960, 401090])
plt.ylim([-8.5, -1])
plt.xlabel("UTM x (m)")
plt.ylabel("Water Depth (m)")
cbar = plt.colorbar(contour_plot_L1)
cbar.set_label("$L_1$ Velocity Error")
plt.legend(loc=7)

ax.get_xaxis().set_major_formatter(
    matplotlib.ticker.FuncFormatter(lambda x, p: format(int(x), ","))
)
plt.xticks(rotation=45)

(array([400960., 400980., 401000., 401020., 401040., 401060., 401080.,
        401100.]),
 [Text(400960.0, 0, '400,960'),
  Text(400980.0, 0, '400,980'),
  Text(401000.0, 0, '401,000'),
  Text(401020.0, 0, '401,020'),
  Text(401040.0, 0, '401,040'),
  Text(401060.0, 0, '401,060'),
  Text(401080.0, 0, '401,080'),
  Text(401100.0, 0, '401,100')])

Error: L2

\(L_2=\left(\frac{D3D-ADCP}{ADCP}\right)^2\)

\(L_2\) is similar to \(L_1\) except \(L_2\) squares the difference between the modeled and measured values.

# L2
L2_Magnitude = (
    (ADCP_ideal_downsamples.magnitude - D3D.magnitude)
    / ADCP_ideal_downsamples.magnitude
) ** 2

The mean squared error (MSE) is the average of \(L_2\) error over the transect area. \(N\) is equation is equle to the number of points in the river cross-sectional area.

\(MSE= \frac{\sum{L_2}}{N}\)

MSE = np.sum(L2_Magnitude * error_filter) / np.size(L2_Magnitude[L2_Magnitude < 1000])
MSE

np.float64(0.15431169343691561)

# Create a contour plot of the error
# Plotting
fig, ax = plt.subplots(figsize=(10, 4.4))
contour_plot_L2 = plt.tripcolor(
    D3D.x,
    -D3D.waterdepth,
    L2_Magnitude * error_filter,
    vmin=min_plot_error,
    vmax=max_plot_error,
)

# Plot the river bottom for comparison
plt.plot(ideal_points.utm_x, -bottom_avg, "k", label="river bottom")

# Plot settings
plt.xlim([400960, 401090])
plt.ylim([-8.5, -1])
plt.xlabel("UTM x (m)")
plt.ylabel("Water Depth (m)")
cbar = plt.colorbar(contour_plot_L1)
cbar.set_label("$L_2$ Velocity Error")
plt.legend(loc=7)

ax.get_xaxis().set_major_formatter(
    matplotlib.ticker.FuncFormatter(lambda x, p: format(int(x), ","))
)
plt.xticks(rotation=45)

C:\Users\sterl\AppData\Local\Temp\ipykernel_4716\3877048929.py:20: UserWarning: Adding colorbar to a different Figure <Figure size 1000x440 with 3 Axes> than <Figure size 1000x440 with 1 Axes> which fig.colorbar is called on.
  cbar = plt.colorbar(contour_plot_L1)

(array([400960., 400980., 401000., 401020., 401040., 401060., 401080.,
        401100.]),
 [Text(400960.0, 0, '400,960'),
  Text(400980.0, 0, '400,980'),
  Text(401000.0, 0, '401,000'),
  Text(401020.0, 0, '401,020'),
  Text(401040.0, 0, '401,040'),
  Text(401060.0, 0, '401,060'),
  Text(401080.0, 0, '401,080'),
  Text(401100.0, 0, '401,100')])

Error: \(L_\infty\)

\(L_\infty = max(L_1)\)

\(L_\infty\) is the maximum L1 error

# L inf
L_inf = np.nanmax(L1_Magnitude * error_filter)
L_inf

np.float64(0.888520413692259)

4. Discussion

A model verification using simulated data of TRTS was compared to field data collected in 2010 using MHKiT. This example showcased the functionality of the MHKiT river and DOLfYN modules using a real-world study. The findings demonstrated how MHKiT modules may be combined to analyze simulated data and collected field data in a simple and quick manner. Other transect plotting software run separately from data analysis software. MHKiT’s ability to process data from multiple common marine energy sources will help ME practitioners reach unique and novel results quickly by eliminating much of the start-up processing scripts most analysts must perform and maintain today. Although this is a river example, this framework of analysing modeled data and comparing it to wave and tidal experimental data is possible as well.

The high error values found between the model and field data were expected for a first attempt on a coarse numerical model domain. The coarse model domain allows MHKiT to host the data in this repository for public use and replication of this study’s results. This example further allows users to swap in their own data and build on this working example. Using this template the TRTS Delft3D model may now be tuned and easily compared to the error metrics calculated to improve our Delft3D Tanana River model. Looking forward, a primary focus of model improvement would be the resolution of the Delft3D model, especially on the TRTS section of the Tanana river. Improving the Delft3D model is of particular interest to the TRTS staff as it can be used to predict the impact of the wake of a CEC on the surrounding environment such as the river bottom before it’s deployed long-term. This can save time and funding when implementing CECs onto an electrical grid while also mitigating environmental impact.